## Matrix inverse - trouble with SGETRI call

Post here if you have a question about LAPACK or ScaLAPACK algorithm or data format

### Matrix inverse - trouble with SGETRI call

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.000000A[1][0]:3.000000        A[1][1]:4.000000Matrix A after sgetrfA[0][0]:2.000000        A[0][1]:0.500000A[1][0]:4.000000        A[1][1]:1.000000info for sgetrf_:0Parameter 3 to routine SGETRI was incorrectMac 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
susanunit1

Posts: 2
Joined: Thu Oct 28, 2010 6:34 pm
Location: San Diego

### Re: Matrix inverse - trouble with SGETRI call

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).
cottrell

Posts: 73
Joined: Thu Jan 15, 2009 1:40 pm

### Re: Matrix inverse - trouble with SGETRI call

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 sgetrfA[0]:3.000000   A[1]:0.333333   A[2]:4.000000   A[3]:0.666667info for sgetrf_:0Parameter 3 to routine SGETRI was incorrectMac OS BLAS parameter error in SGETRI, parameter #0, (unavailable), is 0sunit@susan-lindseys-macbook-pro ~/GLM-LAPACK \$ `

Any help is again appreciated!
susanunit1

Posts: 2
Joined: Thu Oct 28, 2010 6:34 pm
Location: San Diego

### Re: Matrix inverse - trouble with SGETRI call

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.
cottrell

Posts: 73
Joined: Thu Jan 15, 2009 1:40 pm

### Re: Matrix inverse - trouble with SGETRI call

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.
redgee34

Posts: 1
Joined: Mon Nov 24, 2014 7:19 am