Арифметика высокой точности в 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