#include "machdefs.h"

#include <math.h>
#include <iostream.h>

#include "nbrandom.h"
#include "tarrinit.h"
#include "tarray.h"
#include "tvector.h"
#include "timing.h"
#include "intflapack.h"


void lpk_tarr(int n);
void lpk_tmtx(int n);

int main(int narg, char* arg[])
{

  SophyaInit();
  InitTim();   // Initializing the CPU timer



  if (narg < 2) {
    cout << " lpk - LinAlg/LapackServer test - Usage lpk t/s [size=5] [prtlev=0] [nprtmax=100]\n" 
	 << " t:  lpk_tarr() rztest_lapack with TArray<r_4> \n" 
	 << " s:  lpk_tmtx() LapackServer with TMatrix<r_8> \n" << endl;
    exit(0);
  }
  int n = 5;
  int opt = 's';
  opt = *(arg[1]);
  if (narg > 2) n = atoi(arg[2]);
  int nprt = 100;
  int prtlev = 1;
  if (narg > 3) prtlev = atoi(arg[3]);
  if (narg > 4) nprt = atoi(arg[4]);

  BaseArray::SetMaxPrint(nprt, prtlev);
  try {
    if (opt == 's') lpk_tmtx(n);
    else lpk_tarr(n);
  }
  catch (PThrowable exc) {
    cerr << " catched Exception (lpk.cc) " << exc.Msg() << endl;
  }  
  catch (...) {
    cerr << " catched unknown (...) exception (lpk.cc) " << endl; 
  } 
  
  PrtTim(" End of lpk LinAlg/Lapack test ");
  cout << " ---------------  END of Programme -------------- " << endl;
}


void lpk_tmtx(int n)
{
  int i,j,k;
  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
  Matrix a(n,n);
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)  a(j,i) = GauRnd(0., 1.);
  cout << " ------------ Matrix A = \n " << a << "\n" << endl; 
  
  Vector x(n), b;
  //  Matrix  x(n,1), b;
  cout << " ------------ Vector X = \n " << x << "\n" << endl; 
  for(i=0; i<n; i++) x(i) = GauRnd(2., 1.5);
  b = a*x;

  //  cout << ":::::::: rztest_lapack - Size=" << n << "  ::::::::: " << endl;
  cout << " ------- lpk_tmtx() LapackServerTest Using TMatrix<r_8> -------- " << endl; 
  cout << " ------------ Matrix A = \n " << a << "\n" << endl; 
  cout << " ------------ Matrix X = \n " << x << "\n" << endl; 
  cout << " ------------ Matrix B = \n " << b << "\n" << endl; 

  cout << "\n   Calling LapackLinSolve(a,b) .... " << endl;
  LapackLinSolve(a,b);

  cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl; 
  Vector diff = b-x;
  cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl; 

}

void lpk_tarr(int n)
{
  int i,j,k;
  TArray<r_4> a(n,n);
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)  a(i,j,0) = GauRnd(0., 1.);
  
  TArray<r_4> x(n,1), b(n,1);
  r_4 sum ;
  for(i=0; i<n; i++) x(i,0,0) = GauRnd(2., 1.5);
  for(i=0; i<n; i++) {
    sum = 0.;
    for(j=0; j<n; j++) sum +=  a(i,j,0)*x(j,0,0);
    b(i,0,0) = sum;
  }

  //  cout << ":::::::: rztest_lapack - Size=" << n << "  ::::::::: " << endl;
  cout << " ------- lpk_tarr() LapaackTest Using TArray<r_4> -------- " << endl; 
  cout << " ------------ Array A = \n " << a << "\n" << endl; 
  cout << " ------------ Array X = \n " << x << "\n" << endl; 
  cout << " ------------ Array B = \n " << b << "\n" << endl; 

  cout << "\n     Calling rztest_lapack ... " << endl;

  rztest_lapack(a, b);

  cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl; 
}