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