listas - Calcular el número de condición recíproca con lapack(es decir, rcond(x))

listas en racket

Question

Deseo hacer exactamente lo que hace en MATLAB / Octave usando LAPACK desde C. El manual de MATLAB me dice que se usa dgecon, y que usa una norma basada en 1.

Escribí un programa de prueba simple para un caso extremadamente simple; [1,1; 1,0] Para esta entrada, matlab y octava me dan 0.25 utilizando rcond y 1 / cond (x, 1), pero en el caso de usar LAPACK, este programa de ejemplo imprime 0.0. Para otros casos, como identidad, imprime el valor correcto.

Ya que MATLAB supuestamente está usando esta rutina con éxito, ¿qué estoy haciendo mal? Intento descifrar lo que hace Octave, con poco éxito ya que está envuelto en

#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;
}

Compilado con cc testdgecon.c -o testdgecon -llapack; ./testdgecon

Answer 1

Encontré la respuesta a mi propia pregunta.

La matriz debe descomponerse en LU antes de enviarse a dgecon. Esto parece muy lógico ya que a menudo uno quiere resolver el sistema después de verificar la condición, en cuyo caso no hay necesidad de descomponer la matriz dos veces. La misma idea se aplica a la norma que se calcula por separado.

El siguiente código es todas las partes necesarias para calcular el número de condición recíproca con 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;
}