这是一个lapack问题，或者可能是我的代码中的错误？

Question

我正在写c代码来获得使用R中的lapack库获得矩阵A的QR分解。 我的结果与使用R语言命令得到的结果不同。

这是一个lapack的东西，或者可能是我的代码中的错误？

矩阵（行主要）：

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结果是：

   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 

我的结果：

-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; 
} 

Answer 1

首先，R使用LINPACK常规DQRDC2为默认。在R中可以使用qr()命令与LAPACK=TRUE选项使用LAPACK例行：

> 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 

但是，你应该知道，在这种情况下，功能qr()使用LAPACK例行DGEQP3。与您正在使用的DGEQRF例程相反，DGEQP3通过列旋转来计算矩阵的QR分解。

非常正常，你会得到不同的结果，因为你没有使用相同的方法。你应该记住QR分解是而不是一个独特的解决方案。要知道您的QR分解是否正确，您可以简单地检查Q和R矩阵是否满足要求。例如在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 

另请参见?qr。