Matrix inverse - trouble with SGETRI call

[susanunit1]

Hi all,

Thanks in advance for any and all help.

I'm attempting to get the inverse of a matrix using SGETRF, then SGETRI, but I'm running into an error that I can't decipher. My C code is simple, I initialize matrix A[2][2], then call sgetrf_, then sgetri_ but sgetri_ chokes on the third parameter.

Here's the code in entirety:

Code: Select all

#include <string.h>
#include <stdio.h>
#include <stdlib.h>

extern void sgetri_ ( int* n, const void* A, int* lda,
    int* ipiv, const void* work, int* lwork, int* info );

extern void sgetrf_ ( int* m, int* n, const void* a, int* lda,
    int* ipiv, int* info );

int main (void) {
    int i,j,n;
    float A[2][2]; //input matrix
    float alpha, beta;
    int lda, arows, acols;
    int ipiv[4];
    int lwork;
    int info;
    float work[100]; //dimension totally made up - what should this be?

    // silly matrix A is 2x2
    A[0][0]=1.0; A[0][1]=2.0; A[1][0]=3.0; A[1][1]=4.0;

    printf ("Matrix A before any calls:\n");
    for (i=0;i<2;i++) {
        for (j=0;j<2;j++)
            printf ("A[%d][%d]:%f\t", i,j, A[i][j]);
        printf ("\n");
    }
    printf ("\n");

    arows=acols=2;
    lda=2;
    alpha = 1.0; beta=1.0;
    memset (ipiv, 0, sizeof(int)*4);
    memset (work, 0, sizeof(float)*100);

    sgetrf_ ( &arows, &acols, A, &lda, ipiv, &info );
    printf ("Matrix A after sgetrf\n");
    for (i=0;i<2;i++) {
        for (j=0;j<2;j++)
            printf ("A[%d][%d]:%f\t", i,j, A[i][j]);
        printf ("\n");
    }
    printf ("info for sgetrf_:%d\n\n", info);

    n=acols*arows;
    lda=2;
    lwork=100; //magic number - dimension of work

    sgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);
    printf ("Matrix A after sgetri\n");
    for (i=0;i<2;i++) {
        for (j=0;j<2;j++)
            printf ("A[%d][%d]:%f\t", i,j, A[i][j]);
        printf ("\n");
    }
    printf ("info for sgetri_:%d\n\n", info);

    return(0);
}

I'm compiling this as so:

Code: Select all

/usr/bin/gcc -c -g -Wall  inverse.c
/usr/bin/gcc -o inverse inverse.o -framework vecLib

and the output is this:

Code: Select all

Matrix A before any calls:
A[0][0]:1.000000        A[0][1]:2.000000
A[1][0]:3.000000        A[1][1]:4.000000

Matrix A after sgetrf
A[0][0]:2.000000        A[0][1]:0.500000
A[1][0]:4.000000        A[1][1]:1.000000
info for sgetrf_:0

Parameter 3 to routine SGETRI was incorrect
Mac OS BLAS parameter error in SGETRI, parameter #0, (unavailable), is 0

What the heck am I doing wrong? Why is parameter 3, lda, incorrect?

Susan

[cottrell]

You're defining your A matrix as a two-dimensional array:

Code: Select all

float A[2][2]; //input matrix


This won't work for Lapack: you need a one dimensional
array (with the elements in column-major order). For example,

Code: Select all

float A[] = {1, 3, 2, 4};


for a 2x2 matrix with first row (1, 2) and second row (3, 4).

[susanunit1]

Hi, thanks very much for your reply and sorry for my late response (on vacation).

I'm confused as to why I need a one dimensional array. The man pages for sgetri_ & sgetrf_ both say that the input matrix A is:

REAL A( LDA, * )

which looks pretty two-dimensional to me. (although Fortran is not my first language). Anyway I tried it your way and I'm still getting the same error (AAARGH)

Here's the newish code:

Code: Select all

#include <stdlib.h>

extern void sgetri_ ( int* n, const void* A, int* lda,
    int* ipiv, const void* work, int* lwork, int* info );

extern void sgetrf_ ( int* m, int* n, const void* a, int* lda,
    int* ipiv, int* info );

float A[] = {1.0,3.0,2.0,4.0};

int main (void) {
    int i,n;
    float alpha, beta;
    int lda, arows, acols;
    int ipiv[4];
    int lwork;
    int info;
    float work[100];

    arows=acols=2;
    lda=2;
    alpha = 1.0; beta=1.0;
    memset (ipiv, 0, sizeof(int)*4);
    memset (work, 0, sizeof(float)*100);

    sgetrf_ ( &arows, &acols, A, &lda, ipiv, &info );
    printf ("Matrix A after sgetrf\n");
    for (i=0;i<4;i++) {
        printf ("A[%d]:%f\t", i,A[i]);
    }
    printf ("\ninfo for sgetrf_:%d\n\n", info);

    n=acols*arows;
    lda=2;
    lwork=100; //magic number - dimension of work
   
    sgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);
    printf ("Matrix A after sgetri\n");
    for (i=0;i<4;i++) {
            printf ("A[%d]:%f\t", i,A[i]);
    }
    printf ("\ninfo for sgetri_:%d\n\n", info);

    return(0);
}

Here's how I compiled it:

Code: Select all

/usr/bin/gcc -c -g -Wall  2inverse.c
/usr/bin/gcc -o 2inverse 2inverse.o -framework vecLib

and here's the output:

Code: Select all

Matrix A after sgetrf
A[0]:3.000000   A[1]:0.333333   A[2]:4.000000   A[3]:0.666667
info for sgetrf_:0

Parameter 3 to routine SGETRI was incorrect
Mac OS BLAS parameter error in SGETRI, parameter #0, (unavailable), is 0
sunit@susan-lindseys-macbook-pro ~/GLM-LAPACK $

Any help is again appreciated!

[cottrell]

Code: Select all

    n=acols*arows;
    lda=2;
    lwork=100; //magic number - dimension of work

    sgetri_ (&n, A, &lda, ipiv, work, &lwork, &info);


Your n value is wrong. It's supposed to be the order of the
matrix, which is 2 in your example, not the number of
elements in the matrix. Lapack is choking on the lda
argument, because if the order of the input matrix were
really 4 then an lda of 2 would not be valid.

[redgee34]

You are supplying LWORK = -1. This means that you are doing workspace query. The function just writes recommended length of WORK in the first element of YM and immediately exits. For the workspace query case the result of unchanged matrices is expected. For your case the minimal length of YM and value of LWORK is 4.