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

Question

我希望使用来自C的LAPACK来完成MATLAB/Octave中的rcond操作.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;
}