Арифметика высокой точности в ALGLIB, как ее использовать? Опции

[MSP430F]

Всем доброго времени суток. Надеюсь, что мой вопрос не останется без ответа.
Использую функцию rmatrixsolvels из пакета alglib. Подключил просто - закинул в папку с проектом все хедеры и сишники alglib. Из 13 файлов пришлось подключить 6 для вызова лишь одной функции. exe подрос заметно в размерах, время компиляции увеличилось заметно, ну да ладно, главное цель достигнута. Решается система линейных уравнений с прямоугольной матрицей - число уравнений превышает число неизвестных. Минимизируется невязка. Но чувствую, не хватает точности. В alglib есть возможность работы с арифметикой высокой точности http://alglib.sources.ru/equations/linear.php Но как ее подключить ? Что-то я не соображу... Так все запутано в этом alglib. Использую компилятор MinGW-W64. Буду признателен за любые подсказки.

[thermit]

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:
Ошибка ~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

[MSP430F]

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:

Это Вы мне что, прямо таки код рабочий дали ? Спасибо большое, буду проверять. Если все 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

[thermit]

Это Вы мне что, прямо таки код рабочий дали ? Спасибо большое, буду проверять. Если все 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]; ?

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

_pv:
если то же самое в чём радость по сравнению с Гауссом?

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

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

[Guest_TSerg_*]

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

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