Changeset 2647 in Sophya
- Timestamp:
- Feb 7, 2005, 5:43:49 PM (21 years ago)
- Location:
- trunk/SophyaProg/Tests
- Files:
-
- 1 added
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaProg/Tests/lpk.cc
r2615 r2647 12 12 13 13 14 int lpk_ rzt_arr(int n);15 int lpk_ linsolve_mtx(int n);16 int lpk_ svd_mtx(int l, int c, int wsf, bool covu=true);17 int lpk_ leastsquare_mtx(int l, int c);14 int lpk_linsolve_mtx(int n); // ErrorCode = 8 15 int lpk_svd_mtx(int l, int c, int wsf, bool covu=true); // ErrorCode = 16 16 int lpk_leastsquare_mtx(int l, int c); // ErrorCode = 32 17 int lpk_inverse_mtx(int l); // ErrorCode = 64 18 18 19 19 static double TOLERANCE = 1.e-6; 20 static int nprt = 100; 21 static int prtlev = 1; 20 22 21 23 int main(int narg, char* arg[]) … … 29 31 if (narg < 2) { 30 32 cout << " lpk - LinAlg/LapackServer test \n" 31 << " Usage: lpk linsolve/svd/rzt [sizeL,C=5,5] [prtlev=0] \n"33 << " Usage: lpk select [sizeL,C=5,5] [prtlev=1] \n" 32 34 << " [nprtmax=100] [WorkSpaceSizeFactor=2] \n" 35 << " select= linsolve svd svds lss inverse / all= \n" 33 36 << " linsolve: lpk_linsolve_mtx() LapackServer::LinSolve with TMatrix<r_8> \n" 34 37 << " svd: lpk_svd_mtx() LapackServer::SVD(a,s,u,vt) with TMatrix<r_8> \n" 35 38 << " svds: lpk_svd_mtx() LapackServer::SVD(a,s); with TMatrix<r_8> \n" 36 39 << " lss: lpk_leastsquare_mtx() LapackServer::LeastSquareSolve() r_8 \n" 37 << " rzt: lpk_rzt_tarr() rztest_lapack with TArray<r_4>\n" << endl;40 << " inverse: LapackServer::ComputeInverse() r_8 \n" << endl; 38 41 exit(0); 39 42 } … … 43 46 string opt = arg[1]; 44 47 if (narg > 2) sscanf(arg[2], "%d,%d", &l, &c); 45 int nprt = 100;46 int prtlev = 1;47 48 if (narg > 3) prtlev = atoi(arg[3]); 48 49 if (narg > 4) nprt = atoi(arg[4]); 49 50 if (narg > 5) wsf = atoi(arg[5]); 50 51 52 cout << " lpk - LinAlg/LapackServer_Test sizeL,C=" << l << "," << c 53 << " opt= " << opt << endl; 51 54 int rc = 0; 52 55 BaseArray::SetMaxPrint(nprt, prtlev); 53 56 try { 54 if (opt == "linsolve") rc = lpk_linsolve_mtx(l); 57 if (opt == "all") { 58 rc += lpk_linsolve_mtx(l); 59 rc += lpk_svd_mtx(l,c,wsf,true); 60 rc += lpk_leastsquare_mtx(l,c); 61 rc += lpk_inverse_mtx(l); 62 } 63 else if (opt == "linsolve") rc = lpk_linsolve_mtx(l); 55 64 else if (opt == "svd") rc = lpk_svd_mtx(l,c,wsf,true); 56 65 else if (opt == "svds") rc = lpk_svd_mtx(l,c,wsf,false); 57 66 else if (opt == "lss") rc = lpk_leastsquare_mtx(l,c); 58 else if (opt == " rzt") rc = lpk_rzt_arr(l);67 else if (opt == "inverse") rc = lpk_inverse_mtx(l); 59 68 else { cout << " Unknown option " << opt << " ! " << endl; rc = 66; } 60 69 } … … 87 96 Vector x(n), b; 88 97 // Matrix x(n,1), b; 89 cout << " ------------ Vector X = \n " << x << "\n" << endl; 98 if (prtlev > 0) 99 cout << " ------------ Vector X = \n " << x << "\n" << endl; 90 100 for(i=0; i<n; i++) x(i) = GauRnd(2., 1.5); 91 101 b = a*x; 92 102 93 cout << " ---- lpk_tmtx() LapackServer::LinSolve Test Using TMatrix<r_8> ----- " << endl; 94 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 95 cout << " ------------ Matrix X = \n " << x << "\n" << endl; 96 cout << " ------------ Matrix B = \n " << b << "\n" << endl; 97 103 if (prtlev > 0) { 104 cout << " ---- lpk_tmtx() LapackServer::LinSolve Test Using TMatrix<r_8> ----- " << endl; 105 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 106 cout << " ------------ Matrix X = \n " << x << "\n" << endl; 107 cout << " ------------ Matrix B = \n " << b << "\n" << endl; 108 } 98 109 cout << "\n Calling LapackLinSolve(a,b) .... " << endl; 99 110 PrtTim(" Calling LapackLinSolve(a,b) "); … … 101 112 PrtTim(" End LapackLinSolve(a,b) "); 102 113 103 cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl; 114 if (prtlev > 0) 115 cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl; 104 116 Vector diff = b-x; 105 117 PrtTim(" End of Compute(diff)"); 106 cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl; 118 if (prtlev > 0) 119 cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl; 107 120 double min,max; 108 121 diff.MinMax(min, max); … … 112 125 cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 113 126 << endl; 114 return( 99);127 return(8); 115 128 } 116 129 return(0); … … 131 144 Matrix u, vt; 132 145 cout << " ---- lpk_svd_tmtx() LapackServer::SVD Test Using TMatrix<r_8> ---- " << endl; 133 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 146 if (prtlev > 0) 147 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 134 148 135 149 cout << "\n Calling LapackSVD(a,s,u,vt) .... " << endl; … … 141 155 PrtTim(" End LapackSVD() "); 142 156 143 cout << " ------------ Result S = \n " << s << "\n" << endl; 157 if (prtlev > 0) 158 cout << " ------------ Result S = \n " << s << "\n" << endl; 144 159 if (!covu) return(0); 145 cout << " ------------ Result U = \n " << u << "\n" << endl; 146 cout << " ------------ Result VT = \n " << vt << "\n" << endl; 160 if (prtlev > 0) { 161 cout << " ------------ Result U = \n " << u << "\n" << endl; 162 cout << " ------------ Result VT = \n " << vt << "\n" << endl; 163 } 147 164 Matrix sm(m,n); 148 165 int minmn = (m<n) ? m : n ; … … 150 167 Matrix diff = u*(sm*vt) - a; 151 168 PrtTim(" End of Compute(diff)"); 152 cout << " ------------ Matrix diff U*S*Vt - A = \n " << diff << "\n" << endl; 169 if (prtlev > 0) 170 cout << " ------------ Matrix diff U*S*Vt - A = \n " << diff << "\n" << endl; 153 171 double min,max; 154 172 diff.MinMax(min, max); … … 158 176 cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 159 177 << endl; 160 return( 99);178 return(16); 161 179 } 162 180 return(0); … … 182 200 183 201 cout << " ---- lpk_leastsquare_tmtx() LapackLeastSquareSolve<r_8> ---- " << endl; 184 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 185 cout << " ------------ Vector B = \n " << bx << "\n" << endl; 186 cout << " ------------ Vector X = \n " << x << "\n" << endl; 187 202 if (prtlev > 0) { 203 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 204 cout << " ------------ Vector B = \n " << bx << "\n" << endl; 205 cout << " ------------ Vector X = \n " << x << "\n" << endl; 206 } 188 207 cout << "\n Calling LapackLeastSquareSolve(a,bx) .... " << endl; 189 208 PrtTim(" Calling LapackLeastSquareSolve() "); … … 191 210 PrtTim(" End LapackLeastSquareSolve() "); 192 211 bx.Share(bx(Range(0,0,n))); 193 cout << " ------------ Result X = \n " << bx << "\n" << endl; 194 cout << " ------------ X-X_Real = \n " << bx-x << "\n" << endl; 212 if (prtlev > 0) { 213 cout << " ------------ Result X = \n " << bx << "\n" << endl; 214 cout << " ------------ X-X_Real = \n " << bx-x << "\n" << endl; 215 } 195 216 Vector diff = b-a*bx; 196 217 PrtTim(" End of Compute(diff)"); 197 cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl; 218 if (prtlev > 0) 219 cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl; 198 220 double mean,sigma; 199 221 MeanSigma(diff, mean, sigma); … … 203 225 cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 204 226 << endl; 205 return(99); 206 } 207 return(0); 208 209 } 210 211 212 void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb); 213 214 int lpk_rzt_arr(int n) 215 { 216 int i,j; 217 TArray<r_4> a(n,n); 218 for(i=0; i<n; i++) 219 for(j=0; j<n; j++) a(i,j,0) = GauRnd(0., 1.); 220 221 TArray<r_4> x(n,1), b(n,1); 222 r_4 sum ; 223 for(i=0; i<n; i++) x(i,0,0) = GauRnd(2., 1.5); 224 for(i=0; i<n; i++) { 225 sum = 0.; 226 for(j=0; j<n; j++) sum += a(i,j,0)*x(j,0,0); 227 b(i,0,0) = sum; 228 } 229 230 // cout << ":::::::: rztest_lapack - Size=" << n << " ::::::::: " << endl; 231 cout << " ------- lpk_tarr() LapaackTest Using TArray<r_4> -------- " << endl; 232 cout << " ------------ Array A = \n " << a << "\n" << endl; 233 cout << " ------------ Array X = \n " << x << "\n" << endl; 234 cout << " ------------ Array B = \n " << b << "\n" << endl; 235 236 cout << "\n Calling rztest_lapack ... " << endl; 237 238 rztest_lapack(a, b); 239 240 cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl; 241 return(0); 242 } 227 return 32; 228 } 229 return 0; 230 231 } 232 // ----------------------------------------------------------------------- 233 /* Nouvelle-Fonction */ 234 int lpk_inverse_mtx(int m) 235 { 236 BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping); 237 cout << " lpk_inverse_mtx() - Test of LapackComputeInverse<r_8> " << endl; 238 239 Matrix a(m , m), aa, ainv, idmx(m,m), mxprod, diff; 240 a = RandomSequence(RandomSequence::Gaussian, 0., 4.); 241 aa = a; // We make a copy of a, as it is modified by LapackServer 242 PrtTim(" End of a_inint"); 243 244 cout << "\n Calling LapackInverse() .... " << endl; 245 ainv = LapackInverse(aa); 246 PrtTim(" End of LapackInverse(a)"); 247 mxprod = a*ainv; 248 PrtTim(" End of mxprod = a*ainv"); 249 idmx = IdentityMatrix(); 250 PrtTim(" End of idmx = IdentityMatrix()"); 251 diff = mxprod-idmx; 252 PrtTim(" End of Compute(diff)"); 253 254 if (prtlev > 0) { 255 cout << " ------------ Matrix A = \n " << a << "\n" << endl; 256 cout << " ------------ Matrix Inverse(A) = \n " << ainv << "\n" << endl; 257 cout << " ------------ Matrix mxprod = A*Inverse(A) = \n " << mxprod << "\n" << endl; 258 } 259 double min,max; 260 diff.MinMax(min, max); 261 PrtTim(" End of diffCheck"); 262 cout << " Min/Max difference Matrix (?=0) , Min= " << min 263 << " Max= " << max << endl; 264 if ((fabs(min) > TOLERANCE) || (fabs(max) > TOLERANCE)) { 265 cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!" 266 << endl; 267 return 64; 268 } 269 return 0; 270 271 } 272
Note:
See TracChangeset
for help on using the changeset viewer.