#include <math.h>
#include <stdio.h>

int main(int argc, char* argv[])
{
/***********************************************************************
 *                                                                     *
 * Kadin Tseng                                                         *
 * Scientific Computing and Visualization                              *
 * Boston University                                                   *
 * 1998 (modified 10/5/2006)                                           *
 *                                                                     *
 * This program solves                                                 * 
 *                                                                     *
 *  3x - 2y + 2z = 10                                                  *
 *   x + 2y - 3z = -1                                                  *
 *  4x +  y + 2z =  3                                                  *
 *                                                                     *
 * using the LAPACK routines sgemv, sgetrf and sgetrs.                 *
 * Since LAPACK routines are written in fortran, for applications in a *
 * C code, three key facts must be noted:                              *
 * 1) Fortran uses "passed by reference." Hence, EVERYTHING that is    *
 *    passed in the C calling program must be treated as pointers      *
 * 2) Fortran's arrays are column-based while C's are row-based        *
 *    Hence, C arrays must be transposed before passing into a LAPACK  *
 *    routine in order to have the correct matrix passed.              *
 * 3) All references to the LAPACK fortran routines must be appended   *
 *    with "_" on many systems (except IBM)                            *
 *                                                                     *
 ***********************************************************************/
      int i, info, n=3;
      int p, myid;
      int ipvt[n];
      float x[n], b[n];   /* solution and RHS */
      float A[n][n];      /* square matrix A */
      int one=1;
      float zerof=0.0, onef=1.0;
      char N='N';

      A[0][0] =  3.0; A[1][0] = -2.0; A[2][0] =  2.0; x[0] =  2.0;
      A[0][1] =  1.0; A[1][1] =  2.0; A[2][1] = -3.0; x[1] = -3.0;
      A[0][2] =  4.0; A[1][2] =  1.0; A[2][2] =  2.0; x[2] = -1.0;

      sgemv_(&N,&n,&n,&onef,A,&n,x,&one,&zerof,b,&one);  /* b = Ax */
      printf("The solution, x, and the RHS, b = Ax :\n");
      for (i=0; i<n; i++) {
        printf("%f %f\n", x[i], b[i]);
      }
      printf("\n");
      sgetrf_(&n,&n,A,&n,ipvt,&info);             /* factor A */
      sgetrs_(&N,&n,&one,A,&n,ipvt,b,&n,&info);     /* solve for x */

      printf("Solution, x, as obtained by sgetrf/sgetrs :\n");
      for (i=0; i<n; i++) {
        printf("%f ",b[i]);
      }
      printf("\n");
}