Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I am trying to use LAPACK geqp3 function for QR decomposition with pivoting but result is wrong. I tried almost two days but can't figure out the problem. I compared the result with Matlab and Python CvxOpt both are same and both uses the geqp3 LPACK routines as internally .

Here is input matrix and code. Any one can tell me what i am doing wrong ?

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

When i perform QR decomposition using geqp3 routines it gives me this result.

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

While in python and matlab i am getting this result.

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]
share|improve this question
2  
It's not clear to me from your post; what makes the result wrong? –  Geoff Oxberry Jan 13 '13 at 8:17
    
Also, exactly what output did your code produce? –  Brian Borchers Jan 14 '13 at 6:16
add comment

1 Answer

Interpreting the results

To those of you confused as to the difference, he is reporting the compact representation of the pivoted QR factorization, where the upper triangle gives the factor $R$ from the pivoted QR decomposition and each column of the strictly lower triangle represents a Householder reflection, say $$ H_j = I - \tau_j \left[\begin{array}{c}1\\ u_j\end{array}\right] \left[\begin{array}{cc}1 & u_j^H\end{array}\right],$$ where $u_j$ is the vector stored below the $j$'th diagonal entry of the result.

I recommend loading his input matrix into MATLAB or Octave in order to follow along at home. You will notice that the first result is essentially the same as an unpivoted QR factorization, e.g., via

X=qr(A)

but with signs of the upper triangle appropriately changed so that the diagonal entries are positive. This is a strong clue as to the problem.

Fixing the problem

If you look through the comments at the beginning of the dgeqp3 source, you will notice that if an entry of the JPVT array is nonzero, then that column is not allowed to be a pivot candidate. Notice that your code is simply allocating an integer array, which does not ensure that its contents are zero. If you explicitly zero the JPVT array before calling dgeqp3 then I suspect that your problem will go away.

share|improve this answer
1  
You are right it was the problem of JPVT. it working fine now thanx again for suggestion –  alex Jan 18 '13 at 8:15
    
Wonderful! Would you mind accepting my answer then? –  Jack Poulson Jan 18 '13 at 15:22
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.