load "lapack"
// test of lapack load file ... 
// .......................
// load "fflapack" obsolete (F. Hecht version 3.8)
// use  load "lapack"  
int nerr=0; // nomber of err .. 
//  to set a full matrix .. 
macro SETM(A,i,j,fij) 
{ 
    for(int i=0;i<A.n;++i) 
      for(int j=0;j<A.m;++j) 
        A(i,j)= (fij) ;
}//EOM 

// The of Eigen Problem ..
NewMacro TEST(Type,Typevp,vtype,ComputeEV)
{     
cout << "\n\n *****************   VP \n\n\n"; 
int n=5;
Type[int,int] A(n,n),A1(n,n),B(n,n),Id(n,n);
SETM(A,i,j,(i==j) ? n+1 : 1); 
SETM(Id,i,j,real(i==j)); 
A(0,n-1)=vtype;

cout << A << endl;
A1=A^-1;
cout << A1 << endl;

Typevp[int] vp(n);
Typevp[int,int] VP(n,n),KK(n,n);

int nn= ComputeEV(A,vp,VP);
 cout << " vp = " << vp << endl;
 cout << " VP = " << VP << endl;
 
 // verification ... 
 KK =0.;
 for(int i=0;i<n;++i)
   for(int j=0;j<n;++j)
     for(int k=0;k<n;++k)
       //KK(i,j) += (A(i,k) - vp[j]* real(i==k) ) *VP(k,j);
       KK(i,j) += (A(i,k) - ((i==k)?vp[j]:0.0) ) *VP(k,j);
 cout <<" ||KK|| " <<  KK.linfty << endl; 
 nerr += KK.linfty > 1e-9; 
B=0;
B = A*A1; // version 3.13
B -= Id; 
cout <<" ||A*A1-Id|| " <<  B.linfty << endl; 
nerr += B.linfty > 1e-9; 
inv(A1);
A1 -= A; 
cout <<  "|| inv(A1) - A ||" << A1.linfty << endl;
nerr += A1.linfty > 1e-9; 
}
EndMacro 

cout << "Testing real complex dgeev..." << endl;
TEST(real,complex,-100.,dgeev)
cout << nerr << endl;
cout << "Testing complex complex zgeev..." << endl;
TEST(complex,complex,100i,zgeev)
cout << nerr << endl;
cout << "Testing real real dsyev..." << endl;
TEST(real,real,1,dsyev)
cout << nerr << endl;
assert(nerr==0); 


// FFCS - value for regression checks
real regtest=0;

{

int n=5;
real [int,int] A(n,n),  B(n,n),  C(n,n);
SETM(C,i,j,1./(1+i+j)) ; 
SETM(B,i,j,i==j?2.:1./n) ;
 
A = B*C;
cout << A << " " << endl;
A = B + C;
real[int]  b(n),c(n);
real[int,int] AA=A;
AA=A;
real [int,int] At=A';

b = A*c;

At =A';
{

real [int,int] A=[[0,-1,2],[4,11,2],[0,-1,2],[4,11,2]];
real [int,int] B=[[3,-1],[1,2],[6,1]];
real [int,int] E=[[11,0],[35,20],[11,0],[35,20]];
real [int,int]  At=A';
real [int,int] C(A.n,B.m);
C=A*B;
cout << " C = " <<  C << endl;
cout << " E = " <<  E << endl;

C -= E;
assert( C.linfty < 1e-10);
C = A*B;
C -=E; 
assert( C.linfty < 1e-10);
}
{

complex [int,int] A=[[0,-1,2],[4,11,2],[0,-1,2],[4,11,2]];
complex [int,int] B=[[3,-1],[1,2],[6,1]];
complex [int,int] E=[[11,0],[35,20],[11,0],[35,20]];

complex [int,int] C(A.n,B.m);
complex [int,int]  At=A';
cout << "At = " << At << endl;
C=A*B;
cout << " C = " <<  C << endl;
cout << " E = " <<  E << endl;

C -= E;
assert( C.linfty < 1e-10);
C = A*B;
C -=E; 
assert( C.linfty < 1e-10);

// FFCS - value for regression checks
regtest=C.linfty;
}
}
