//
// File:    ~/numerical/test_la.cpp
//
// Author:  Chuck Gartland (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.[so|a]
//     BLAS:    /usr/lib/libblas.[so|a]
// 
// Compiling:  There are several different ways to accomplish this.
//   One can either work from the C++ compiler (g++) or compile
//   the C++ driver (and any other C++ programs) separately
//   and then invoke the linker from Fortran (g77).  Either way,
//   one must include various libraries (mostly for I/O)
//   for the "other" language.
//
//   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).
//
//     $ g++ -c test_la.cpp
//     $ g77 test_la.o -llapack -lstdc++
//
//   or
//
//     $ g++ test_la.cpp -llapack -lg2c
//
//   See man pages for "c++/g++" and "f77/g77".
//
//   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 "C" { extern void sgetrf_(int *, int *, float (*)[5], int *,
                                 int [], int *) ; }

extern "C" { extern void sgetrs_(unsigned char *, int *, int *, float (*)[5],
                                 int *, int [], float [], int *, int *) ; }

int 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 solve

  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] ) ;

  return 1 ;
}

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