当前位置: 动力学知识库 > 问答 > 编程问答 >

c - Is this a lapack problem or maybe a bug in my code?

问题描述:

I am writing c code to get get the QR decomposition for a matrix A using lapack lib in R.

My result is different from the one i get using R language command.

Is this a lapack thing or maybe a bug in my code ?

for matrix (row major) :

1, 19.5, 43.1, 29.1,

1, 24.7, 49.8, 28.2,

1, 30.7, 51.9, 37.0,

1, 29.8, 54.3, 31.1,

1, 19.1, 42.2, 30.9,

1, 25.6, 53.9, 23.7,

1, 31.4, 58.5, 27.6,

1, 27.9, 52.1, 30.6,

1, 22.1, 49.9, 23.2,

1, 25.5, 53.5, 24.8,

1, 31.1, 56.6, 30.0,

1, 30.4, 56.7, 28.3,

1, 18.7, 46.5, 23.0,

1, 19.7, 44.2, 28.6,

1, 14.6, 42.7, 21.3,

1, 29.5, 54.4, 30.1,

1, 27.7, 55.3, 25.7,

1, 30.2, 58.6, 24.6,

1, 22.7, 48.2, 27.1,

1, 25.2, 51.0, 27.5

r result is :

 V1 V2 V3 V4

[1,] -4.4721360 -113.16740034 -2.288392e+02 -123.52039508

[2,] 0.2236068 -21.89587861 -2.107945e+01 -7.27753395

[3,] 0.2236068 0.29484219 8.733781e+00 -14.04825478

[4,] 0.2236068 0.25373857 7.566965e-02 -1.55436071

[5,] 0.2236068 -0.23493787 2.999600e-01 0.26555995

[6,] 0.2236068 0.06192165 -3.343037e-01 0.12188660

[7,] 0.2236068 0.32681168 -2.315941e-01 0.40765540

[8,] 0.2236068 0.16696425 1.213823e-01 -0.18580207

[9,] 0.2236068 -0.09792578 -2.561224e-01 -0.16369010

[10,] 0.2236068 0.05735458 -2.993562e-01 0.52588892

[11,] 0.2236068 0.31311047 -4.660317e-02 0.29409317

[12,] 0.2236068 0.28114099 -1.340150e-01 0.16746961

[13,] 0.2236068 -0.25320615 -2.357881e-01 0.26072358

[14,] 0.2236068 -0.20753545 1.360745e-01 0.18135493

[15,] 0.2236068 -0.44045600 -2.456167e-01 0.15393180

[16,] 0.2236068 0.24003736 3.166468e-02 -0.02119950

[17,] 0.2236068 0.15783011 -2.667146e-01 0.32042553

[18,] 0.2236068 0.27200685 -3.732646e-01 0.05926317

[19,] 0.2236068 -0.07052337 3.634501e-03 -0.20518296

[20,] 0.2236068 0.04365337 -4.566657e-02 -0.03457051

my result:

-4.4721, -113.1674, -228.8392, -123.5204,

0.1827, -21.8959, -21.0794, -7.2775,

0.1827, 0.2888, 8.7338, -14.0483,

0.1827, 0.2486, 0.0523, -1.5544,

0.1827, -0.2301, 0.2071, 0.2316,

0.1827, 0.0607, -0.2309, 0.1063,

0.1827, 0.3201, -0.1599, 0.3555,

0.1827, 0.1636, 0.0838, -0.1620,

0.1827, -0.0959, -0.1769, -0.1427,

0.1827, 0.0562, -0.2067, 0.4586,

0.1827, 0.3067, -0.0322, 0.2565,

0.1827, 0.2754, -0.0925, 0.1460,

0.1827, -0.2480, -0.1628, 0.2274,

0.1827, -0.2033, 0.0940, 0.1581,

0.1827, -0.4315, -0.1696, 0.1342,

0.1827, 0.2351, 0.0219, -0.0185,

0.1827, 0.1546, -0.1842, 0.2794,

0.1827, 0.2665, -0.2578, 0.0517,

0.1827, -0.0691, 0.0025, -0.1789,

0.1827, 0.0428, -0.0315, -0.0301,

#include <stdio.h>

#include <R.h>

#include <R_ext/BLAS.h>

#include <R_ext/Lapack.h>

int min(int x, int y) {

if (x <= y)

return x;

else

return y;

}

int max(int x, int y) {

if (x >= y)

return x;

else

return y;

}

void transpose(int* nrow, int* ncol, double* a) {

int i, j, index, k = 0;

double* atransp = malloc(*nrow**ncol * sizeof (double));

//compute transpose

for (i = 0; i<*ncol; i++) {

index = i;

atransp[k] = a[index];

k++;

for (j = 0; j<*nrow - 1; j++) {

index += *ncol;

atransp[k] = a[index];

k++;

}

}

//copy transpose in array a

for (i = 0; i<*nrow**ncol; i++)

a[i] = atransp[i];

//free memory

free(atransp);

}

