CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // $Id: testMatrix.cc,v 1.3 2003/08/13 20:00:12 garren Exp $ 00003 // 00004 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00005 // 00006 // This is a small program for testing the classes from the HEP Matrix module. 00007 // 00008 00009 #include "CLHEP/Matrix/defs.h" 00010 #include "CLHEP/Matrix/Matrix.h" 00011 #include "CLHEP/Matrix/SymMatrix.h" 00012 #include "CLHEP/Matrix/DiagMatrix.h" 00013 #include "CLHEP/Matrix/Vector.h" 00014 #include "CLHEP/Random/Random.h" 00015 #include "CLHEP/Random/JamesRandom.h" 00016 #include "CLHEP/Random/RandFlat.h" 00017 #include <iostream> 00018 #include <iomanip> 00019 00020 using std::cout; 00021 using std::endl; 00022 00023 using namespace CLHEP; 00024 00025 #define PRINTOUT 00026 00027 int matrix_test1(const HepGenMatrix&m) { 00028 // 00029 // test of virtual finctions. 00030 // 00031 static int dum = 0; 00032 dum += m.num_col(); 00033 return 0; 00034 } 00035 00036 // 00037 // test function 00038 // 00039 double neg(double f, int, int) { 00040 return -f; 00041 } 00042 00043 double absf(double f, int, int) { 00044 return fabs(f); 00045 } 00046 00047 double negv(double f, int) { 00048 return -f; 00049 } 00050 00051 void matrix_test2(const HepSymMatrix&m, HepSymMatrix &n) { 00052 n = m.sub(2,5); 00053 } 00054 00055 00056 void matrix_test() { 00057 // 00058 // complete test of all functions in Matrix.cc 00059 // 00060 cout << "Starting Matrix tests" << endl; 00061 { 00062 // Test of HepMatrix() 00063 HepMatrix a; 00064 } 00065 00066 { 00067 // Test of HepMatrix(p,q) 00068 HepMatrix a(3,5); 00069 HepMatrix b(4,5,0); 00070 cout << "4x5 matrix initialized to zero " << b; 00071 HepMatrix c(3,3,1); 00072 cout << "3x3 matrix initialized to identity " << c; 00073 } 00074 00075 { 00076 // Test of HepMatrix(p,q,Random) 00077 HepRandom r; 00078 r.setTheSeed(31); 00079 HepMatrix a(3,3,r); 00080 cout << "3x3 matrix initialized with Random " << a; 00081 // Test of HepMatrix(const HepMatrix&); 00082 HepMatrix b(a); 00083 cout << "matrix initialized to the previous matrix " << b; 00084 } 00085 00086 { 00087 // Test of sub matrix 00088 HepRandom r; 00089 r.setTheSeed(31); 00090 HepMatrix a(8,5,r); 00091 HepMatrix b = a.sub(2,6,3,5); 00092 cout << "8x5 matrix" << a; 00093 cout << "sub matrix (2-6 x 3-5 of the previous matrix." << b; 00094 HepMatrix c(2,3,0); 00095 a.sub(4,3,c); 00096 cout << "embedding sub matrix at 4,3" << a; 00097 // Test of dsum 00098 HepMatrix d(3,2,r); 00099 cout << "dsum" << dsum(a,d); 00100 } 00101 00102 { 00103 // m += m1; 00104 HepRandom r; 00105 r.setTheSeed(31); 00106 HepMatrix a(5,5,r); 00107 HepMatrix b(5,5,r); 00108 cout << "a" << a; 00109 cout << "b" << b; 00110 cout << "a += b" << (a+=b); 00111 HepMatrix c; 00112 c = a + b; 00113 cout << "a + b" << c; 00114 } 00115 00116 { 00117 // test of T(); 00118 HepRandom r; 00119 r.setTheSeed(31); 00120 HepMatrix a(5,3,r); 00121 HepMatrix b = a.T(); 00122 cout << "a" << a; 00123 cout << "a.T" << b; 00124 } 00125 00126 cout << "End of Matrix tests" << endl; 00127 00128 } 00129 00130 void symmatrix_test() { 00131 { 00132 // Test of HepSymMatrix() 00133 HepSymMatrix a; 00134 } 00135 00136 { 00137 // Test of HepSymMatrix(p) 00138 HepSymMatrix a(3); 00139 HepSymMatrix b(4,0); 00140 cout << "4x4 Symmetric matrix initialuzed to zero " << b; 00141 HepSymMatrix c(3,1); 00142 cout << "3x3 Symmetric matrix initialized to identity "<< c; 00143 } 00144 00145 { 00146 // Test of HepSymMatrix(p,Random) 00147 HepRandom r; 00148 r.setTheSeed(31); 00149 HepSymMatrix a(3,r); 00150 cout << "3x3 symmetric matrix initialized with Random " << a << endl; 00151 // Test of HepSymMatrix(const HepSymMatrix&); 00152 HepSymMatrix b(a); 00153 cout << "symmetric matrix initialized to the previous matrix " << b; 00154 HepMatrix c(b); 00155 cout << "matrix initialized to the previous symmetric matrix " 00156 << c; 00157 // 00158 // Access to elements 00159 // 00160 double f = 3.8; 00161 double g = 22.5; 00162 cout << c(1,1) << " " << c[0][0] << endl; 00163 c(1,2) = f; 00164 c[2][1] = g; 00165 c(2,3) = f; 00166 c[1][2] = g; 00167 cout << c << endl; 00168 HepMatrix &d = c; 00169 cout << d(1,1) << " " << d[0][0] << endl; 00170 00171 } 00172 00173 { 00174 // Test of sub symmatrix 00175 HepRandom r; 00176 r.setTheSeed(31); 00177 HepSymMatrix a(5,r); 00178 HepSymMatrix b = a.sub(2,5); 00179 cout << "5x5 sym matrix" << a; 00180 cout << "sub sym matrix (2-5 x 2-5 of the previous sub matrix." << b; 00181 HepSymMatrix c(3,0); 00182 a.sub(2,c); 00183 cout << "embedding sub matrix at 2" << a; 00184 } 00185 { 00186 // m = m1 + s; 00187 HepRandom r; 00188 r.setTheSeed(31); 00189 HepMatrix a(5,5,r); 00190 HepSymMatrix b(5,r); 00191 cout << "a" << a; 00192 cout << "b(sym)" << b; 00193 cout << "a += b" << (a+=b); 00194 HepMatrix c = a + b; 00195 cout << "a + b" << c; 00196 } 00197 { 00198 // test of similarity(Matrix) 00199 HepRandom r; 00200 r.setTheSeed(31); 00201 HepMatrix a(5,3,r); 00202 HepSymMatrix b(3,r); 00203 HepSymMatrix c = b.similarity(a); 00204 cout << "a" << a; 00205 cout << "b" << b; 00206 cout << "c" << c; 00207 } 00208 { 00209 // test of similarityT(Matrix) 00210 HepRandom r; 00211 r.setTheSeed(31); 00212 HepMatrix a(3,5,r); 00213 HepSymMatrix b(3,r); 00214 HepSymMatrix c = b.similarityT(a); 00215 cout << "a" << a; 00216 cout << "b" << b; 00217 cout << "c" << c; 00218 } 00219 { 00220 // test of similarity(SymMatrix) 00221 HepRandom r; 00222 r.setTheSeed(31); 00223 HepSymMatrix a(3,r); 00224 HepSymMatrix b(3,r); 00225 HepSymMatrix c = b.similarity(a); 00226 cout << "a" << a; 00227 cout << "b" << b; 00228 cout << "c" << c; 00229 } 00230 { 00231 // test of similarity(Vector) 00232 HepRandom r; 00233 r.setTheSeed(31); 00234 HepVector a(3,r); 00235 HepSymMatrix b(3,r); 00236 double c = b.similarity(a); 00237 HepSymMatrix cc = b.similarity(HepMatrix(a.T())); 00238 cout << "a" << a; 00239 cout << "b" << b; 00240 cout << "c\t" << c << endl; 00241 cout << "c" << cc; 00242 } 00243 cout << "End of SymMatrix tests" << endl; 00244 00245 } 00246 00247 void diagmatrix_test() { 00248 { 00249 // Test of HepDiagMatrix() 00250 HepDiagMatrix a; 00251 } 00252 00253 { 00254 // Test of HepDiagMatrix(p) 00255 HepDiagMatrix a(3); 00256 HepDiagMatrix b(4,0); 00257 cout << "4x4 diagonal matrix initialized to zero " << b; 00258 HepDiagMatrix c(3,1); 00259 cout << "3x3 diagonal matrix initialized to identity " << c; 00260 } 00261 00262 { 00263 // Test of HepDiagMatrix(p,Random) 00264 HepRandom r; 00265 r.setTheSeed(31); 00266 HepDiagMatrix a(3,r); 00267 cout << "3x3 diagonal matrix initialized to Random " << a; 00268 // Test of HepDiagMatrix(const HepDiagMatrix&); 00269 HepDiagMatrix b(a); 00270 cout << "diagonal matrix initialized to the previous matrix " << b; 00271 HepMatrix c(b); 00272 cout << "matrix initialized to the previous diagonal matrix " 00273 << c; 00274 HepSymMatrix d(b); 00275 cout << "Symmetric matrix initialized to the previous diagonal matrix " 00276 << d; 00277 } 00278 00279 { 00280 // Test of sub diagmatrix 00281 HepRandom r; 00282 r.setTheSeed(31); 00283 HepDiagMatrix a(8,r); 00284 HepDiagMatrix b = a.sub(2,5); 00285 cout << "8x8 diag matrix" << a; 00286 cout << "sub diag matrix (2-5 x 2-5 of the previous diag matrix." << b; 00287 HepDiagMatrix c(3,0); 00288 a.sub(2,c); 00289 cout << "embedding sub matrix at 2" << a; 00290 } 00291 { 00292 // m = m1 + s; 00293 HepRandom r; 00294 r.setTheSeed(31); 00295 HepMatrix a(5,5,r); 00296 HepDiagMatrix b(5,r); 00297 cout << "a" << a; 00298 cout << "b(diag)" << b; 00299 cout << "a += b" << (a+=b); 00300 HepMatrix c = a + b; 00301 cout << "a + b" << c; 00302 } 00303 00304 cout << "End of DiagMatrix tests" << endl; 00305 } 00306 00307 void vector_test() { 00308 { 00309 // Test of HepVector() 00310 HepVector a; 00311 } 00312 00313 { 00314 // Test of HepVector(p) 00315 HepVector a(3); 00316 HepVector b(4,0); 00317 cout << "4 vector initialized to zero "<< b; 00318 HepVector c(3,1); 00319 cout << "3 vector initialized to identity " << c; 00320 } 00321 00322 { 00323 // Test of HepVector(p,Random) 00324 HepRandom r; 00325 r.setTheSeed(31); 00326 HepVector a(3,r); 00327 cout << "3 vector initialized to Random " << a; 00328 // Test of HepVector(const HepVector&); 00329 HepVector b(a); 00330 cout << "Vector initialized to the previous vector " << b; 00331 HepMatrix c(b); 00332 cout << "matrix initialized to the previous vector " 00333 << c; 00334 HepVector d(c); 00335 cout << "Vector initialized to the previous matrix " 00336 << d; 00337 } 00338 00339 { 00340 // Test of sub diagmatrix 00341 HepRandom r; 00342 r.setTheSeed(31); 00343 HepVector a(8,r); 00344 HepVector b = a.sub(2,5); 00345 cout << "8 vector" << a; 00346 cout << "sub vector (2-5 of the previous vector." << b; 00347 HepVector c(3,0); 00348 a.sub(2,c); 00349 cout << "embedding sub vector at 2 " << a; 00350 } 00351 { 00352 // m = m1 + s; 00353 HepRandom r; 00354 r.setTheSeed(31); 00355 HepMatrix a(5,1,r); 00356 HepVector b(5,r); 00357 cout << "a" << a; 00358 cout << "b(vector)" << b; 00359 cout << "a += b" << (a+=b); 00360 HepMatrix c = a + b; 00361 cout << "a + b" << c; 00362 } 00363 00364 cout << "End of Vector tests" << endl; 00365 } 00366 00367 int main() { 00368 // 00369 // Test of HepMatrix 00370 // 00371 cout << std::setiosflags(std::ios::fixed) << std::setw(10); 00372 matrix_test(); 00373 symmatrix_test(); 00374 diagmatrix_test(); 00375 vector_test(); 00376 00377 // 00378 // Test of default constructor 00379 // 00380 HepMatrix m; 00381 // 00382 // Test of constructors. 00383 // 00384 HepMatrix d(3,3,1); 00385 HepMatrix e(3,3,0); 00386 HepMatrix f(5,3); 00387 HepMatrix g(f); 00388 matrix_test1(g); 00389 // 00390 // Test of constructor with a Random object. 00391 // 00392 HepRandom r; 00393 r.setTheSeed(31); 00394 HepMatrix a(3,3,r); 00395 #if defined(PRINTOUT) 00396 cout << a << endl; 00397 #endif 00398 HepMatrix b(3,5,r); 00399 matrix_test1(b); 00400 #if defined(PRINTOUT) 00401 cout << b << endl; 00402 #endif 00403 00404 HepSymMatrix dds(3,r); 00405 HepMatrix ddm(dds); 00406 HepDiagMatrix ddd(3,r); 00407 HepMatrix ddmd(ddd); 00408 HepVector ddv(3,r); 00409 HepMatrix ddmv(ddv); 00410 #if defined(PRINTOUT) 00411 cout << "nraw=" << b.num_row() << " ncol=" << b.num_col() << endl; 00412 #endif 00413 #if defined(PRINTOUT) 00414 cout << "b(1,1)=" << b(1,1) << " b(2,1)=" << b(2,1) << endl; 00415 #endif 00416 b(2,3) = 1.0; 00417 00418 b /= 3.0; 00419 b *= 6.0; 00420 00421 HepSymMatrix sm(3,r); 00422 HepDiagMatrix dm(3,r); 00423 HepVector vvm(3,r); 00424 00425 d += e; 00426 d += sm; 00427 d += dm; 00428 ddmv += vvm; 00429 00430 a -= e; 00431 a -= sm; 00432 a -= dm; 00433 ddmv -= vvm; 00434 00435 d = sm; 00436 d = dm; 00437 d = a; 00438 ddmv = vvm; 00439 00440 e = d.apply(neg); 00441 00442 a = d.T(); 00443 00444 e = b.sub(1,2,3,2); 00445 00446 b.sub(1,3,e); 00447 00448 swap(a,d); 00449 00450 m = d * a; 00451 m = d * sm; 00452 m = d * dm; 00453 m = sm * d; 00454 #if defined(PRINTOUT) 00455 cout << "Dm=" << dm << " d = " << d << endl; 00456 #endif 00457 m = dm * d; 00458 m = sm * sm; 00459 m = dm * dm; 00460 00461 // 00462 // SymMatrix 00463 // 00464 HepSymMatrix s; 00465 HepSymMatrix ss(6); 00466 HepSymMatrix si(6,1); 00467 HepSymMatrix sz(6,0); 00468 HepSymMatrix sr(6,r); 00469 HepSymMatrix sc(sr); 00470 HepSymMatrix scd(dm); 00471 00472 matrix_test1(si); 00473 #if defined(PRINTOUT) 00474 cout << "nraw=" << sr.num_row() << " ncol=" << sr.num_col() << endl; 00475 #endif 00476 #if defined(PRINTOUT) 00477 cout << "sr(1,1)=" << sr(1,1) << " sr(2,1)=" << sr(2,1) << endl; 00478 #endif 00479 00480 sr(4,3) = r(); 00481 sr.fast(3,1) = r(); 00482 00483 s.assign(d); 00484 #if defined(PRINTOUT) 00485 cout << "d=" << d << "s=" << s << endl; 00486 #endif 00487 ss.assign(sc); 00488 00489 ss *= 0.5; 00490 ss /= 0.2; 00491 00492 HepDiagMatrix ddds(ss.num_row(),r); 00493 ss += si; 00494 ss += ddds; 00495 ss -= sr; 00496 ss -= ddds; 00497 ss = ddds; 00498 00499 s = sr; 00500 s = ss.T(); 00501 s.apply(neg); 00502 HepMatrix mm(6,4,r); 00503 ss = s.similarityT(mm); 00504 00505 HepMatrix mt(8,6,r); 00506 ss = s.similarity(mt); 00507 00508 s.assign(d); 00509 00510 s = sr.sub(2,5); 00511 00512 ss.sub(2,s); 00513 s = ss.sub(2,7); 00514 00515 m = s * sr; 00516 m = s * m; 00517 m = m * s; 00518 00519 HepVector v(6,r); 00520 s = vT_times_v(v); 00521 00522 // 00523 // Test of HepDiagMatrix 00524 // 00525 00526 HepDiagMatrix dt; 00527 HepDiagMatrix dn(6); 00528 HepDiagMatrix di(6,1); 00529 HepDiagMatrix dz(6,0); 00530 HepDiagMatrix dr(6,r); 00531 HepDiagMatrix dc(dr); 00532 00533 matrix_test1(dr); 00534 #if defined(PRINTOUT) 00535 cout << "Nr=" << di.num_row() << " Nc=" << dz.num_col() << endl; 00536 #endif 00537 #if defined(PRINTOUT) 00538 cout << "dr(1,1)=" << dr(1,1) << endl; 00539 #endif 00540 00541 dr(4,4) = r(); 00542 dr.fast(2,2) = r(); 00543 00544 dt.assign(m); 00545 dt.assign(s); 00546 dt.assign(dc); 00547 00548 dr *= 3.0; 00549 dr /= 3.0; 00550 00551 dr += dt; 00552 dr -= di; 00553 00554 dt = dz; 00555 00556 dz = dt.T(); 00557 di = dt.apply(neg); 00558 00559 s = dr.similarityT(mm); 00560 s = dr.similarity(mt); 00561 #if defined(PRINTOUT) 00562 cout << "sim=" << dr.similarity(v); 00563 #endif 00564 00565 dt = dr.sub(2,5); 00566 dz.sub(3,dt); 00567 00568 // 00569 // Test of HepVector 00570 // 00571 00572 HepVector vt; 00573 HepVector vp(6); 00574 HepVector vz(6,0); 00575 HepVector vi(6,1); 00576 HepVector vr(6,r); 00577 HepVector vc(vr); 00578 #if defined(PRINTOUT) 00579 cout << "vz=" << vz << "vi=" << vi << endl; 00580 #endif 00581 00582 HepVector vm(HepMatrix(6,1,r)); 00583 HepVector vs(HepSymMatrix(1,r)); 00584 HepVector vd(HepDiagMatrix(1,r)); 00585 #if defined(PRINTOUT) 00586 cout << "Nr=" << vr.num_row() << endl; 00587 #endif 00588 #if defined(PRINTOUT) 00589 cout << "vr(3)=" << vr(3) << endl; 00590 #endif 00591 00592 vr(2) = r(); 00593 00594 vi *= 2.5; 00595 vi /= 2.5; 00596 00597 vm += HepMatrix(6,1,r); 00598 vs += HepSymMatrix(1,r); 00599 vd += HepDiagMatrix(1,r); 00600 vi += vr; 00601 00602 vm -= HepMatrix(6,1,r); 00603 vs -= HepSymMatrix(1,r); 00604 vd -= HepDiagMatrix(1,r); 00605 vi -= vc; 00606 00607 vm = HepMatrix(6,1,r); 00608 vs = HepSymMatrix(1,r); 00609 vd = HepDiagMatrix(1,r); 00610 vt = vc; 00611 00612 vt = vc.apply(negv); 00613 vt = vc.sub(2,5); 00614 vc.sub(3,vt); 00615 #if defined(PRINTOUT) 00616 cout << "Normsq=" << vc.normsq() << "Norm=" << vc.norm() << endl; 00617 #endif 00618 00619 m = vc.T(); 00620 00621 swap(vc,vr); 00622 00623 vt = sr * vr; 00624 vt = dr * vr; 00625 vt = mt * vr; 00626 f = vr * m; 00627 00628 // 00629 // Test of other operators 00630 // 00631 00632 f = 3.5 * m; 00633 s = 3.5 * ss; 00634 dt = 3.5 * dr; 00635 vt = 3.5 * vr; 00636 00637 f = m * 3.6; 00638 s = ss * 3.6; 00639 dt = dr * 3.6; 00640 vt = vr * 3.6; 00641 00642 f = m / 3.6; 00643 s = ss / 3.6; 00644 dt = dr / 3.6; 00645 vt = vr / 3.6; 00646 00647 d = HepMatrix(6,3,r); 00648 e = HepMatrix(6,3,r); 00649 m = d + e; 00650 d = HepMatrix(6,6,r); 00651 m = d + sr; 00652 m = d + dr; 00653 00654 m = sr + d; 00655 m = dr + d; 00656 00657 s = sr + si; 00658 dt = dr + di; 00659 00660 s = dr + sr; 00661 s = sr + dr; 00662 00663 e = HepMatrix(6,1,r); 00664 00665 v = e + vr; 00666 00667 v = vr + e; 00668 00669 v = vr + vz; 00670 00671 00672 d = HepMatrix(6,3,r); 00673 e = HepMatrix(6,3,r); 00674 m = d - e; 00675 d = HepMatrix(6,6,r); 00676 dr = HepDiagMatrix(6,r); 00677 m = d - sr; 00678 m = d - dr; 00679 00680 m = sr - d; 00681 m = dr - d; 00682 00683 s = sr - si; 00684 dt = dr - di; 00685 00686 s = dr - sr; 00687 s = sr - dr; 00688 00689 e = HepMatrix(6,1,r); 00690 00691 v = e - vr; 00692 00693 v = vr - e; 00694 00695 v = vr - vz; 00696 00697 // 00698 // Test of Matrix inversion 00699 // 00700 00701 m = HepMatrix(6,6,r); 00702 // a = m; 00703 int inv; 00704 a = m.inverse(inv); 00705 #if defined(PRINTOUT) 00706 cout << "Inv=" << inv << endl; 00707 #endif 00708 if(inv==0) { 00709 #if defined(PRINTOUT) 00710 HepMatrix am(a * m), ma(m * a); 00711 cout << "a*m=" << am.apply(absf) << "m*a=" << ma.apply(absf) << endl; 00712 #endif 00713 } 00714 00715 v = HepVector(6,r); 00716 00717 vt = solve(m, v); 00718 #if defined(PRINTOUT) 00719 cout << "m*vt=" << m*vt << "v=" << v << endl; 00720 #endif 00721 00722 00723 ss = HepSymMatrix(6,r); 00724 // s = ss; 00725 s = ss.inverse(inv); 00726 #if defined(PRINTOUT) 00727 cout << "Inv=" << inv << endl; 00728 #endif 00729 if(inv==0) { 00730 #if defined(PRINTOUT) 00731 HepMatrix sss(s*ss), ssss(ss*s); 00732 cout << sss.apply(absf) << ssss.apply(absf) << endl; 00733 #endif 00734 } 00735 00736 #if defined(PRINTOUT) 00737 cout << "Before diag:ss=" << ss << endl; 00738 s = ss; 00739 #endif 00740 m = diagonalize(&ss); 00741 #if defined(PRINTOUT) 00742 cout << "m=" << m << endl; 00743 #endif 00744 #if defined(PRINTOUT) 00745 cout << "ss=" << ss << endl; 00746 cout << "Sim" << ss.similarity(m) << endl; 00747 HepSymMatrix diff(ss.similarity(m) - s); 00748 cout << "Diff" << diff.apply(absf) << endl; 00749 #endif 00750 00751 m = HepMatrix(6,6,r); 00752 a = qr_inverse(m); 00753 #if defined(PRINTOUT) 00754 HepMatrix am(a * m), ma(m * a); 00755 cout << am.apply(absf) << ma.apply(absf) << endl; 00756 #endif 00757 #if defined(PRINTOUT) 00758 cout << "Norm 1 =" << norm1(m) << endl; 00759 #endif 00760 #if defined(PRINTOUT) 00761 cout << "Norm 2 =" << norm(m) << endl; 00762 #endif 00763 #if defined(PRINTOUT) 00764 cout << "Norm i =" << norm_infinity(m) << endl; 00765 #endif 00766 00767 m = HepMatrix(8,6,r); 00768 #if defined(PRINTOUT) 00769 cout << "m=" << m << endl; 00770 #endif 00771 a = qr_decomp(&m); 00772 #if defined(PRINTOUT) 00773 cout << "a=" << a << "m=" << m << endl; 00774 #endif 00775 #if defined(PRINTOUT) 00776 cout << "a*m=" << a*m << endl; 00777 #endif 00778 00779 m = HepMatrix(8,8,r); 00780 v = HepVector(8,r); 00781 vt = qr_solve(m,v); 00782 #if defined(PRINTOUT) 00783 cout << "v=" << v << " m*vt=" << m*vt << endl; 00784 #endif 00785 00786 m = HepMatrix(8,8,r); 00787 a = HepMatrix(8,3,r); 00788 b = qr_solve(m,a); 00789 #if defined(PRINTOUT) 00790 cout << "v=" << a << " m*b=" << m*b << endl; 00791 #endif 00792 #if defined(PRINTOUT) 00793 cout << "Test was successful" << endl; 00794 #endif 00795 00796 // 00797 // Some useful things it can do. 00798 // 00799 // Spherisity 00800 // HepUniformRandomAB ab(-1,1); 00801 HepJamesRandom Myengine; 00802 RandFlat ab(Myengine); 00803 HepVector p[10]; 00804 int i; 00805 for(i=0;i<10;i++) { 00806 p[i] = HepVector(3,ab); 00807 p[i] *= 2.0; 00808 p[i](1) -= 1; 00809 p[i](2) -= 1; 00810 p[i](3) -= 1; 00811 } 00812 HepSymMatrix sp(3,0); 00813 for(i=0;i<10;i++) { 00814 double psq = p[i].normsq(); 00815 sp(1,1) += psq - p[i](1)*p[i](1); 00816 sp(2,2) += psq - p[i](2)*p[i](2); 00817 sp(3,3) += psq - p[i](3)*p[i](3); 00818 sp(2,1) -= p[i](2)*p[i](1); 00819 sp(3,1) -= p[i](3)*p[i](1); 00820 sp(3,2) -= p[i](3)*p[i](2); 00821 } 00822 HepMatrix spd = diagonalize(&sp); 00823 cout << "sp=" << sp << " spd=" << spd << endl; 00824 HepSymMatrix spps(7,r); 00825 HepSymMatrix sppss(spps); 00826 cout << "condition = " << condition(spps) << endl; 00827 HepMatrix sb = diagonalize(&spps); 00828 HepSymMatrix sppssdiff(sppss - spps.similarity(sb)); 00829 cout << "spps=" << spps << "sppss - sim = " << sppssdiff.apply(absf) 00830 << endl; 00831 return 0; 00832 }