使用 lapack 计算倒数条件数(即 rcond(x))

Question

我希望使用 C 中的 LAPACK 来完成 rcond 在 MATLAB/Octave 中所做的事情。MATLAB 手册告诉我使用了 dgecon，即使用基于 1 的范数。

我为一个极其简单的案例写了一个简单的测试程序； [1,1; 1,0]对于这个输入 matlab 和 octave 使用 rcond 和 1/cond(x,1) 给我 0.25，但是在使用 LAPACK 的情况下，这个示例程序打印 0.0。对于其他情况，例如身份，它会打印正确的值。

由于假设 MATLAB 实际上成功地使用了这个例程，那么我做错了什么？我试图破译 Octave 的作用，但收效甚微，因为它包含在

#include <stdio.h>

extern void dgecon_(const char *norm, const int *n, const double *a, 
     const int *lda, const double *anorm, double *rcond, double *work, 
     int *iwork, int *info, int len_norm);

int main()
{
    int i, info, n, lda;
    double anorm, rcond;

    double w[8] = { 0,0,0,0,0,0,0,0 };

    int iw[2] = { 0,0 };

    double x[4] = { 1, 1, 1, 0 };
    anorm = 2.0; /* maximum column sum, computed manually */
    n = 2;
    lda = 2;

    dgecon_("1", &n, x, &lda, &anorm, &rcond, w, iw, &info, 1);

    if (info != 0) fprintf(stderr, "failure with error %d\n", info);
    printf("%.5e\n", rcond);
    return 0;
}

使用 cc testdgecon.c -o testdgecon -llapack 编译； ./testdgecon

Answer 1

我找到了自己问题的答案。

矩阵在发送到 dgecon 之前必须进行 LU 分解。这看起来非常合乎逻辑，因为人们经常希望在检查条件后求解系统，在这种情况下不需要将矩阵分解两次。同样的想法也适用于单独计算的范数。

以下代码是使用 LAPACK 计算倒数条件数的所有必要部分。

#include "stdio.h"

extern int dgecon_(const char *norm, const int *n, double *a, const int *lda, const double *anorm, double *rcond, double *work, int *iwork, int *info, int len);
extern int dgetrf_(const int *m, const int *n, double *a, const int *lda, int *lpiv, int *info);
extern double dlange_(const char *norm, const int *m, const int *n, const double *a, const int *lda, double *work, const int norm_len);

int main()
{
    int i, info, n, lda;
    double anorm, rcond;

    int iw[2];
    double w[8];
    double x[4] = {7,3,-9,2 };
    n = 2;
    lda = 2;

    /* Computes the norm of x */
    anorm = dlange_("1", &n, &n, x, &lda, w, 1);

    /* Modifies x in place with a LU decomposition */
    dgetrf_(&n, &n, x, &lda, iw, &info);
    if (info != 0) fprintf(stderr, "failure with error %d\n", info);

    /* Computes the reciprocal norm */
    dgecon_("1", &n, x, &lda, &anorm, &rcond, w, iw, &info, 1);
    if (info != 0) fprintf(stderr, "failure with error %d\n", info);

    printf("%.5e\n", rcond);
    return 0;
}