Это Вы мне что, прямо таки код рабочий дали ? Спасибо большое, буду проверять. Если все OK, то куда donation перевести ?

Так как же должно быть в BackwardSubstitution():

if(A[(C-1)*R+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*R+C-1]; или

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1]; ?


Сообщение отредактировал MSP430F - Dec 11 2013, 13:27

5*10^-10
если то же самое в чём радость по сравнению с Гауссом?

Последняя редакция верна. Первый раз скопировал с ошибкой.




На плохо обусловленных системах с числом уравнений большим числа неизвестных ошибка гаусса будет на порядок больше. Нужно будет принимать специальные меры.
Кроме того, применение ортогонального преобразования гарантирует состоятельность и численную устойчивость алгоритма. Про гаусса такого сказать нельзя.
Если число операций не является приоритетом, лучше использовать вращения или отражения.

зы
старик уилкинсон утверждает, что ему всегда хватало гаусса. Но это было в эпоху медленных компьютеров...

Не только.
Все чаще дело имеем с разреженными системами, а Гаусс разреженность лучше сохраняет, в отличии от методов отражений и вращений.

Я правильно догадался, что в этом месте ошибка:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];

Должно быть:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[i][j];

Первый индекс - строка, второй - столбец. Это же одинаково для всех массивов ?

Ошибки здесь нет. Исходная матрица лежит в памяти постолбцово (1-я строка состоит из единиц). Во всяком случае я так считал. Для работы алгоритмов матрица должна лежать в памяти построчно. for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i]; Это и есть такое преобразование (транспонирование).

Сообщение отредактировал thermit - Dec 12 2013, 09:19

Если считать, что первый индекс - строка, а второй - столбец, то первая строка как раз и не состоит из 1.0, а равна
{1.000000, 1.229402, 1.229566, 1.230016, 1.229648, 1.229771, 1.229975, 1.229361, 1.229525, 1.229648, 1.229484, 1.229238, 1.229156, 1.228992, 1.228214, 1.227518, 1.227928, 1.227108, 1.227231, 1.226453, 1.226576},
как она и проинициализирована. Зачем же здесь тогда транспонирование ?

Можете проверить:

i = 0;
for(j=0; j<21; j++) printf("a = %f\n", a[i][j]);

Сообщение отредактировал MSP430F - Dec 12 2013, 10:21

тогда простой Гаусс даст ошибку в 10^-12.
а для прямоугольной матрицы тем же Гауссом постойте псевдообратную и будет вам счастье и без ALGLIBа.

Для данных, которые вы прислали погрешность вычислялась (и с вашей совпадала) для матрицы, где строки являются столбцами из вашей таблицы. Поскольку я вставил ее себе в код копипастом, пришлось ввести транспонирование. С индексами все путем: 1-й - номер строки, 2-й - номер столбца.
Честно говоря, не понял проблемы. Разместите свои данные построчно да и все дела.



Собственно, то, что алглиб тут лишний, давно понятно.

Нет, Вы не правы. В моих данных все на месте - строки это строки, столбцы - столбцы. Первый элемент каждой строки - это 1.0 (так и должно быть). Транспонирование не нужно.

Это непринципиально. Кстати, что за данные? Откуда?

Проверяю метод ближайших соседей для прогнозирования значения временного ряда хаотического типа. Для примера взял котировки EUR/USD с шагом 1 час.

Ваш код хорошо работает. Правда, пришлось немного попотеть, чтобы можно было в качестве аргумента передавать матрицу в виде нормального двумерного массива, так как размерность матрицы может меняться динамически. Ваш код дает ошибку не хуже, чем alglib. При небольшой прямоугольности разница заметна в пользу Вашего кода. Если матрица 21 * 25 или 21 * 31, то разница в ошибке незаметна. Зато код получается небольшой и нет кучи warning-ов при компиляции, как при использовании alglib.

Зачем передавать в виде двумерного массива? Решателю достаточно передать адрес начала матрицы и ее размерность.

У топик-стартера матрица аллокируется не одним блоком памяти, а построчно:

Однако это обстоятельство ничуть не мешает ее передавать через указатель double**, а еще лучше, как &double**.
Хотя на C++ можно было бы создать для матрицы класс, где кроме указателя double**, хранились бы и обе ее размерности. А конструктор и деструктор память под эту матрицу аллокировали и деаллокировали. В этом случае ссылку (здесь она удобнее указателя) на матрицу можно передавать одним единственным парамером. Кроме того, данный класс можно снабдить фукциями членами класса дял матричного умножения, сложения и т.п.