/*
 * File:    ~/numerical/test_la.c
 *
 * Author:  Chuck Gartland (6-99, 10-01, 10-05)
 *
 * C test program for using LAPACK on the Linux machines on the
 * Kent Math/CS Network.
 *
 * Library locations:
 *
 *     LAPACK:  /usr/lib/liblapack.[a|so|...]
 *     BLAS:    /usr/lib/libblas.[a|so|...]
 *
 * Compiling:  There are several different ways to accomplish this,
 *   but I believe the easiest path is to compile the C driver
 *   (and any other C programs) separately and then invoke the
 *   linker/loader from Fortran (in order to get the proper Fortran
 *   runtime libraries and load path).
 *
 *   Note the appended underscore in the names of the LAPACK routines 
 *   in the source code below: sgetrf_() (vs sgetrf() from Fortran) and
 *   sgetrs_() (vs sgetrs() from Fortran).
 *
 *     $ cc -c test_la.c
 *     $ f77 test_la.o -llapack
 *
 *   See man pages for "cc" and "f77".
 *
 *   NOTE - In calling a Fortran sub-program from C, one must
 *          be attentive to certain language differences.
 *          In particular, Fortran stores arrays by columns
 *          ("column major" order); whereas C stores them by rows
 *          ("row major" order).  Also, Fortran passes all arguments
 *          by reference (pointer).
 *
 *          To see the difference, compare this program with
 *          the companion Fortran version "test_la.f".
 *
 * Test prob: [ 1 2 3 ][x]   [ 6] ,  w/ true sol'n  [1]
 *            [ 4 5 6 ][y] = [15]                   [1]
 *            [ 7 8 0 ][z]   [15]                   [1]
 *
 */

#include <stdio.h>

float A[4][5] = { { 1., 2., 3., 0., 0. },          /* deliberately chosen    */
                  { 4., 5., 6., 0., 0. },          /* larger dimensions      */
                  { 7., 8., 0., 0., 0. },          /* to illustrate the role */
                  { 0., 0., 0., 0., 0. } } ;       /* of 'lda'               */

float b[3] = { 6., 15., 15. } ;

int ipiv[3] ;

extern void sgetrf_() ;

extern void sgetrs_() ;

main()
{
  int M = 3, N = 3, LDA = 5, LDB = 3, NRHS = 1, INFO ;

  unsigned char TRANS = 'T' ;

  sgetrf_( &M, &N, A, &LDA, ipiv, &INFO ) ;                     /* LU factor */

  printf( " info = %2d \n", INFO ) ;

  printf( " ipiv(1) = %2d \n", ipiv[0] ) ;
  printf( " ipiv(2) = %2d \n", ipiv[1] ) ;
  printf( " ipiv(3) = %2d \n", ipiv[2] ) ;

  sgetrs_( &TRANS, &N, &NRHS, A, &LDA, ipiv, b, &LDB, &INFO ) ;  /* back slv */

  printf( " info = %2d \n", INFO ) ;

  printf( " b(1) = %.6e \n", b[0] ) ;
  printf( " b(2) = %.6e \n", b[1] ) ;
  printf( " b(3) = %.6e \n", b[2] ) ;
}

/*
 * Appended sample output:
 *
 *  info =  0 
 *  ipiv(1) =  3 
 *  ipiv(2) =  3 
 *  ipiv(3) =  3 
 *  info =  0 
 *  b(1) = 1.000000e+00 
 *  b(2) = 9.999999e-01 
 *  b(3) = 1.000000e+00 
 *
 */