#include "sopnamsp.h"
#include "machdefs.h"

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

#include "srandgen.h"
#include "tarrinit.h"
#include "tvector.h"
#include "sopemtx.h"
#include "timing.h"

/********************************************************
  Programme test pour les fonctionalites de calcul
  sur matrices et vecteurs dans SOPHYA/TArray
              Reza - fev 2005 
*********************************************************/ 

int mxv_inverse(int l);                    // ErrorCode = 64

static double TOLERANCE = 1.e-6;
static int nprt = 100;
static int prtlev = 1;

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

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



  if (narg < 2) {
    cout << " tmxv : Test SOPHYA Matrix/Vector operations\n" 
	 << " Usage: tmxv select [sizeL,C=5,5] [prtlev=1] [nprtmax=100]\n" 
	 << "   select= inverse / all= \n"
	 << "   inverse:  ComputeInverse() r_8 \n" << endl;
    exit(0);
  }
  int l,c;
  l = c = 5;
  int wsf = 2;
  string opt = arg[1];
  if (narg > 2) sscanf(arg[2], "%d,%d", &l, &c); 
  if (narg > 3) prtlev = atoi(arg[3]);
  if (narg > 4) nprt = atoi(arg[4]);

  cout << " tmxv - TMatrix/TVector operations tests sizeL,C=" << l << "," << c 
       << " opt= " << opt << endl;
  int rc = 0;
  BaseArray::SetMaxPrint(nprt, prtlev);
  try {
    if (opt == "inverse") rc = mxv_inverse(l);
    else { cout << "tmxv/ Unknown option " << opt << " ! " << endl; rc = 66; }
  }
  catch (PThrowable exc) {
    cerr << " tmxv catched Exception (tmxv.cc) " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (...) {
    cerr << " tmxv catched unknown (...) exception (tmxv.cc) " << endl; 
    rc = 78; 
  } 
  
  PrtTim(" End of tmxv matrix inversion test ");
  cout << " ---------------  END of Programme -------- (Rc= " 
       << rc << ") --- " << endl;
  return(rc);
}

// -----------------------------------------------------------------------
/* Nouvelle-Fonction */
int mxv_inverse(int m)
{
  cout << " mxv_inverse() - Test SimpleMatrixOperation<r_8> " << endl;
  
  Matrix a(m , m), ainv, idmx(m,m), mxprod, diff;
  a = RandomSequence(RandomSequence::Gaussian, 0., 4.);
  PrtTim(" End of a_inint");

  cout << "\n   Calling  SimpleMatrixOperation::Inverse .... " << endl;
  SimpleMatrixOperation<r_8> smo;
  ainv = smo.Inverse(a);
  PrtTim(" End of smo.Inverse()");
  mxprod = a*ainv;
  PrtTim(" End of mxprod = a*ainv");
  idmx = IdentityMatrix();
  PrtTim(" End of idmx = IdentityMatrix()");
  diff = mxprod-idmx;
  PrtTim(" End of Compute(diff)");

  if (prtlev > 0) {
    cout << " ------------ Matrix A = \n " << a << "\n" << endl; 
    cout << " ------------ Matrix Inverse(A) = \n " << ainv << "\n" << endl; 
    cout << " ------------ Matrix mxprod = A*Inverse(A) = \n " << mxprod << "\n" << endl; 
  }
  double min,max;
  diff.MinMax(min, max);
  PrtTim(" End of diffCheck");
  cout << " Min/Max difference Matrix (?=0) , Min= " << min 
       << " Max= " << max << endl;
  if ((fabs(min) > TOLERANCE) || (fabs(max) > TOLERANCE)) {
    cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 
	 << endl;
    return 64;
  }
  return 0;

}