7

C から LAPACK を使用して、MATLAB/Octave で rcond が行うことを正確に実行したいと考えています。

非常に単純なケース用の単純なテスト プログラムを作成しました。[1,1; 1,0] この入力の matlab と octave では、rcond と 1/cond(x,1) を使用すると 0.25 が得られますが、LAPACK を使用する場合、このサンプル プログラムは 0.0 を出力します。ID などのその他のケースでは、正しい値が出力されます。

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

4

1 に答える 1

11

私自身の質問に対する答えを見つけました。

行列は、dgecon に送信する前に LU 分解する必要があります。条件を確認した後にシステムを解きたいことがよくあるため、これは非常に論理的に思えます。この場合、行列を 2 回分解する必要はありません。同じ考え方が、別々に計算されるノルムにも当てはまります。

次のコードは、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;
}
于 2011-02-12T01:45:49.660 に答える