| 1 | #include "sopnamsp.h"
|
|---|
| 2 | #include "machdefs.h"
|
|---|
| 3 |
|
|---|
| 4 | #include <math.h>
|
|---|
| 5 | #include <iostream>
|
|---|
| 6 |
|
|---|
| 7 | #include "srandgen.h"
|
|---|
| 8 | #include "tarrinit.h"
|
|---|
| 9 | #include "array.h"
|
|---|
| 10 | #include "timing.h"
|
|---|
| 11 | #include "intflapack.h"
|
|---|
| 12 |
|
|---|
| 13 |
|
|---|
| 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 |
|
|---|
| 19 | static double TOLERANCE = 1.e-6;
|
|---|
| 20 | static int nprt = 100;
|
|---|
| 21 | static int prtlev = 1;
|
|---|
| 22 |
|
|---|
| 23 | int main(int narg, char* arg[])
|
|---|
| 24 | {
|
|---|
| 25 |
|
|---|
| 26 | SophyaInit();
|
|---|
| 27 | InitTim(); // Initializing the CPU timer
|
|---|
| 28 |
|
|---|
| 29 |
|
|---|
| 30 |
|
|---|
| 31 | if (narg < 2) {
|
|---|
| 32 | cout << " lpk - LinAlg/LapackServer test \n"
|
|---|
| 33 | << " Usage: lpk select [sizeL,C=5,5] [prtlev=1] \n"
|
|---|
| 34 | << " [nprtmax=100] [WorkSpaceSizeFactor=2] \n"
|
|---|
| 35 | << " select= linsolve svd svds lss inverse / all= \n"
|
|---|
| 36 | << " linsolve: lpk_linsolve_mtx() LapackServer::LinSolve with TMatrix<r_8> \n"
|
|---|
| 37 | << " svd: lpk_svd_mtx() LapackServer::SVD(a,s,u,vt) with TMatrix<r_8> \n"
|
|---|
| 38 | << " svds: lpk_svd_mtx() LapackServer::SVD(a,s); with TMatrix<r_8> \n"
|
|---|
| 39 | << " lss: lpk_leastsquare_mtx() LapackServer::LeastSquareSolve() r_8 \n"
|
|---|
| 40 | << " inverse: LapackServer::ComputeInverse() r_8 \n" << endl;
|
|---|
| 41 | exit(0);
|
|---|
| 42 | }
|
|---|
| 43 | int l,c;
|
|---|
| 44 | l = c = 5;
|
|---|
| 45 | int wsf = 2;
|
|---|
| 46 | string opt = arg[1];
|
|---|
| 47 | if (narg > 2) sscanf(arg[2], "%d,%d", &l, &c);
|
|---|
| 48 | if (narg > 3) prtlev = atoi(arg[3]);
|
|---|
| 49 | if (narg > 4) nprt = atoi(arg[4]);
|
|---|
| 50 | if (narg > 5) wsf = atoi(arg[5]);
|
|---|
| 51 |
|
|---|
| 52 | cout << " lpk - LinAlg/LapackServer_Test sizeL,C=" << l << "," << c
|
|---|
| 53 | << " opt= " << opt << endl;
|
|---|
| 54 | int rc = 0;
|
|---|
| 55 | BaseArray::SetMaxPrint(nprt, prtlev);
|
|---|
| 56 | try {
|
|---|
| 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);
|
|---|
| 64 | else if (opt == "svd") rc = lpk_svd_mtx(l,c,wsf,true);
|
|---|
| 65 | else if (opt == "svds") rc = lpk_svd_mtx(l,c,wsf,false);
|
|---|
| 66 | else if (opt == "lss") rc = lpk_leastsquare_mtx(l,c);
|
|---|
| 67 | else if (opt == "inverse") rc = lpk_inverse_mtx(l);
|
|---|
| 68 | else { cout << " Unknown option " << opt << " ! " << endl; rc = 66; }
|
|---|
| 69 | }
|
|---|
| 70 | catch (PThrowable exc) {
|
|---|
| 71 | cerr << " catched Exception (lpk.cc) " << exc.Msg() << endl;
|
|---|
| 72 | rc = 77;
|
|---|
| 73 | }
|
|---|
| 74 | catch (...) {
|
|---|
| 75 | cerr << " catched unknown (...) exception (lpk.cc) " << endl;
|
|---|
| 76 | rc = 78;
|
|---|
| 77 | }
|
|---|
| 78 |
|
|---|
| 79 | PrtTim(" End of lpk LinAlg/Lapack test ");
|
|---|
| 80 | cout << " --------------- END of Programme -------- (Rc= "
|
|---|
| 81 | << rc << ") --- " << endl;
|
|---|
| 82 | return(rc);
|
|---|
| 83 | }
|
|---|
| 84 |
|
|---|
| 85 | // -----------------------------------------------------------------------
|
|---|
| 86 | /* Nouvelle-Fonction */
|
|---|
| 87 | int lpk_linsolve_mtx(int n)
|
|---|
| 88 | {
|
|---|
| 89 | int i,j;
|
|---|
| 90 | BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
|
|---|
| 91 | cout << " lpk_linsolve_mtx() - Test of LapackServer::LinSolve() " << endl;
|
|---|
| 92 | Matrix a(n,n);
|
|---|
| 93 | for(i=0; i<n; i++)
|
|---|
| 94 | for(j=0; j<n; j++) a(j,i) = GaussianRand(1.,0.);
|
|---|
| 95 |
|
|---|
| 96 | Vector x(n), b;
|
|---|
| 97 | // Matrix x(n,1), b;
|
|---|
| 98 | if (prtlev > 0)
|
|---|
| 99 | cout << " ------------ Vector X = \n " << x << "\n" << endl;
|
|---|
| 100 | for(i=0; i<n; i++) x(i) = GaussianRand(1.5,2.);
|
|---|
| 101 | b = a*x;
|
|---|
| 102 |
|
|---|
| 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 | }
|
|---|
| 109 | cout << "\n Calling LapackLinSolve(a,b) .... " << endl;
|
|---|
| 110 | PrtTim(" Calling LapackLinSolve(a,b) ");
|
|---|
| 111 | LapackLinSolve(a,b);
|
|---|
| 112 | PrtTim(" End LapackLinSolve(a,b) ");
|
|---|
| 113 |
|
|---|
| 114 | if (prtlev > 0)
|
|---|
| 115 | cout << " ------------ Result B(=X ?) = \n " << b << "\n" << endl;
|
|---|
| 116 | Vector diff = b-x;
|
|---|
| 117 | PrtTim(" End of Compute(diff)");
|
|---|
| 118 | if (prtlev > 0)
|
|---|
| 119 | cout << " ------------ Vector diff B-X = \n " << diff << "\n" << endl;
|
|---|
| 120 | double min,max;
|
|---|
| 121 | diff.MinMax(min, max);
|
|---|
| 122 | cout << " Min/Max difference Vector (?=0) , Min= " << min
|
|---|
| 123 | << " Max= " << max << endl;
|
|---|
| 124 | if ((fabs(min) > TOLERANCE) || (fabs(max) > TOLERANCE)) {
|
|---|
| 125 | cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
|
|---|
| 126 | << endl;
|
|---|
| 127 | return(8);
|
|---|
| 128 | }
|
|---|
| 129 | return(0);
|
|---|
| 130 |
|
|---|
| 131 | }
|
|---|
| 132 |
|
|---|
| 133 | // -----------------------------------------------------------------------
|
|---|
| 134 | /* Nouvelle-Fonction */
|
|---|
| 135 | int lpk_svd_mtx(int m, int n, int wsf, bool covu)
|
|---|
| 136 | {
|
|---|
| 137 | BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
|
|---|
| 138 | cout << " lpk_svd_mtx() - Test of LapackServer::SVD " << endl;
|
|---|
| 139 |
|
|---|
| 140 | Matrix a(m , n), aa;
|
|---|
| 141 | a = RandomSequence(RandomSequence::Gaussian, 0., 4.);
|
|---|
| 142 | aa = a;
|
|---|
| 143 | Vector s;
|
|---|
| 144 | Matrix u, vt;
|
|---|
| 145 | cout << " ---- lpk_svd_tmtx() LapackServer::SVD Test Using TMatrix<r_8> ---- " << endl;
|
|---|
| 146 | if (prtlev > 0)
|
|---|
| 147 | cout << " ------------ Matrix A = \n " << a << "\n" << endl;
|
|---|
| 148 |
|
|---|
| 149 | cout << "\n Calling LapackSVD(a,s,u,vt) .... " << endl;
|
|---|
| 150 | PrtTim(" Calling LapackSVD() ");
|
|---|
| 151 | LapackServer<r_8> lpks;
|
|---|
| 152 | lpks.SetWorkSpaceSizeFactor(wsf);
|
|---|
| 153 | if (covu) lpks.SVD(aa, s, u, vt);
|
|---|
| 154 | else lpks.SVD(aa, s);
|
|---|
| 155 | PrtTim(" End LapackSVD() ");
|
|---|
| 156 |
|
|---|
| 157 | if (prtlev > 0)
|
|---|
| 158 | cout << " ------------ Result S = \n " << s << "\n" << endl;
|
|---|
| 159 | if (!covu) return(0);
|
|---|
| 160 | if (prtlev > 0) {
|
|---|
| 161 | cout << " ------------ Result U = \n " << u << "\n" << endl;
|
|---|
| 162 | cout << " ------------ Result VT = \n " << vt << "\n" << endl;
|
|---|
| 163 | }
|
|---|
| 164 | Matrix sm(m,n);
|
|---|
| 165 | int minmn = (m<n) ? m : n ;
|
|---|
| 166 | for(int k=0; k< minmn ; k++) sm(k,k) = s(k);
|
|---|
| 167 | Matrix diff = u*(sm*vt) - a;
|
|---|
| 168 | PrtTim(" End of Compute(diff)");
|
|---|
| 169 | if (prtlev > 0)
|
|---|
| 170 | cout << " ------------ Matrix diff U*S*Vt - A = \n " << diff << "\n" << endl;
|
|---|
| 171 | double min,max;
|
|---|
| 172 | diff.MinMax(min, max);
|
|---|
| 173 | cout << " Min/Max difference Matrix (?=0) , Min= " << min
|
|---|
| 174 | << " Max= " << max << endl;
|
|---|
| 175 | if ((fabs(min) > TOLERANCE) || (fabs(max) > TOLERANCE)) {
|
|---|
| 176 | cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
|
|---|
| 177 | << endl;
|
|---|
| 178 | return(16);
|
|---|
| 179 | }
|
|---|
| 180 | return(0);
|
|---|
| 181 |
|
|---|
| 182 | }
|
|---|
| 183 |
|
|---|
| 184 | // -----------------------------------------------------------------------
|
|---|
| 185 | /* Nouvelle-Fonction */
|
|---|
| 186 | int lpk_leastsquare_mtx(int m, int n)
|
|---|
| 187 | {
|
|---|
| 188 | BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
|
|---|
| 189 | cout << " lpk_leastsquare_mtx() - Test of LapackLeastSquareSolve " << endl;
|
|---|
| 190 |
|
|---|
| 191 | Matrix a(m , n), aa;
|
|---|
| 192 | a = RandomSequence(RandomSequence::Gaussian, 0., 4.);
|
|---|
| 193 | aa = a;
|
|---|
| 194 | Vector x(n),noise(m),b,bx;
|
|---|
| 195 | x = RandomSequence(RandomSequence::Gaussian, 2., 1.);
|
|---|
| 196 | double signoise = 0.1;
|
|---|
| 197 | noise = RandomSequence(RandomSequence::Gaussian, 0., 0.1);
|
|---|
| 198 | b = a*x;
|
|---|
| 199 | bx = b+noise;
|
|---|
| 200 |
|
|---|
| 201 | cout << " ---- lpk_leastsquare_tmtx() LapackLeastSquareSolve<r_8> ---- " << endl;
|
|---|
| 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 | }
|
|---|
| 207 | cout << "\n Calling LapackLeastSquareSolve(a,bx) .... " << endl;
|
|---|
| 208 | PrtTim(" Calling LapackLeastSquareSolve() ");
|
|---|
| 209 | LapackLeastSquareSolve(aa, bx);
|
|---|
| 210 | PrtTim(" End LapackLeastSquareSolve() ");
|
|---|
| 211 | bx.Share(bx(Range(0,0,n)));
|
|---|
| 212 | if (prtlev > 0) {
|
|---|
| 213 | cout << " ------------ Result X = \n " << bx << "\n" << endl;
|
|---|
| 214 | cout << " ------------ X-X_Real = \n " << bx-x << "\n" << endl;
|
|---|
| 215 | }
|
|---|
| 216 | Vector diff = b-a*bx;
|
|---|
| 217 | PrtTim(" End of Compute(diff)");
|
|---|
| 218 | if (prtlev > 0)
|
|---|
| 219 | cout << " ------------ Matrix diff b-a*x = \n " << diff << "\n" << endl;
|
|---|
| 220 | double mean,sigma;
|
|---|
| 221 | MeanSigma(diff, mean, sigma);
|
|---|
| 222 | cout << " MeanSigma(diff) , Mean= " << mean << " Sigma=" << sigma
|
|---|
| 223 | << " SigNoise= " << signoise << endl;
|
|---|
| 224 | if ((fabs(mean) > signoise) || (fabs(sigma-signoise) > signoise)) {
|
|---|
| 225 | cout << " !!! Difference exceeding tolerance (=" << TOLERANCE << ") !!!"
|
|---|
| 226 | << endl;
|
|---|
| 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 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
