Вычисление числа взаимного условия с помощью lapack (то есть, rcond (x)) - PullRequest
Новая
       32

Вычисление числа взаимного условия с помощью lapack (то есть, rcond (x))

7 голосов
/ 12 февраля 2011

Я хочу сделать именно то, что делает rcond в MATLAB / Octave, используя LAPACK из C. В руководстве MATLAB сказано, что dgecon используется, и это основано на норме, основанной на 1.

Я написал простой тестпрограмма для предельно простого случая;[1,1;1,0] Для этого ввода Matlab и октава дают мне 0,25 с использованием rcond и 1 / cond (x, 1), но в случае с использованием 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

1 Ответ

11 голосов
/ 12 февраля 2011

Я нашел ответ на свой вопрос.

Матрица должна быть разложена LU перед отправкой в ​​dgecon. Это кажется очень логичным, поскольку часто требуется решить систему после проверки условия, и в этом случае нет необходимости дважды разлагать матрицу. Та же идея относится и к норме, которая рассчитывается отдельно.

Следующий код - это все необходимые части, которые вычисляют номер взаимного условия с помощью 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;
}
...