// Test de l'inversion de matrice (cmv 14/04/00)
// Pb numeriques: dapax22 et daplxa115
//                tsttminv -v 55,100,130,1000000000 -l 0 -a 1 

// if defined COMPLEX , if not REAL
//#define COMPLEX
// if defined 32 bits precision, if not 64 bits
// #define PRECIS32
#define ALSO_LAPACK

#include "machdefs.h"
#include <iostream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "timing.h"
#include "ntoolsinit.h"
#include "pexceptions.h"
#include "array.h"
#include "srandgen.h"

#ifdef ALSO_LAPACK
#include "intflapack.h"
#endif


#if defined(COMPLEX)
  #define ABS_VAL(_x_) sqrt((double)((_x_).real()*(_x_).real() + (_x_).imag()*(_x_).imag()))
  #if defined(PRECIS32)
    #define TYPE complex<r_4>
  #else
    #define TYPE complex<r_8>
  #endif
#else
  #define ABS_VAL(_x_) fabs((double)_x_)
  #if defined(PRECIS32)
    #define TYPE r_4
  #else
    #define TYPE r_8
  #endif
#endif

int main(int narg,char *arg[])
{
//--------------------------------------------------------
// number of lines/columns
uint_4 N = 5;
// scale of the value (if =1 values between -1 and 1)
r_8 scale = 1.;
// number of values change by +/- vbig
uint_4 nbig = N;
r_8 vbig = 1000.;
// Nombre de lignes de matrice a imprimer
uint_4 nprline = 10;
// Initialisation du pauvre de l'aleatoire
uint_4 nalea = 0;
// data scaling for gauss pivoting and determinant
int tscal = 1;
bool detok=false;
bool timok = false;
//--------------------------------------------------------

//-- Decodage arguments
char c;
while((c = getopt(narg,arg,"hdtv:l:a:n:")) != -1) {
  switch (c) {
  case 'v' :
    sscanf(optarg,"%d,%lf,%d,%lf",&N,&scale,&nbig,&vbig);
    break;
  case 'l' :
    sscanf(optarg,"%d",&nprline);
    break;
  case 'a' :
    sscanf(optarg,"%d",&nalea);
    break;
  case 'n' :
    sscanf(optarg,"%d",&tscal);
    break;
  case 'd' :
    detok = true;
    break;
  case 't' :
    timok = true;
    break;
  case 'h' :
    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline] [-a nalea]"<<endl
        <<"         [-n typs] [-d] [-t]"<<endl
        <<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<endl
	<<"-n 0/1/2 : data scaling 0=no, 1=global (def), 2=row-by-row"<<endl
	<<"-d : also compute determinant"<<endl
	<<"-t : do timing"<<endl;
    exit(-1);
  }
}
if(N<=1) N = 1;
cout<<"Taille matrice NxN, N = "<<N<<endl;
cout<<"Elements entre +/- scale = "<<scale<<endl;
cout<<"Nombre de valeurs hors standard nbig = "<<nbig<<endl;
cout<<"Valeurs hors standard (+/- vbig = "<<vbig<<" ) * scale"<<endl;
cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<endl;
cout<<"Initialisation de l aleatoire par "<<nalea<<" tirages"<<endl;
 cout<<"Data scaling "<<tscal<<" determinant="<<detok<<" timing="<<timok<<endl;
cout<<endl;

//-- Initialization arrays
SophyaInit();
if(timok) InitTim();
#ifdef ALSO_LAPACK
  BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
#endif

//-- Definition arrays
TMatrix< TYPE > A(N,N);
TMatrix< TYPE > InvA(N,N);
TMatrix< TYPE > AiA(N,N);
TMatrix< TYPE > B(N,N);
TMatrix< TYPE > Adum(N,N);
TVector< int >  Vind(N*N);
A.Show();

SimpleMatrixOperation< TYPE >::SetGausPivScal(tscal);

