/* solving the matrix equation A*x=b using LAPACK */
 
#include <stdio.h>
#define SIZE 3				/* dimension of matrix */

int sgesv_(int *n, int *nrhs, float *a, int *lda,
           int *ipiv, float *b, int *ldb, int *info);

int main(int argc, char **argv)
{
   int i, j , pivot[SIZE], ok;
   float A[SIZE][SIZE], b[SIZE], AT[SIZE*SIZE];	/* single precision!!! */


   A[0][0]=3.1;  A[0][1]=1.3;  A[0][2]=-5.7;	/* matrix A */
   A[1][0]=1.0;  A[1][1]=-6.9; A[1][2]=5.8;	
   A[2][0]=3.4;  A[2][1]=7.2;  A[2][2]=-8.8;	

   b[0]=-1.3;			/* if you define b as a matrix then you */
   b[1]=-0.1;			/* can solve multiple equations with */
   b[2]=1.8;			/* the same A but different b */ 	
 
   /* ******************************************************************* */
   /* IMPORTANT CONVERSION !!!						  */
   /* ******************************************************************* */
   for (i=0; i<SIZE; i++)	/* to call a Fortran routine from C we */
   {				/* have to transform the matrix */
      for(j=0; j<SIZE; j++) AT[j+SIZE*i]=A[j][i];		
   }						


   /* find solution using LAPACK routine SGESV, all the arguments have to */
   /* be pointers and you have to add an underscore to the routine name */

   sgesv_(3, 1, AT, 3, pivot, b, 3, &ok);      

   /*
      parameters in the order as they appear in the function call
      order of matrix A, number of right hand sides (b), matrix A,
      leading dimension of A, array that records pivoting, 
      result vector b on entry, x on exit, leading dimension of b
      return value 
   */ 
     						
   for (j=0; j<SIZE; j++) 
      printf("%f\n", b[j]);	/* print solution vector */
}