lapack 子例程 dgeqrf 的 QR 因式分解在 C 中不起作用

Question

我想使用 lapack 从 A 矩阵的 QR 分解为 A=QR 中找到 R 矩阵的对角线元素。我尝试了 lapack dgeqrf 子例程，但它返回相同的矩阵 A，即输入和输出矩阵相同。如何找到 R 矩阵及其对角线？我不知道这段代码出了什么问题。我正在用 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;
}

Answer 1

矩阵 matA 未被修改，因为调用 LAPACK 的 dgeqrf() 时使用参数 LWORK 的值 -1。这对应于工作区查询：

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.

实际上，使用 dgeqrf() 和 LAPACK 中的许多其他例程的通常方法是调用它们两次：一次用于工作区查询，一次用于结果的实际计算。 例如，LAPACK 的 C 接口将 dgeqrf() 包装在 lapacke__dgeqrf() 中，正是出于这个原因，它调用了两次 lapacke__dgeqrf_work()。

修改代码的方法如下：

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

由gcc main.c -o main -llapack -lblas -lm.

编译

鉴于您提出的问题，标题为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，以便于检查结果的正确性。

lapack 子例程 dgeqrf 的 QR 因式分解在 C 中不起作用

QR factorization from lapack subroutine dgeqrf not working in C

lapack