void getQR(int* rowX, int* colX, double* X, double* Tau) {

const int m = *rowX;

const int n = *colX;

double* a = X;

const int lda = max(1, m);

double* tau = malloc(min(m, n) * sizeof (double));

const int lwork = max(1, n);

double* work = malloc(max(1, lwork) * sizeof (double));

int info;

F77_NAME(dgeqrf)(&m, &n, a, &lda, tau, work, &lwork, &info);

printf("\n dgeqrf() ended with : %d\n",info);

copyTo(min(m, n), tau, Tau);

free(work);

free(tau);

}

int main() {

int rX = 20, cX = 4;

double X[] = {

1, 19.5, 43.1, 29.1,

1, 24.7, 49.8, 28.2,

1, 30.7, 51.9, 37.0,

1, 29.8, 54.3, 31.1,

1, 19.1, 42.2, 30.9,

1, 25.6, 53.9, 23.7,

1, 31.4, 58.5, 27.6,

1, 27.9, 52.1, 30.6,

1, 22.1, 49.9, 23.2,

1, 25.5, 53.5, 24.8,

1, 31.1, 56.6, 30.0,

1, 30.4, 56.7, 28.3,

1, 18.7, 46.5, 23.0,

1, 19.7, 44.2, 28.6,

1, 14.6, 42.7, 21.3,

1, 29.5, 54.4, 30.1,

1, 27.7, 55.3, 25.7,

1, 30.2, 58.6, 24.6,

1, 22.7, 48.2, 27.1,

1, 25.2, 51.0, 27.5

};

//column major

transpose(&rX,&cX,X);

//tau is needed to extract Q later

double* tau = malloc(min(rX, cX) * sizeof (double));

double* QR = malloc(rX*cX*sizeof(double));

copyTo(rX*cX,X,QR);

getQR(&rX, &cX, QR, tau);

//printmat(cX,rX,QR,"QR");

return 0;

}

网友答案:

First of all, R uses the LINPACK routine DQRDC2 as a default. You can use the qr() command in R with the LAPACK=TRUE option to use a LAPACK routine :

> QR <- qr(Mat,LAPACK=TRUE)
> QR
$qr
                V3            V4            V2          V1
 [1,] -229.9739116 -123.04447843 -114.61600066 -4.45006998
 [2,]    0.1823682  -19.23736803   -1.81750327 -0.25177355
 [3,]    0.1900584    0.41052447  -12.08962672  0.36451274
 [4,]    0.1988473    0.04298840    0.18492760  0.02485389
 [5,]    0.1545369    0.37519893   -0.14576039  0.48992952
 [6,]    0.1973825   -0.32149888   -0.01277324 -0.02631498
 [7,]    0.2142277   -0.25359559    0.19391630  0.15368409
 [8,]    0.1907908    0.07984478    0.13051129 -0.21692785
 [9,]    0.1827344   -0.23371181   -0.11707414 -0.19513972
[10,]    0.1959177   -0.25431801   -0.01531797  0.38150533
[11,]    0.2072699   -0.07795264    0.21041668  0.11596815
[12,]    0.2076361   -0.16711575    0.17604756 -0.02208373
[13,]    0.1702836   -0.14766629   -0.23298857  0.30390902
[14,]    0.1618609    0.20180495   -0.14729488  0.33577876
[15,]    0.1563679   -0.12647957   -0.37138938  0.29164143
[16,]    0.1992135   -0.01062557    0.17041958 -0.12582222
[17,]    0.2025093   -0.25954266    0.06531149  0.14217017
[18,]    0.2145939   -0.40877853    0.13742162 -0.20783552
[19,]    0.1765090    0.01244890   -0.06067115 -0.15702073
[20,]    0.1867626   -0.04646282    0.01506042 -0.06657167

However, you should be aware that the function qr() in this case uses the LAPACK routine DGEQP3. Contrary to the DGEQRF routine you're using, DGEQP3 calculates the QR decomposition of a matrix with column pivoting.

So pretty normal you get different results, as you're not using the same methods. You should remember that a QR decomposition is not a unique solution. To know whether your QR decomposition is correct, you can simply check if the Q and R matrices fulfill the requirements. eg in R :

> Q <- qr.Q(QR)
> round( t(Q) %*% Q , 10 )
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

> all.equal(Q %*% qr.R(QR),Mat)
[1] TRUE

See also ?qr.

分享给朋友:
您可能感兴趣的文章:
随机阅读: