计算科学 - 在 LAPACK 中使用 Pivoting 的错误结果 QR 分解 - 吾爱随笔录

我正在尝试使用 LAPACK geqp3 函数进行带有旋转的 QR 分解，但结果是错误的。我尝试了将近两天，但无法弄清楚问题所在。我将结果与 Matlab 和 Python CvxOpt 进行了比较，两者都相同，并且都在内部使用 geqp3 LPACK 例程。

这是输入矩阵和代码。任何人都可以告诉我我做错了什么？

> Blockquote A =
   -0.9586   -0.8171   -0.6033    1.1895    1.3031   -0.0759
    1.9877    1.0979    0.7795   -0.7164   -0.0505   -0.2412
    0.2328    0.2210    0.9401    0.0691    0.1082   -0.8700
    0.2808    0.9539    0.0198   -0.2278   -1.1498   -0.3127
    0.1335    1.2238    0.5449    0.2712   -0.7684   -2.5430
    0.3996   -0.5725    0.7090   -1.2492    0.5508   -1.0513
   -0.3073    0.2512   -1.2444    0.1169   -1.1271    0.6893
    0.3238   -0.8187    2.7957    0.3582    1.8493   -1.5166
    0.6406   -0.4428   -0.1864   -1.0097    0.6713    0.7002


void fReadMtrx(string f, double *mtrx, int nrow, int ncol){

    ifstream data(f.c_str(), ios::in);
    if (!data) {
        cerr << "File " << f << " could not be opened." << endl;
        exit(1);
    }

    for(int i = 0; i < nrow; i++) 
        for(int j = 0; j < ncol; j++) 
            data >> mtrx[j*nrow+i]; 
}

void showMtrx(double *mtrx, int nrow, int ncol){

    for (int i = 0; i < nrow; i++) {
        for (int j = 0; j < ncol; j++) { 
            cout <<  mtrx[j*nrow+i] << "  ";
        }
        cout << endl;
    }
}


int main(int argc, char *argv[])
{
    int m = 9;
    int n = 6;

    //allocate and read in a
    double* a = new double[m*n];

    fReadMtrx("A.dat", a, m, n);

    cout << endl << "a:" << endl;
    showMtrx(a, m, n);

    int ldA = m;
    double wl;
    int lwork;
    int info;


    // QR decomposition
    lwork = -1;
    dgeqp3_(&m, &n, NULL, &ldA, NULL, NULL, &wl, &lwork, &info);

    lwork = (int)wl;
    double* work  = new double[lwork];
    int* jpvt_ptr = new int[n];
    double* tau = new double[n];

    dgeqp3_(&m, &n, a, &ldA, jpvt_ptr,tau,(double *)work, &lwork, &info);

    cout << endl << "QR:" << endl;
    showMtrx(a, m, n);


   return 0;
}

当我使用 geqp3 例程执行 QR 分解时，它给了我这个结果。

QR :
 2.4063   1.0778   1.6110  -1.5135  -0.0642  -0.7111 
-0.5907   2.1035  -1.1498   0.4283  -2.6023  -0.6335 
-0.0692  -0.0423   2.8526   0.7086   0.9759  -2.7723
-0.0835  -0.3745  -0.0332   1.3431   0.8472   0.5267
-0.0397  -0.5405  -0.3887  -0.3816   0.7893  -1.2573
-0.1188   0.3753  -0.1044   0.9602  -0.0370   1.1425
 0.0913  -0.1996   0.4502  -0.1306   2.1257  -0.1959
-0.0962   0.4710  -1.1228   0.0682   1.4831  -0.4424
-0.1904   0.3782   0.4133   0.4090  -1.0372   0.2011

在 python 和 matlab 中，我得到了这个结果。

QR :
   [ 3.47196268, 1.63382425, -2.39788322, -0.26189627,  1.11652819,  -0.1965108 ],
   [-0.19128521,  2.52976731,  0.93296076,  0.49012156, -0.78217128, -2.06427021],
   [-0.23067443, -0.0120604 ,  2.26624645, -0.67846333,  0.74839188,  -0.28431637],
   [-0.00486702,  0.43558037, -0.00634783, -2.00278565,  1.22747977,  -0.01853067],
   [-0.13370508,  0.3074231 ,  0.76958444, -0.0335343 ,  1.36566733,   1.04995806],
   [-0.17396713, -0.18659867,  0.30029599, -0.6800568 ,  0.37513732,   0.3212295 ],
   [ 0.30534072,  0.38815978, -0.10102664,  0.06497129,  0.35993755,   0.72318025],
   [-0.68601098, -0.61373672,  0.14243916,  0.28777571, -0.68376182,  0.3557272 ],
   [ 0.04573902, -0.25965053, -0.15498792, -0.52485689, -0.0986465 ,  -0.38370839]