Цикл и от 1 до н
п1 = п1 + п2*п3/п4
Метод Гауса
Вобщем меня интересует, кто реализовал на С++ метод Гаусса. И сколько требуется времени для того чтоб свести матрицу long double *dArrey[1000*1000]; к диагональной. У меня это получается за 5.273 с. Я не знаю сколько это. Много или мало? Что скажете товарищи программисты? :)
Забыл дописать. Скорость то зависит от машины. ЦП-1700Гц
Количество арифметических операций в методе Гаусса примерно равно (2*n^2)/3. В твоём случае их будет примерно 1 333 333. Возьмём каждую операцию, скажем, за шесть тактов процессора и увидим, что даже при таком грубом преуменьшении возможностей современных компьютеров, твоя задача должна быть решена менее чем за секунду. Полагай...
Если в методе Гаусса нужо выбрать ведущий элемент - то операций больше. Пример во вложении. На моей машине (C-333 160 Мб ОЗУ, VC6 Debug) это заняло более 3 минут
Я вообще то пошел по такому принципу делим первый елемент матрици второй строки на первый елемент первой строки. Тогда умножаем всю вторую строку на этот елемнт и отнимаем первую строку. И так дальше по диагонали. Тоесть чем ближе к последней строке тем менше оперераций умножения и вычитания и на конец в последней строке вообще их две.
for(UINT iMeter0X=0; iMeter0X<iSize0X; iMeter0X++)
{
uiSizeY_iMeter0X=uiSizeY*iMeter0X;
for(UINT iMeterY=iMeter0X+1;iMeterY<uiSizeY;iMeterY++)
{
uiSizeY_iMeterY=uiSizeY*iMeterY;
dMultiplauer=m_dArrey[uiSizeY_iMeter0X+iMeter0X]/m_dArrey[uiSizeY_iMeterY+iMeter0X];
for(UINT iMeterXX=iMeter0X; iMeterXX<uiSizeY; iMeterXX++)
{
m_dArrey[uiSizeY_iMeterY+iMeterXX]=
m_dArrey[uiSizeY_iMeterY+iMeterXX]*dMultiplauer-m_dArrey[uiSizeY_iMeter0X+iMeterXX];
}
m_dArreyB[iMeterY]=m_dArreyB[iMeterY]-m_dArreyB[iMeter0X]*dMultiplauer;
}
}
for(UINT iMeter0X=0; iMeter0X<iSize0X; iMeter0X++)
{
uiSizeY_iMeter0X=uiSizeY*iMeter0X;
for(UINT iMeterY=iMeter0X+1;iMeterY<uiSizeY;iMeterY++)
{
uiSizeY_iMeterY=uiSizeY*iMeterY;
dMultiplauer=m_dArrey[uiSizeY_iMeter0X+iMeter0X]/m_dArrey[uiSizeY_iMeterY+iMeter0X];
for(UINT iMeterXX=iMeter0X; iMeterXX<uiSizeY; iMeterXX++)
{
m_dArrey[uiSizeY_iMeterY+iMeterXX]=
m_dArrey[uiSizeY_iMeterY+iMeterXX]*dMultiplauer-m_dArrey[uiSizeY_iMeter0X+iMeterXX];
}
m_dArreyB[iMeterY]=m_dArreyB[iMeterY]-m_dArreyB[iMeter0X]*dMultiplauer;
}
}
Это верный принцип. В схеме с выбором ведущего элемента на каждом шаге i (i=1,2,3,..n n-размер матрицы) предварительно находится наибольший (по абсолютной величине) элемент по всей матрице начиная со строки i. Этот элемент путем перестановок строк и столбцов ставится на первое место строки i.
В твоем же случае это просто всегда первый элемент первой строки, что может сказаться на точности
В твоем же случае это просто всегда первый элемент первой строки, что может сказаться на точности
Это я заметил, максимальная точность была 99.99999997%. Просто если искать максимальный элемент по абсолютной величине тогда много сравнений а насколько я помню сравнение занимает у процессора около 4-х тактов, это выше чем, допустим, вычитание и добавление.
Обычно в мат.алгоритмах поступают наоборот - идут в ущерб скорости и повышают качество (точность). Иначе нужно обоснование, что заданная точность допустима и т.п. ... К тому же оптимизировать можно и другие шаги алгоритма.
в твоем случае перестановки и поиск ведущего элемента не займут много времени - посмотри вложение в предыдущем моем посте.Это программа, которая вначале генерирует заведомо невырожденную (не сингулярную) квадратную матрицу A (детерминант отличет от нуля, иначе не решить систему), затем считает матрицу B (при условии что все x1=x2=...=xn=1). После этого находятся корни уравнения
Ax=B в процедуре asg2d (кстати взятой с сервера численного анализа МГУ)
На моем 333 Celeron выполнение программы решения методом Гаусса заняло около 3-4 минут при условии что размер матрицы A равен 1000*1000. На более мощном компьютере скорость возрастет.
в твоем случае перестановки и поиск ведущего элемента не займут много времени - посмотри вложение в предыдущем моем посте.Это программа, которая вначале генерирует заведомо невырожденную (не сингулярную) квадратную матрицу A (детерминант отличет от нуля, иначе не решить систему), затем считает матрицу B (при условии что все x1=x2=...=xn=1). После этого находятся корни уравнения
Ax=B в процедуре asg2d (кстати взятой с сервера численного анализа МГУ)
На моем 333 Celeron выполнение программы решения методом Гаусса заняло около 3-4 минут при условии что размер матрицы A равен 1000*1000. На более мощном компьютере скорость возрастет.
В принципи, я решил вставить один и другой алгоритм в программу, и попробивать сравнить результаты. Чтото реализую поделюсь успехом.
Вобщем написал я код, немного долго исправлял мат ошибки, но теперь все работает и етот код сводит матрицу к диагональной за 8.8с. Но и точность н очень то заметна.
void CMetodGausaDlg::OnCalculateMaxEl()
{
UINT uiTimeCalcMaxEl=GetTickCount();
long double dMaxEl=fabs(m_dArrey[0]);
long double dMultiplauer=0;
int iSaveIndex=0;
int iSizeY=0;
int iSizeX=0;
int buf=0;
/*заполняем масив индексов, чтоб не делать физической перестановки елементов*.
for(int i=0; i<m_DialogSize.m_EDIT_SIZE1; i++)
{
m_iIndex=i;
}
//Цикл которий двигаетса по столбцам
for(int x=0; x<m_DialogSize.m_EDIT_SIZE1-1; x++)
{
//Присваиваем сменной dMaxEl первый елемент сравнения dMaxEl=fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[x]+x]);
//В цикле который перемещается по строкам
for(int j=x+1; j<m_DialogSize.m_EDIT_SIZE1; j++)
{
//Находим максимальний елемент if(dMaxEl<fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[j]+x]))
{
MaxEl=fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[j]+x]);
iSaveIndex=j;
}
}
Меняем индекси строк
buf=m_iIndex[x];
m_iIndex[x]=m_iIndex[iSaveIndex];
m_iIndex[iSaveIndex]=buf;
iSizeX=m_DialogSize.m_EDIT_SIZE1*m_iIndex[x];
//Делим строку на максимальний елемент
dMaxEl=m_dArrey[iSizeX+x];
for(int c=x; c<m_DialogSize.m_EDIT_SIZE1; c++)
{
m_dArrey[iSizeX+c]=m_dArrey[iSizeX+c]/dMaxEl;
}
m_dArreyB[m_iIndex[x]]=m_dArreyB[m_iIndex[x]]/dMaxEl;
/*Умножаем строку с главним елементом на множитель и вычитаем ёё из текущей*/
for(int z=x+1; z<m_DialogSize.m_EDIT_SIZE1; z++)
{
iSizeY=m_DialogSize.m_EDIT_SIZE1*m_iIndex[z];
dMultiplauer=m_dArrey[iSizeY+x];
for(int k=0; k<m_DialogSize.m_EDIT_SIZE1; k++)
{
m_dArrey[iSizeY+k]=m_dArrey[iSizeY+k]-m_dArrey[iSizeX+k]*dMultiplauer;
}
m_dArreyB[m_iIndex[z]]=m_dArreyB[m_iIndex[z]]-m_dArreyB[m_iIndex[x]]*dMultiplauer;
}
}
}
void CMetodGausaDlg::OnCalculateMaxEl()
{
UINT uiTimeCalcMaxEl=GetTickCount();
long double dMaxEl=fabs(m_dArrey[0]);
long double dMultiplauer=0;
int iSaveIndex=0;
int iSizeY=0;
int iSizeX=0;
int buf=0;
/*заполняем масив индексов, чтоб не делать физической перестановки елементов*.
for(int i=0; i<m_DialogSize.m_EDIT_SIZE1; i++)
{
m_iIndex=i;
}
//Цикл которий двигаетса по столбцам
for(int x=0; x<m_DialogSize.m_EDIT_SIZE1-1; x++)
{
//Присваиваем сменной dMaxEl первый елемент сравнения dMaxEl=fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[x]+x]);
//В цикле который перемещается по строкам
for(int j=x+1; j<m_DialogSize.m_EDIT_SIZE1; j++)
{
//Находим максимальний елемент if(dMaxEl<fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[j]+x]))
{
MaxEl=fabs(m_dArrey[m_DialogSize.m_EDIT_SIZE1*m_iIndex[j]+x]);
iSaveIndex=j;
}
}
Меняем индекси строк
buf=m_iIndex[x];
m_iIndex[x]=m_iIndex[iSaveIndex];
m_iIndex[iSaveIndex]=buf;
iSizeX=m_DialogSize.m_EDIT_SIZE1*m_iIndex[x];
//Делим строку на максимальний елемент
dMaxEl=m_dArrey[iSizeX+x];
for(int c=x; c<m_DialogSize.m_EDIT_SIZE1; c++)
{
m_dArrey[iSizeX+c]=m_dArrey[iSizeX+c]/dMaxEl;
}
m_dArreyB[m_iIndex[x]]=m_dArreyB[m_iIndex[x]]/dMaxEl;
/*Умножаем строку с главним елементом на множитель и вычитаем ёё из текущей*/
for(int z=x+1; z<m_DialogSize.m_EDIT_SIZE1; z++)
{
iSizeY=m_DialogSize.m_EDIT_SIZE1*m_iIndex[z];
dMultiplauer=m_dArrey[iSizeY+x];
for(int k=0; k<m_DialogSize.m_EDIT_SIZE1; k++)
{
m_dArrey[iSizeY+k]=m_dArrey[iSizeY+k]-m_dArrey[iSizeX+k]*dMultiplauer;
}
m_dArreyB[m_iIndex[z]]=m_dArreyB[m_iIndex[z]]-m_dArreyB[m_iIndex[x]]*dMultiplauer;
}
}
}
Цитата: AVDEY
Вобщем написал я код, немного долго исправлял мат ошибки, но теперь все работает и етот код сводит матрицу к диагональной за 8.8с. Но и точность н очень то заметна.
Заметь, что деление на элемент равносильно умножению на 1.0/значение.
Умножение быстрее деления.
Кроме того, если используешь метод Гаусса с выбором ведущего элемента максимальный по модулю элемент в матрице нужно искать каждый раз после сведения очередного максимального элемента к единице. По-моему, ты этого не делаешь.
К тому же стоит ли держать массив индексов, когда проще (и понятней) поменять местами строки и столбцы?
по части реализации по-моему метод Гаусса лучше бы смотрелся в отдельной процедуре с независимыми параметрами, чем в методе класса диалогового окна...
Цитата: michael_is_98
Кроме того, если используешь метод Гаусса с выбором ведущего элемента максимальный по модулю элемент в матрице нужно искать каждый раз после сведения очередного максимального элемента к единице. По-моему, ты этого не делаешь.
Ну как же не делаю.
Первый цикл передвигаетса по столбцам. =0
Тоисть находим максимальний елемент в первом столбце.
Меняем индекси.
Делим строку с главным элементом на главный элемент.
Умножаем строку с главным элементом на множитель и вычитаем с текущей .
Потом возвращаемся к первому циклу =1.
Находим максимальний элемент. (Проверка производитса без первой строки так как етот елемент будет менше других).
Извини, по диагонали просмотрел код. Как ты запоминаешь перестановки столбцов, строк? Ведь для этого нужно завести массив перестановок, который должен быть равен 2*n (n - размер матрицы). Если перестановок не было, в массив поставить текущие значения строки и столбца.
Лучше бы в отдельную процедуру вынести, имена переменных сократить.
Лучше бы в отдельную процедуру вынести, имена переменных сократить.
Зачем 2*n достаточно масива Arrey[n] n-размер матрици. Так как переставляем строки, тоисть [Y]. Тогда елмент b[0][1] что b[0][n] будет иметь тот же индекс.
Создаем масив индексов.
UINT uiArreyIndex[n];
Инициализируем елементи масива по умолчанию тоисть
uiArreyIndex[0]=0;
uiArreyIndex[n]=n;
for(int i=0; i<n; i++)
{
uiArreyIndex=i;
}
Тогда доступ к масиву будет через индекс например первий елемент третей строки
b[uiArreyIndex[2]][0];
когда надо поменять строки ми просто меняем значения масива индексов
Например надо поменять первую и третью строки.
uiArreyIndex[0]=2;
uiArreyIndex[2]=0;
И тогда физически строка 2 остаетса на своем же месте но доступ к ней уже будет:
b[uiArreyIndex[0]][0];
так как uiArreyIndex[0] у нас уже 2.
Я вот не понял на счет того что умножение быстрее деления но 1.0/значение ето и есть деление?
Создаем масив индексов.
UINT uiArreyIndex[n];
Инициализируем елементи масива по умолчанию тоисть
uiArreyIndex[0]=0;
uiArreyIndex[n]=n;
for(int i=0; i<n; i++)
{
uiArreyIndex=i;
}
Тогда доступ к масиву будет через индекс например первий елемент третей строки
b[uiArreyIndex[2]][0];
когда надо поменять строки ми просто меняем значения масива индексов
Например надо поменять первую и третью строки.
uiArreyIndex[0]=2;
uiArreyIndex[2]=0;
И тогда физически строка 2 остаетса на своем же месте но доступ к ней уже будет:
b[uiArreyIndex[0]][0];
так как uiArreyIndex[0] у нас уже 2.
Я вот не понял на счет того что умножение быстрее деления но 1.0/значение ето и есть деление?
Насчет второй части здесь все просто - если в цикле нужно разделить на какое-то значение (знач.), то перед входом в цикл лучше получить значение знач2=1./знач., а затем в теле цикла умножать на него (знач2). Где-то у тебя это можно использовать.
По поводу перестановок - ты делаешь лишние телодвижения. Зачем тебе инициализировать массив перестановок, это можно делать в процессе вычилений - это раз. Второе, самое главное - переставлять (переставлять физически) лучше не только строки, но и столбцы. В этом случае у тебя на диагонали матрицы всегда будут стоять 1, а по значению массива размерности 2*n ты всегда узнаешь, какая перестановка столбца и строки на каждом шаге была сделана.
По поводу перестановок - ты делаешь лишние телодвижения. Зачем тебе инициализировать массив перестановок, это можно делать в процессе вычилений - это раз. Второе, самое главное - переставлять (переставлять физически) лучше не только строки, но и столбцы. В этом случае у тебя на диагонали матрицы всегда будут стоять 1, а по значению массива размерности 2*n ты всегда узнаешь, какая перестановка столбца и строки на каждом шаге была сделана.
Но сколько займет времени перестановка двух строк матрицы размером 1000х1000, это надо 3000 присваиваний или 2000 прмсваиваний и 2000 сложений. А здесь на перестановку 1000 элементов уходит всего 3 пере присваивания.
А на счет 1/значение, упадет точность вычислений лишняя операция с возможной потерей младшего разряда.
А на счет 1/значение, упадет точность вычислений лишняя операция с возможной потерей младшего разряда.
Цитата: AVDEY
Но сколько займет времени перестановка двух строк матрицы размером 1000х1000, это надо 3000 присваиваний или 2000 прмсваиваний и 2000 сложений. А здесь на перестановку 1000 элементов уходит всего 3 пере присваивания.
тебе главное - решить систему методом ГаусСа. А чтобы ее решить, нужно чтобы на диагонали стояли единичные элементы. Будешь переставлять строки и столбцы - сможешь этого достичь. Если будешь переставлять строки - не уверен. А сколько это займет времени - не так важно по сравнению с логичной структурой алгоритма. Так , по-моему, хотя тебе решать ...
Цитата:
А на счет 1/значение, упадет точность вычислений лишняя операция с возможной потерей младшего разряда.
Лишняя операция (деление) будет накапливаться в цикле. К тому же деление выполняется всегда медленнее, чем умножение.
Т.е.
Код:
медленнее чем
Код:
знач = п3/п4
Цикл и от 1 до н
п1 = п1 + п2*знач
Цикл и от 1 до н
п1 = п1 + п2*знач
Ну даже если переставлять и строки и столбцы то нужно 2 массива UINT uiIndex_X[n]; и UINT uiIndex_Y[n]; так как мы переставляем весь столбец и всю строку а не отдельный элемент.
Но у меня в коде нет лишнего деления.
Так как у меня отдельный цикл на деление.
И потом умножение идет по строках а не по столбцах, и даже если по столбцах то деление должно быть одно на весь столбец. А все остальные операции это умножение с вычитанием, так как после деления первый элемент равен 1 тогда делитель находится в строке где должен быть 0 но с противоположным знаком.
dMaxEl=m_dArrey[iSizeX+x];
for(int c=x; c<m_DialogSize.m_EDIT_SIZE1; c++)
{
m_dArrey[iSizeX+c]=m_dArrey[iSizeX+c]/dMaxEl;
}
m_dArreyB[m_iIndex[x]]=m_dArreyB[m_iIndex[x]]/dMaxEl;
/*Умножаем строку с главним елементом на множитель и вычитаем ёё из текущей*/
for(int z=x+1; z<m_DialogSize.m_EDIT_SIZE1; z++)
{
iSizeY=m_DialogSize.m_EDIT_SIZE1*m_iIndex[z];
dMultiplauer=m_dArrey[iSizeY+x];
for(int k=0; k<m_DialogSize.m_EDIT_SIZE1; k++)
{
m_dArrey[iSizeY+k]=m_dArrey[iSizeY+k]-m_dArrey[iSizeX+k]*dMultiplauer;
}
m_dArreyB[m_iIndex[z]]=m_dArreyB[m_iIndex[z]]-m_dArreyB[m_iIndex[x]]*dMultiplauer;
}
Но у меня в коде нет лишнего деления.
Так как у меня отдельный цикл на деление.
И потом умножение идет по строках а не по столбцах, и даже если по столбцах то деление должно быть одно на весь столбец. А все остальные операции это умножение с вычитанием, так как после деления первый элемент равен 1 тогда делитель находится в строке где должен быть 0 но с противоположным знаком.
dMaxEl=m_dArrey[iSizeX+x];
for(int c=x; c<m_DialogSize.m_EDIT_SIZE1; c++)
{
m_dArrey[iSizeX+c]=m_dArrey[iSizeX+c]/dMaxEl;
}
m_dArreyB[m_iIndex[x]]=m_dArreyB[m_iIndex[x]]/dMaxEl;
/*Умножаем строку с главним елементом на множитель и вычитаем ёё из текущей*/
for(int z=x+1; z<m_DialogSize.m_EDIT_SIZE1; z++)
{
iSizeY=m_DialogSize.m_EDIT_SIZE1*m_iIndex[z];
dMultiplauer=m_dArrey[iSizeY+x];
for(int k=0; k<m_DialogSize.m_EDIT_SIZE1; k++)
{
m_dArrey[iSizeY+k]=m_dArrey[iSizeY+k]-m_dArrey[iSizeX+k]*dMultiplauer;
}
m_dArreyB[m_iIndex[z]]=m_dArreyB[m_iIndex[z]]-m_dArreyB[m_iIndex[x]]*dMultiplauer;
}
Но все же лучше в виде отдельной процедуры метод Гаусса реализовать ...