20 mai 2011

Calcularea inversei unei matrici

EXEMPLU

#include "lapack.h"

int main( int argc, char *argv[] )
{
int info = 0;
double *a;
int *ipiv;
int n;
int i,j;

if( argc < 2 )
{
fprintf(stderr, "%s N\n", argv[0]);
return 0;
}
n = atoi(argv[1]);

posix_memalign((void **)((void*)(&ipiv)), 128, sizeof(int)*n);
posix_memalign((void **)((void*)(&a)), 128, sizeof(double)*n*n);
if( a == NULL || ipiv == NULL )
{
fprintf(stderr, "a/ipiv malloc error\n");
return 0;
}

for( i = 0; i < n; i++ )
for( j = 0; j < n; j++ )
a[i+j*n]= (drand48() - 0.5f)*4;

for( i = 0; i < n; i++ )
for( j = 0; j < n; j++ )
printf("%2.1f ",a[i+j*n]);
printf("\n");

/*---------Call Cell LAPACK library---------*/
dgetrf_(&n, &n, a, &n, ipiv, &info);
if( info != 0 )
{
fprintf(stderr, "Call dgetrf error\n");
goto end;
}
printf("getrf_ \n");
for( i = 0; i < n; i++ )
{
for( j = 0; j < n; j++ )
printf("%2.1f ",a[i+j*n]);
printf("\n");
}

/*---------Query workspace-------*/
double workspace;
int tmp=-1;
int lwork;
double *work;
dgetri_(&n, a, &n, ipiv, &workspace, &tmp, &info);
lwork = (int)workspace;
work = malloc(sizeof(double)*lwork);
if(work == NULL)
{
printf("work malloc error\n");
goto end;
}

printf("getri_ (1)\n");
for( i = 0; i < n; i++ )
{
for( j = 0; j < n; j++ )
printf("%2.1f ",a[i+j*n]);
printf("\n");
}

/*---------Call Cell LAPACK library---------*/
dgetri_(&n, a, &n, ipiv, work, &lwork, &info);
if( info != 0 )
{
fprintf(stderr, "Call dgetri error\n");
free(work);
goto end;
}

printf("Inverse matrix completed!\n");

printf("getri_ (2)\n");
for( i = 0; i < n; i++ )
{
for( j = 0; j < n; j++ )
printf("%2.1f ",a[i+j*n]);
printf("\n");
}

end:
free(ipiv);
free(a);
return 0;
}


COMPILARE

gcc inverse.c -o inverse -llapack -lblas -lm

Niciun comentariu: