// 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
int_4 N = 5;
// scale of the value (if =1 values between -1 and 1)
r_8 scale = 1.;
// number of values change by +/- vbig
int_4 nbig = N;
r_8 vbig = 1000.;
// Nombre de ligne de matrice a imprimer
 int_4 nprline = 10;
//--------------------------------------------------------

//-- Decodage arguments
char c;
while((c = getopt(narg,arg,"hv:l:")) != -1) {
  switch (c) {
  case 'v' :
    sscanf(optarg,"%d,%lf,%d,%lf",&N,&scale,&nbig,&vbig);
    break;
  case 'l' :
    sscanf(optarg,"%d",&nprline);
    break;
  case 'h' :
    cout<<"tsttminv [-h] [-v N,scale,nbig,vbig] [-l nprline]"<<endl;
    cout<<"matrix filled with : {[-1,1] +/- vbig(nbig time)}*scale"<<endl;
    exit(-1);
  }
}
if(N<=1) N = 1;
cout<<"Taille matrice N = "<<N<<endl;
cout<<"Elements entre +/- "<<scale<<endl;
cout<<"Nombre de valeurs hors standard "<<nbig<<endl;
cout<<"Valeurs hors standard v = v*(+/- "<<vbig<<" )"<<endl;
cout<<"Nombre de lignes de matrice a imprimer "<<nprline<<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);
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;
int_4 i,j; double s;
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.; A[k] += (TYPE) s*vbig;}
A *= (TYPE) scale;
#if defined(COMPLEX)
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.; B[k] += (TYPE) s*vbig;}
B *= (TYPE) scale;
A += TYPE(0.,1.)*B;
#endif

//-- 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 vmin,vmax;
k = 0;
for(int i=0;i<N;i++) {
  double absv = ABS_VAL( AiA(i,i) );
  if(k==0) {vmin = vmax = absv; k++; continue;}
  if( absv > vmax ) vmax = absv;
  if( absv < vmin ) vmin = absv;
  k++;
}
cout<<"Limites sur la diagonale  "
    <<vmin<<" , "<<vmax<<"  n="<<k<<endl;
k = 0;
for(int i=0;i<N;i++) for(int j=0;j<N;j++) {
  if( i == j ) continue;
  double absv = ABS_VAL( AiA(i,j) );
  if(k==0) {vmin = vmax = absv; k++; continue;}
  if( absv > vmax ) vmax = absv;
  if( absv < vmin ) vmin = absv;
  k++;
}
cout<<"Limites hors diagonale  "
    <<vmin<<" , "<<vmax<<"  n="<<k<<endl;

exit(0);
}