1load "lapack"
2// test of lapack load file ...
3// .......................
4// load "fflapack" obsolete (F. Hecht version 3.8)
5// use  load "lapack"
6int nerr=0; // nomber of err ..
7//  to set a full matrix ..
8macro SETM(A,i,j,fij)
9{
10    for(int i=0;i<A.n;++i)
11      for(int j=0;j<A.m;++j)
12        A(i,j)= (fij) ;
13}//EOM
14
15// The of Eigen Problem ..
16NewMacro TEST(Type,Typevp,vtype,ComputeEV)
17{
18cout << "\n\n *****************   VP \n\n\n";
19int n=5;
20Type[int,int] A(n,n),A1(n,n),B(n,n),Id(n,n);
21SETM(A,i,j,(i==j) ? n+1 : 1);
22SETM(Id,i,j,real(i==j));
23A(0,n-1)=vtype;
24
25cout << A << endl;
26A1=A^-1;
27cout << A1 << endl;
28
29Typevp[int] vp(n);
30Typevp[int,int] VP(n,n),KK(n,n);
31
32int nn= ComputeEV(A,vp,VP);
33 cout << " vp = " << vp << endl;
34 cout << " VP = " << VP << endl;
35
36 // verification ...
37 KK =0.;
38 for(int i=0;i<n;++i)
39   for(int j=0;j<n;++j)
40     for(int k=0;k<n;++k)
41       //KK(i,j) += (A(i,k) - vp[j]* real(i==k) ) *VP(k,j);
42       KK(i,j) += (A(i,k) - ((i==k)?vp[j]:0.0) ) *VP(k,j);
43 cout <<" ||KK|| " <<  KK.linfty << endl;
44 nerr += KK.linfty > 1e-9;
45B=0;
46B = A*A1; // version 3.13
47B -= Id;
48cout <<" ||A*A1-Id|| " <<  B.linfty << endl;
49nerr += B.linfty > 1e-9;
50inv(A1);
51A1 -= A;
52cout <<  "|| inv(A1) - A ||" << A1.linfty << endl;
53nerr += A1.linfty > 1e-9;
54}
55EndMacro
56
57cout << "Testing real complex dgeev..." << endl;
58TEST(real,complex,-100.,dgeev)
59cout << nerr << endl;
60cout << "Testing complex complex zgeev..." << endl;
61TEST(complex,complex,100i,zgeev)
62cout << nerr << endl;
63cout << "Testing real real dsyev..." << endl;
64TEST(real,real,1,dsyev)
65cout << nerr << endl;
66assert(nerr==0);
67
68
69// FFCS - value for regression checks
70real regtest=0;
71
72{
73
74int n=5;
75real [int,int] A(n,n),  B(n,n),  C(n,n);
76SETM(C,i,j,1./(1+i+j)) ;
77SETM(B,i,j,i==j?2.:1./n) ;
78
79A = B*C;
80cout << A << " " << endl;
81A = B + C;
82real[int]  b(n),c(n);
83real[int,int] AA=A;
84AA=A;
85real [int,int] At=A';
86
87b = A*c;
88
89At =A';
90{
91
92real [int,int] A=[[0,-1,2],[4,11,2],[0,-1,2],[4,11,2]];
93real [int,int] B=[[3,-1],[1,2],[6,1]];
94real [int,int] E=[[11,0],[35,20],[11,0],[35,20]];
95real [int,int]  At=A';
96real [int,int] C(A.n,B.m);
97C=A*B;
98cout << " C = " <<  C << endl;
99cout << " E = " <<  E << endl;
100
101C -= E;
102assert( C.linfty < 1e-10);
103C = A*B;
104C -=E;
105assert( C.linfty < 1e-10);
106}
107{
108
109complex [int,int] A=[[0,-1,2],[4,11,2],[0,-1,2],[4,11,2]];
110complex [int,int] B=[[3,-1],[1,2],[6,1]];
111complex [int,int] E=[[11,0],[35,20],[11,0],[35,20]];
112
113complex [int,int] C(A.n,B.m);
114complex [int,int]  At=A';
115cout << "At = " << At << endl;
116C=A*B;
117cout << " C = " <<  C << endl;
118cout << " E = " <<  E << endl;
119
120C -= E;
121assert( C.linfty < 1e-10);
122C = A*B;
123C -=E;
124assert( C.linfty < 1e-10);
125
126// FFCS - value for regression checks
127regtest=C.linfty;
128}
129}
130