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