| |
 Арифметика высокой точности в ALGLIB, как ее использовать? |
|
|
|
|
 Dec 10 2013, 09:06
|
Частый гость
Группа: Участник
Сообщений: 85
Регистрация: 20-05-13
Пользователь №: 76 911
|
Всем доброго времени суток. Надеюсь, что мой вопрос не останется без ответа. Использую функцию rmatrixsolvels из пакета alglib. Подключил просто - закинул в папку с проектом все хедеры и сишники alglib. Из 13 файлов пришлось подключить 6 для вызова лишь одной функции. exe подрос заметно в размерах, время компиляции увеличилось заметно, ну да ладно, главное цель достигнута. Решается система линейных уравнений с прямоугольной матрицей - число уравнений превышает число неизвестных. Минимизируется невязка. Но чувствую, не хватает точности. В alglib есть возможность работы с арифметикой высокой точности http://alglib.sources.ru/equations/linear.php Но как ее подключить ? Что-то я не соображу... Так все запутано в этом alglib. Использую компилятор MinGW-W64. Буду признателен за любые подсказки. 
|
|
|
|
|
|
|
|
Dec 10 2013, 10:26
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Цитата MSP430F:
Но чувствую, не хватает точности. Каким образом почувствовали? Какого размера система? Каким способом решается? Вероятно, необходимость в высокой точности необоснована.
|
|
|
|
|
|
|
Guest_TSerg_*
|
Dec 10 2013, 11:49
|
Guests
|
Проверьте с использованием мат.пакетов Матлаб, Маткад..
|
|
|
|
|
|
|
|
Dec 10 2013, 12:32
|
Частый гость
Группа: Участник
Сообщений: 85
Регистрация: 20-05-13
Пользователь №: 76 911
|
Цитата(thermit @ Dec 10 2013, 14:26)  Каким образом почувствовали? Какого размера система? Каким способом решается?
Вероятно, необходимость в высокой точности необоснована. Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки, то что же ожидать от прямоугольной ?
|
|
|
|
|
|
|
|
Dec 10 2013, 13:47
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Цитата MSP430F:
Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ? Исходные данные можете предоставить?
|
|
|
|
|
|
|
|
Dec 10 2013, 14:33
|
Частый гость
Группа: Участник
Сообщений: 85
Регистрация: 20-05-13
Пользователь №: 76 911
|
Цитата(Xenia @ Dec 10 2013, 16:47) Если у вас на квадратной матрице rmatrixsolvels дает ошибку от 10-5 до 10-3, то при переходе на прямоугольною матрицу хуже уже не будет. Ухудшение на "3-5 порядков" - это плата за метод более общего типа (в нем арифметики больше), которую вы уже заплатили. Спасибо за ответ! Я чувствовал нечто подобное, а Вы подтвердили мои догадки. Да, похоже, что использовать арифметику высокой точности по-любому очень проблемно здесь. Но если бы кто, подсказал, как, то можно было и попробовать. Если бы можно было просто заменить double на long double... Но там иакой финт не пройдет, слишком все в alglib наворочено. Цитата(thermit @ Dec 10 2013, 17:47) Исходные данные можете предоставить? Да. В файле матрица 21x21 элемент и массив B из 21-го элемента (правая часть). Максимальная ошибка для этих данных соответственно у меня получилась 1,1х10-9 и 1,9х10-4.
|
|
|
|
|
|
|
|
Dec 10 2013, 16:27
|
Гуру
Группа: Модератор FTP
Сообщений: 4 479
Регистрация: 20-02-08
Из: Москва
Пользователь №: 35 237
|
Цитата(thermit @ Dec 10 2013, 19:31) Боюсь, алглиб вас не спасет. Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. На вашем примере гивенс дает максимальную погрешность 1.5*10^-9. Alglib действительно не спасет, там алгоритм для этого дела гадкий . За Хаусхолдера же бояться нечего, это ортогональное преобразование, и ошибок оно практически не вносит (лишь на уровне точности арифметики). А вот на стадии обратной подстановки при плохо обусловленной матрице неминуемо возникнут проблемы, т.к. обусловленность полученной треугольной матрицы останется такой же малой, как у исходной. Боюсь, что в данном случае надо считать через сингулярное разложение исходной матрицы, а то и понижая при этом размерность. Т.е. получать псевдорешение. Впрочем, я бы рекомендовала сначала не заморачиваться всем этим, а применить процедуру rmatrixsolvels, как она есть, к прямоугольной матрице, на которой, по моим ощущениям, ошибка должна ученьшиться по сравнению с квадратной.
|
|
|
|
|
|
|
|
Dec 10 2013, 16:44
|
Частый гость
Группа: Участник
Сообщений: 85
Регистрация: 20-05-13
Пользователь №: 76 911
|
Цитата(thermit @ Dec 10 2013, 19:31) Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной... Цитата(Xenia @ Dec 10 2013, 20:27) Боюсь, что в данном случае надо считать через сингулярное разложение исходной матрицы. Но ведь Alglib так и считает...
Сообщение отредактировал MSP430F - Dec 10 2013, 16:45
|
|
|
|
|
|
|
|
Dec 10 2013, 19:03
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Цитата MSP430F:
И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной... Можно легко модифицировать для прямоугольной. Если исходники есть. ps Могу нарисовать за умеренную плату, гы...
|
|
|
|
|
|
|
|
Dec 11 2013, 07:51
|
Частый гость
Группа: Участник
Сообщений: 85
Регистрация: 20-05-13
Пользователь №: 76 911
|
Цитата(thermit @ Dec 10 2013, 23:03) Можно легко модифицировать для прямоугольной. Если исходники есть. Вот гляньте этот файл. Эта функция дает ошибку 2.7*10-8 на тех же исходных данных. Она не оригинальна, чуть модифицирована отсюда Мной добавлено создание копии матрицы в теле функции, а создание массива результатов убрано. Также введена проверка деления на 0.
Сообщение отредактировал MSP430F - Dec 11 2013, 07:54
|
|
|
|
|
|
|
|
Dec 11 2013, 11:00
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Честно говоря, лень ковыряться с чужими исходниками. Должно быть как-то так: Ошибка ~9*10^-10 main.c Код #include <stdio.h>
#include <math.h>
int SolveByGivensRotation(double *A,double *b,int R,int C,double *x);
double a[21][21]=
{
{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},
{1.000000, 1.229220, 1.228551, 1.229852, 1.229774, 1.230227, 1.230539, 1.230034, 1.229774, 1.229279, 1.228967, 1.229019, 1.228978, 1.228677, 1.228413, 1.227573, 1.227919, 1.227521, 1.226982, 1.226752, 1.226759},
{1.000000, 1.229566, 1.229689, 1.229846, 1.230118, 1.229974, 1.229979, 1.229705, 1.229582, 1.229385, 1.229313, 1.229417, 1.228537, 1.228060, 1.228113, 1.227828, 1.227652, 1.227244, 1.227276, 1.226674, 1.226855},
{1.000000, 1.228701, 1.229879, 1.229808, 1.230219, 1.230501, 1.230044, 1.229808, 1.229361, 1.229078, 1.229125, 1.229088, 1.228816, 1.228577, 1.227816, 1.228129, 1.227769, 1.227281, 1.227072, 1.227079, 1.226659},
{1.000000, 1.229696, 1.229674, 1.229618, 1.229730, 1.229786, 1.229718, 1.229696, 1.229607, 1.229741, 1.229528, 1.229237, 1.229349, 1.228286, 1.227962, 1.227481, 1.227738, 1.227279, 1.226708, 1.226932, 1.227044},
{1.000000, 1.229095, 1.229253, 1.230450, 1.229750, 1.229885, 1.229389, 1.229389, 1.229479, 1.229614, 1.229547, 1.229140, 1.229253, 1.228824, 1.228576, 1.227943, 1.227966, 1.227153, 1.226860, 1.226815, 1.226431},
{1.000000, 1.228654, 1.228773, 1.230379, 1.230498, 1.229963, 1.229784, 1.229190, 1.230022, 1.229784, 1.229071, 1.228773, 1.228952, 1.228297, 1.228714, 1.228119, 1.227822, 1.226810, 1.227286, 1.226929, 1.226989},
{1.000000, 1.229305, 1.229601, 1.229552, 1.229848, 1.230292, 1.229453, 1.229552, 1.229700, 1.229798, 1.229650, 1.228713, 1.228368, 1.229108, 1.228516, 1.227480, 1.228072, 1.227185, 1.227037, 1.227234, 1.226346},
{1.000000, 1.228883, 1.228629, 1.230443, 1.230170, 1.230034, 1.230053, 1.230190, 1.229663, 1.228766, 1.228766, 1.228961, 1.229000, 1.228746, 1.228473, 1.227440, 1.227849, 1.228122, 1.227089, 1.226640, 1.226893},
{1.000000, 1.229543, 1.229630, 1.229356, 1.229499, 1.229683, 1.230001, 1.229833, 1.229839, 1.229518, 1.229374, 1.229143, 1.229059, 1.229181, 1.228151, 1.227593, 1.227655, 1.227322, 1.227116, 1.226638, 1.226676},
{1.000000, 1.229429, 1.229542, 1.230312, 1.229936, 1.230049, 1.230105, 1.229336, 1.229148, 1.229054, 1.228716, 1.228491, 1.229636, 1.228359, 1.228453, 1.227590, 1.228266, 1.228040, 1.227609, 1.226276, 1.226464},
{1.000000, 1.229692, 1.229447, 1.229575, 1.229739, 1.230022, 1.229872, 1.229878, 1.229592, 1.229464, 1.229258, 1.229183, 1.229292, 1.228375, 1.227878, 1.227933, 1.227636, 1.227453, 1.227027, 1.227061, 1.226433},
{1.000000, 1.229186, 1.229563, 1.229939, 1.230316, 1.229971, 1.229469, 1.229312, 1.229374, 1.229531, 1.229469, 1.229280, 1.229061, 1.228527, 1.228810, 1.227617, 1.227335, 1.226644, 1.226833, 1.227021, 1.227554},
{1.000000, 1.229820, 1.229803, 1.230004, 1.229990, 1.229898, 1.229716, 1.229279, 1.229198, 1.229403, 1.229509, 1.229152, 1.229086, 1.228881, 1.227502, 1.227404, 1.227758, 1.227441, 1.227067, 1.226886, 1.227012},
{1.000000, 1.229545, 1.229452, 1.229675, 1.229582, 1.229861, 1.229805, 1.229712, 1.229694, 1.229824, 1.229285, 1.229174, 1.229025, 1.228413, 1.228394, 1.228227, 1.227707, 1.227206, 1.227002, 1.226575, 1.226649},
{1.000000, 1.229498, 1.229330, 1.229305, 1.229378, 1.229955, 1.230075, 1.229882, 1.229642, 1.229522, 1.229618, 1.229354, 1.228993, 1.228176, 1.228945, 1.227743, 1.227310, 1.227671, 1.226949, 1.226829, 1.226637},
{1.000000, 1.229748, 1.229748, 1.229999, 1.229971, 1.229776, 1.229608, 1.229385, 1.229469, 1.229497, 1.229329, 1.229134, 1.228799, 1.228687, 1.228576, 1.227516, 1.227237, 1.227432, 1.227460, 1.227153, 1.226288},
{1.000000, 1.228414, 1.229260, 1.229719, 1.229719, 1.230036, 1.230000, 1.229754, 1.229542, 1.229260, 1.229366, 1.229401, 1.229190, 1.228943, 1.228520, 1.228379, 1.228238, 1.226899, 1.226547, 1.226793, 1.226828},
{1.000000, 1.229178, 1.228977, 1.229153, 1.230486, 1.229706, 1.229857, 1.229304, 1.229304, 1.229405, 1.229556, 1.229480, 1.229027, 1.229153, 1.228675, 1.228398, 1.227694, 1.227719, 1.226814, 1.226487, 1.226436},
{1.000000, 1.229249, 1.229674, 1.229702, 1.229759, 1.229957, 1.229759, 1.229702, 1.229787, 1.229702, 1.229475, 1.229192, 1.228399, 1.227946, 1.227776, 1.227890, 1.228314, 1.227691, 1.228059, 1.226360, 1.226417},
{1.000000, 1.228960, 1.229996, 1.230061, 1.230319, 1.230190, 1.230190, 1.229219, 1.229414, 1.228960, 1.228766, 1.228572, 1.229090, 1.229155, 1.227990, 1.227860, 1.227666, 1.227278, 1.226890, 1.227278, 1.226955}
};
double b[21]=
{
1.226985, 1.226294, 1.226775, 1.225659, 1.227156, 1.226476, 1.225978, 1.225508, 1.226484, 1.225971, 1.226426, 1.226622, 1.224824, 1.226897, 1.226687, 1.226372, 1.226260, 1.226441, 1.226009, 1.226643, 1.226760
};
double aorig[21][21],aa[21][21],x[21],bb[21];
int main()
{
int i,j;
double acc,maxerror;
for(i=0;i<21;i++)
{
for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];
bb[i]=b[i];
}
SolveByGivensRotation(&aa[0][0],bb,21,21,x);
maxerror=0.0;
for(i=0;i<21;i++)
{
acc=b[i];
for(j=0;j<21;j++) acc-=aorig[i][j]*x[j];
if(fabs(acc)>maxerror) maxerror=fabs(acc);
}
printf("Max error = %.15g\n",maxerror);
return 0;
} givens.c Код #include <math.h>
#define SQR(x) (x*x)
int Rotation(double *A,double *b,int R,int C,int ri,int ci)
{
int i;
double t1,t2,t3,c,s;
t1=A[C*ci+ci];
t2=A[C*ri+ci];
t3=sqrt(SQR(t1)+SQR(t2));
if(t3==0.0) return 1;
c=t1/t3;
s=t2/t3;
for(i=0;i<C;i++)
{
t1=c*A[ci*C+i]+s*A[ri*C+i];
t2=c*A[ri*C+i]-s*A[ci*C+i];
A[ci*C+i]=t1;
A[ri*C+i]=t2;
}
t1=c*b[ci]+s*b[ri];
t2=c*b[ri]-s*b[ci];
b[ci]=t1;
b[ri]=t2;
return 0;
}
int GivensTriang(double *A,double *b,int R,int C)
{
int i,j;
for(i=0;i<C;i++)//For col
for(j=R-1;j>i;j--) //For row
if(Rotation(A,b,R,C,j,i)) return 1;
return 0;
}
void BackwardSubstitution(double *A,double *b,double *x,int R,int C)
{
int i,j;
double tmp;
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];
for(i=C-2;i>=0;i--)
{
tmp=b[i];
for(j=i+1;j<C;j++) tmp-=x[j]*A[i*C+j];
if(A[i*C+i]==0.0) x[i]=0.0;
else x[i]=tmp/A[i*C+i];
}
}
int SolveByGivensRotation(double *A,double *b,int R,int C,double *x)
{
if(GivensTriang(A,b,R,C)) return 1;
BackwardSubstitution(A,b,x,R,C);
return 0;
}
Сообщение отредактировал thermit - Dec 11 2013, 12:48
|
|
|
|
|
|
|
|
Dec 11 2013, 12:07
|
Знающий
Группа: Участник
Сообщений: 781
Регистрация: 3-08-09
Пользователь №: 51 730
|
Цитата _pv:
обычный Гаусс на этих данных даёт ошибку 10^-12. Не дает. Те же самые ~9*10^-10. Матрицу нужно транспонировать.
|
|
|
|
|
|
1 чел. читают эту тему (гостей: 1, скрытых пользователей: 0)
Пользователей: 0
|
|
|
