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

Postby susanunit1 » Thu Nov 04, 2010 4:11 pm

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
susanunit1
 
Posts: 2
Joined: Thu Oct 28, 2010 6:34 pm
Location: San Diego

Re: Matrix inverse - trouble with SGETRI call

Postby cottrell » Mon Nov 08, 2010 2:05 pm

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: 68
Joined: Thu Jan 15, 2009 1:40 pm

Re: Matrix inverse - trouble with SGETRI call

Postby susanunit1 » Wed Nov 17, 2010 7:41 pm

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!
susanunit1
 
Posts: 2
Joined: Thu Oct 28, 2010 6:34 pm
Location: San Diego

Re: Matrix inverse - trouble with SGETRI call

Postby cottrell » Mon Nov 22, 2010 6:20 pm

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: 68
Joined: Thu Jan 15, 2009 1:40 pm

Re: Matrix inverse - trouble with SGETRI call

Postby redgee34 » Mon Nov 24, 2014 7:21 am

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


Return to Algorithm / Data

Who is online

Users browsing this forum: No registered users and 2 guests