// Test de l'inversion de matrice (cmv 14/04/00)
#include "machdefs.h"
#include <iostream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "ntoolsinit.h"
#include "pexceptions.h"
#include "array.h"
#include "srandgen.h"

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

#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;
//--------------------------------------------------------

//-- Decodage arguments
char c;
while((c = getopt(narg,arg,"hv:l:a:")) != -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 'h' :
    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline] [-a nalea]"<<endl;
    cout<<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<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<<endl;

//-- Initialization arrays
SophyaInit();

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

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

//-- Fill matrices
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;

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

//-- Inversion
cout<<"------------ Inversion"<<endl;
InvA = Inverse(A);
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
double 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;

exit(0);
}