//-- Mise a zero
A    = (TYPE) 0;
InvA = (TYPE) 0;
AiA  = (TYPE) 0;
B    = (TYPE) 0;
Adum = (TYPE) 0;
BaseArray::SetMaxPrint(nprline*N,0);

//-- Fill matrices
double vmax=-1.;
uint_8 k;
uint_4 i,j; double s;
uint_4 nbr=0, nbi=0;
Vind = 0;
if(nalea>0) for(i=0;i<nalea;i++) drand01();
A = Sequence(RandomSequence(RandomSequence::Flat,0.,1.));
if(nbig>0) for(i=0;i<nbig;i++) {
  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
  s=(drand01()>0.5)?1.:-1.;
  if(Vind[k]==0) {A[k] += (TYPE) s*vbig; Vind[k]=1; nbr++;}
}
A *= (TYPE) scale;
#if defined(COMPLEX)
Vind = 0;
B = Sequence(RandomSequence(RandomSequence::Flat,0.,1.));
if(nbig>0) for(i=0;i<nbig;i++) {
  k=(uint_8)(drand01()*N*N); if(k>=N*N) k--;
  s=(drand01()>0.5)?1.:-1.;
  if(Vind[k]==0) {B[k] += (TYPE) s*vbig; Vind[k]=1; nbi++;}
}
B *= (TYPE) scale;
A += TYPE(0.,1.)*B;
#endif
cout<<"Nombre de valeurs BIG reelles = "<<nbr<<", imaginaires = "<<nbi<<endl;
Adum = A;

//-- Print matrice A
cout<<"------------ TMatrix A :"<<endl;
if(nprline>0) {cout<<A; cout<<endl<<endl;}

//-- Inversion
cout<<"------------ Inversion"<<endl;
//InvA = Inverse(A);
InvA = IdentityMatrix(1.,N);
if(timok) PrtTim("--- End of filling ---");
TYPE det = GausPiv(A,InvA,detok);
if(timok) PrtTim("--- End of GausPiv ---");
cout<<"Det = "<<det<<endl;
cout<<"------------ TMatrix InvA = A^(-1):"<<endl;
if(nprline>0) {cout<<InvA; cout<<endl<<endl;}

//-- AiA = A * InvA
cout<<"AiA = A * InvA"<<endl;
AiA = A * InvA;
cout<<"------------ TMatrix AiA = A * InvA:"<<endl;
if(nprline>0) {cout<<AiA; cout<<endl;}

//-- Check
vmax=-1.;
for(i=0;i<N;i++) {
  double absv = ABS_VAL( 1. - AiA(i,i) );
  if( absv > vmax ) vmax = absv;
}
cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
vmax = -1.;
for(i=0;i<N;i++) for(j=0;j<N;j++) {
  if( i == j ) continue;
  double absv = ABS_VAL( AiA(i,j) );
  if( absv > vmax ) vmax = absv;
}
cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;

////////////////////////////////////
///////// Test avec Lapack /////////
////////////////////////////////////
#ifdef ALSO_LAPACK

cout<<endl<<"LAPACK Cross-Check"<<endl;
InvA = IdentityMatrix(1.,N);
if(timok) PrtTim("--- End of check ---");
LapackLinSolve(Adum,InvA);
if(timok) PrtTim("--- End of LapackLinSolve ---");
AiA = A * InvA;
vmax=-1.;
for(i=0;i<N;i++) {
  double absv = ABS_VAL( 1. - AiA(i,i) );
  if( absv > vmax ) vmax = absv;
}
cout<<"Ecart maximum par rapport a 1 sur diagonale = "<<vmax<<endl;
vmax = -1.;
for(i=0;i<N;i++) for(j=0;j<N;j++) {
  if( i == j ) continue;
  double absv = ABS_VAL( AiA(i,j) );
  if( absv > vmax ) vmax = absv;
}
cout<<"Ecart maximum par rapport a 0 hors diagonale = "<<vmax<<endl;

#endif

if(timok) PrtTim("--- End of Job ---");
exit(0);
}