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

////////////////////////////////////////////////////////////////////////////////////////
int main(int narg,char *arg[])
{
SophyaInit();
r_8 v[3][4] = {{1,2,3,4},{4,5,6,5},{7,8,9,10}};

{

TMatrix<r_8> A(2,3);
{for(int i=0;i<2;i++) for(int j=0;j<3;j++) A(i,j) = v[i][j];}
cout<<"A:"<<endl<<A;

{
cout << " Test ecriture PPersist TMatrix " << endl; 
FIO_TArray<r_8> ftm(A);
ftm.Write("tmtx.ppf");
}
{
cout << " Test lecture PPersist TMatrix " << endl; 
FIO_TArray<r_8> ftm("tmtx.ppf");
cout<<"Alue:"<<(TMatrix<r_8>)ftm<<endl;
}

cout<<"Matrices B(2,3)"<<endl;
TMatrix<r_8> B(3,2);
{for(int i=0;i<3;i++) for(int j=0;j<2;j++) B(i,j) = v[i][j];}
cout<<"B:"<<endl<<B;

cout<<"Matrices C(3,3)"<<endl;
TMatrix<r_8> C(3,3);
{for(int i=0;i<3;i++) for(int j=0;j<3;j++) C(i,j) = v[i][j];}
cout<<"C:"<<endl<<C;

cout<<"Matrices D(2,2)"<<endl;
TMatrix<r_8> D(2,2);
{for(int i=0;i<2;i++) for(int j=0;j<2;j++) D(i,j) = v[i][j];}
cout<<"D:"<<endl<<D;

cout<<"Matrices E(3,4)"<<endl;
TMatrix<r_8> E(3,4);
{for(int i=0;i<3;i++) for(int j=0;j<4;j++) E(i,j) = v[i][j];}
cout<<"E:"<<endl<<E;

//-----------------------------------------------
cout<<endl<<"Matrices AA(A)"<<endl;
TMatrix<r_8> AA(A);
cout<<"AA:"<<endl<<AA;

//-----------------------------------------------
cout<<endl<<"Matrices AA2; AA2=A;"<<endl;
TMatrix<r_8> AA2; AA2=A;
cout<<"AA2:"<<endl<<AA2;

//-----------------------------------------------
cout<<endl<<"Matrices AAA.Clone(A)"<<endl;
TMatrix<r_8> AAA; AAA.Clone(A);
cout<<"AAA:"<<endl<<AAA;
cout<<endl<<"Matrices BBB.Clone(B)"<<endl;
TMatrix<r_8> BBB; BBB.Clone(B);
cout<<"BBB:"<<endl<<BBB;
cout<<"Matrices CCC.Clone(C)"<<endl;
TMatrix<r_8> CCC; CCC.Clone(C);
cout<<"CCC:"<<endl<<CCC;

//-----------------------------------------------
cout<<endl<<"Matrices I(5,5)=1.234"<<endl;
TMatrix<r_8> I(5,5); I = 5;
cout<<"I:"<<endl<<I;

//-----------------------------------------------
cout<<endl<<"Matrices AAA+=1"<<endl;
AAA += 1;
cout<<"AAA:"<<endl<<AAA;
cout<<"Matrices AAA-=1"<<endl;
AAA -= 1;
cout<<"AAA:"<<endl<<AAA;
cout<<"Matrices AAA*=10"<<endl;
AAA *= 10;
cout<<"AAA:"<<endl<<AAA;
cout<<"Matrices AAA/=10"<<endl;
AAA /= 10;
cout<<"AAA:"<<endl<<AAA;

//-----------------------------------------------
cout<<endl<<"Matrices AAA+=A"<<endl;
AAA += A;
cout<<"AAA:"<<endl<<AAA;
cout<<"Matrices AAA-=A"<<endl;
AAA -= A;
cout<<"AAA:"<<endl<<AAA;

//-----------------------------------------------
cout<<endl<<"Matrices AAA=A"<<endl;
AAA = A;
cout<<"AAA:"<<endl<<AAA;
cout<<endl<<"Matrices BBB=B"<<endl;
BBB = B;
cout<<"BBB:"<<endl<<BBB;

//---------------------------------------------
TMatrix<r_8> a(A,false);
cout<<endl<<"Matrices AAA = AAA * C"<<endl;
{for(uint_4 i=0;i<AAA.NRows();i++) for(uint_4 j=0;j<C.NCols();j++)
  {a(i,j) = 0; for(uint_4 k=0;k<AAA.NCols();k++) a(i,j) += AAA(i,k)*C(k,j);}}
AAA = AAA * C;
cout<<"AAA:"<<endl<<AAA;
cout<<"a:"<<endl<<a;
TMatrix<r_8> b(B,false);
cout<<"Matrices BBB = BBB * D"<<endl;
{for(uint_4 i=0;i<BBB.NRows();i++) for(uint_4 j=0;j<D.NCols();j++)
  {b(i,j) = 0; for(uint_4 k=0;k<BBB.NCols();k++) b(i,j) += BBB(i,k)*D(k,j);}}
BBB = BBB * D;
cout<<"B:"<<endl<<BBB;
cout<<"b:"<<endl<<b;

//-----------------------------------------------
AA.Clone(A);
cout<<endl<<"Matrices AA = AA + 10"<<endl;
AA = AA + 10.; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = AA - 10"<<endl;
AA = AA - 10.; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = 10 - AA"<<endl;
AA = 10. - AA; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = AA - 10"<<endl;
AA = AA - 10.; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = -1 * AA"<<endl;
AA = -1. * AA; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = 10 * AA"<<endl;
AA = 10. * AA; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = AA * 10"<<endl;
AA = AA * 10.; cout<<"AA"<<endl<<AA;
cout<<endl<<"Matrices AA = AA / 100"<<endl;
AA = AA / 100.; cout<<"AA"<<endl<<AA;

//-----------------------------------------------
TMatrix<r_8> R;
AA.Clone(A); AA += 10.; AAA.Clone(A); AAA += 100.;
cout<<endl<<"Matrices R = A + AA"<<endl;
R = A + AA; cout<<"R"<<endl<<R;
cout<<endl<<"Matrices R = AA - A"<<endl;
R = AA - A; cout<<"R"<<endl<<R;
cout<<endl<<"Matrices R = A + AA + AAA"<<endl;
R = A + AA + AAA; cout<<"R"<<endl<<R;
cout<<endl<<"Matrices R = A + 100 + 2*A"<<endl;
R = A + 100. + 2.*A; cout<<"R"<<endl<<R;

//-----------------------------------------------
TMatrix<r_8> RR(2,4);
{for(uint_4 i=0;i<A.NRows();i++) for(uint_4 j=0;j<E.NCols();j++)
  {RR(i,j) = 0; for(uint_4 k=0;k<A.NCols();k++) RR(i,j) += A(i,k)*E(k,j);}}
cout<<endl<<"Matrices R = A * E"<<endl;
R = A * E;
cout<<"R"<<endl<<R; cout<<"RR"<<endl<<RR;

//-----------------------------------------------
cout<<endl<<"Matrices A = A.Transpose 2x"<<endl;
cout<<"A"<<endl<<A;
A = A.Transpose(); cout<<"trA"<<endl<<A;
A = A.Transpose(); cout<<"trtrA"<<endl<<A;
cout<<endl<<"Matrices AT = A.Transpose()  2x"<<endl;
TMatrix<r_8> AT;
AT = A.Transpose(); cout<<"AT"<<endl<<AT;
AT = AT.Transpose(); cout<<"trAT"<<endl<<AT;

//-----------------------------------------------
cout<<endl<<"Test Inversion avec couplage a M=Matrix(C)"<<endl;
TMatrix<r_8> X(10,10);
{for(int i=0;i<10;i++) for(int j=0;j<10;j++)
  {X(i,j) = drand01(); if(drand01()>0.8) X(i,j) += 1.e+15*drand01();}}
cout<<"TMatrix X"<<endl<<X;
{
Matrix M(X); cout<<"Matrix M"<<endl<<M<<endl;
Matrix MI = Inverse(M); cout<<"MI"<<endl<<MI;
MI = MI * M; cout<<"MI*M"<<endl<<MI;
r_8 xmin=-1.;
for(int i=0;i<10;i++) for(int j=0;j<10;j++)
   if(i!=j && fabs(MI(i,j))>xmin) xmin=fabs(MI(i,j));
cout<<"Biggest off diagonal identity matrix element is "<<xmin<<endl;
} // destruction de M,MI
cout<<"TMatrix X"<<endl<<X;

//-----------------------------------------------
{
cout<<endl<<"Test Inversion directe"<<endl;
cout<<"TMatrix X"<<endl<<X;
TMatrix<r_8> Xinv;
Xinv = Inverse(X);
cout<<"TMatrix Xinv"<<endl<<Xinv;
Xinv = Xinv * X;
cout<<"Xinv = Xinv * X"<<endl<<Xinv;
r_8 xmin=-1.;
for(int i=0;i<10;i++) for(int j=0;j<10;j++)
   if(i!=j && fabs(Xinv(i,j))>xmin) xmin=fabs(Xinv(i,j));
cout<<"Biggest off diagonal identity matrix element is "<<xmin<<endl;
cout<<"TMatrix X"<<endl<<X;
}

//-----------------------------------------------
cout<<endl<<"Test liens TMatrix L avec Matrix MI"<<endl;
Matrix MI(10,5);
{for(int i=0;i<10;i++) for(int j=0;j<5;j++) MI(i,j) = 1000*i+j;}
cout<<"MI"<<endl<<MI;
{
  sa_size_t siz[5];
  siz[0] = MI.NCols();
  siz[1] = MI.NRows();
  TArray<r_8> A(2, siz, MI.Data(), 1, 0, new Bridge);
  TMatrix<r_8> L(A);
  cout<<"L"<<endl<<L;
  L *= 100.; cout<<"L*=100."<<endl<<L;
} // destruction de L
cout<<"MI(4,4)="<<MI(4,4)<<endl;

} // destruction de toutes les matrices
exit(0);
}