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


Return to Algorithm / Data

Who is online

Users browsing this forum: Google [Bot] and 2 guests

cron