来自 lapack 子程序 dgeqrf 的 QR 分解在 C 中不起作用

Question

我想使用 lapack 找到从 A 矩阵的 QR 分解中获得的 R 矩阵的对角元素，如 A=QR。我尝试了 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 调用的。这对应于工作区查询：

如果 LWORK = -1，则假定为工作区查询；例行公事 只计算 WORK 数组的最佳大小，返回 此值作为 WORK 数组的第一个条目，并且没有错误 与 LWORK 相关的消息由 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 的问题，行列式的计算被添加到代码的末尾。矩阵matA现在是Vandermonde matrix，方便检查结果的正确性。