QR factorization from lapack subroutine dgeqrf not working in C

Question

I want to find diagonal elements of R matrix obtained from QR-decomposition of A matrix as A=QR using lapack. I tried lapack dgeqrf subroutine but it is giving back the same matrix A i.e. input and output matrices are same. How to find R matrix and its diagonals ? I can't figure out what is going wrong with this code. I am programming in C.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void dgeqrf_(int *rows, int *cols, double *matA, int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO);
int main()
{
    int rows=3;
    int cols=3;
    double *matA=malloc(sizeof(double)*rows*cols);

    matA[0]=10;
    matA[1]=20;
    matA[2]=10;
    matA[3]=40;
    matA[4]=20;
    matA[5]=50;
    matA[6]=70;
    matA[7]=30;
    matA[8]=20;

    for(int i=0; i<rows*cols; i++)
    {
        printf("%f ",matA[i]);
    }
    printf("\n");
    int LDA=rows;
    int INFO;
    double *WORK=malloc(sizeof(double)*2);
    int LWORK=-1;
    int rowcolmin=rows;
    if(rowcolmin>cols)
    {
        rowcolmin=cols;
    }
    double *TAU=malloc(sizeof(double)*rowcolmin);
    dgeqrf_(&rows, &cols, matA, &LDA, TAU, WORK, &LWORK, &INFO);
    for(int i=0; i<rows*cols; i++)
    {
        printf("%f ",matA[i]);
    }
    printf("\n");
    free(WORK);
    free(TAU);
    free(matA);
    return 0;
}

francis francis · Accepted Answer · 2017-11-17T11:14:43

The matrix matA is not modified because LAPACK's dgeqrf() is called using a value of -1 for the argument LWORK. This correspond to a workspace query:

If LWORK = -1, then a workspace query is assumed; the routine only calculates the optimal size of the WORK array, returns this value as the first entry of the WORK array, and no error message related to LWORK is issued by XERBLA.

Indeed, the usual way to use dgeqrf() and many other routines from LAPACK is to call them twice: once for the workspace query and once for the actual computation of the result. For instance, the C interface to LAPACK wraps dgeqrf() in lapacke__dgeqrf(), which calls lapacke__dgeqrf_work() twice for this very reason.

Here is how your code could be modified:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
void dgeqrf_(int *rows, int *cols, double *matA, int *LDA, double *TAU, double *WORK, int *LWORK, int *INFO);
int main()
{
    int i;
    int rows=3;
    int cols=3;
    double *matA=malloc(sizeof(double)*rows*cols);

    matA[0]=1;
    matA[1]=2;
    matA[2]=4;
    matA[3]=1;
    matA[4]=3;
    matA[5]=9;
    matA[6]=1;
    matA[7]=4;
    matA[8]=16;

    for(i=0; i<rows*cols; i++)
    {
        printf("%f ",matA[i]);
    }
    printf("\n");
    int LDA=rows;
    int INFO;

    int rowcolmin=rows;
    if(rowcolmin>cols)
    {
        rowcolmin=cols;
    }

    double *TAU=malloc(sizeof(double)*rowcolmin);
    int LWORK=-1;
    // since the value of LWORK is -1, this is a workspace query.
    // it only return the optimum size of work in work[0].
    double lwork2;    
    dgeqrf_(&rows, &cols, matA, &LDA, TAU, &lwork2, &LWORK, &INFO);
    // allocate work using the optimal size
    int lwork3=(int)lwork2;
    double *WORK=malloc(sizeof(double)*lwork3);
    // perform the QR factorization
    dgeqrf_(&rows, &cols, matA, &LDA, TAU, WORK, &lwork3, &INFO);
    if(INFO !=0){fprintf(stderr,"QR factorization failed, error code %d\n",INFO);exit(1);}

    for(i=0; i<rows*cols; i++)
    {
        printf("%f ",matA[i]);
    }
    printf("\n");


    // for instance, the determinant is...
    if(cols==rows){
        // det of R
        double det=1;
        for (i=0;i<cols;i++){
            det*=matA[i*cols+i];
        }
        // det of Q, Householder algorithm
        if(cols%2==0){
            det*=-1;
        }
        printf("det is %g\n",det);
    }
    free(WORK);
    free(TAU);
    free(matA);
    return 0;
}

It is compiled by gcc main.c -o main -llapack -lblas -lm.

Given the question you asked, entitled compute determinant from LU decomposition in lapack , the computation of the determinant is added at the end of the code. The matrix matA is now a Vandermonde matrix, so as to easily check the correctness of the result.

QR factorization from lapack subroutine dgeqrf not working in C

1 Answers