Page 1 of 1

Matrix inverse - trouble with SGETRI call

PostPosted: Thu Nov 04, 2010 4:11 pm
by 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

Re: Matrix inverse - trouble with SGETRI call

PostPosted: Mon Nov 08, 2010 2:05 pm
by 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).

Re: Matrix inverse - trouble with SGETRI call

PostPosted: Wed Nov 17, 2010 7:41 pm
by 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!

Re: Matrix inverse - trouble with SGETRI call

PostPosted: Mon Nov 22, 2010 6:20 pm
by 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.

Re: Matrix inverse - trouble with SGETRI call

PostPosted: Mon Nov 24, 2014 7:21 am
by 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.