home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1994 March / Source_Code_CD-ROM_Walnut_Creek_March_1994.iso / compsrcs / misc / volume34 / newmat07 / part07 < prev    next >
Encoding:
Text File  |  1993-01-10  |  67.3 KB  |  2,310 lines
  1. Newsgroups: comp.sources.misc
  2. From: robertd@kauri.vuw.ac.nz (Robert Davies)
  3. Subject: v34i113:  newmat07 - A matrix package in C++, Part07/08
  4. Message-ID: <1993Jan11.153328.2612@sparky.imd.sterling.com>
  5. X-Md4-Signature: 281209e80b5e7b57b84c4fab2ca49406
  6. Date: Mon, 11 Jan 1993 15:33:28 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: robertd@kauri.vuw.ac.nz (Robert Davies)
  10. Posting-number: Volume 34, Issue 113
  11. Archive-name: newmat07/part07
  12. Environment: C++
  13. Supersedes: newmat06: Volume 34, Issue 7-13
  14.  
  15. #! /bin/sh
  16. # This is a shell archive.  Remove anything before this line, then unpack
  17. # it by saving it into a file and typing "sh file".  To overwrite existing
  18. # files, type "sh file -c".  You can also feed this as standard input via
  19. # unshar, or by typing "sh <file", e.g..  If this archive is complete, you
  20. # will see the following message at the end:
  21. #        "End of archive 7 (of 8)."
  22. # Contents:  tmt.cxx tmt1.cxx tmt2.cxx tmt3.cxx tmt4.cxx tmt5.cxx
  23. #   tmt6.cxx tmt7.cxx tmt8.cxx tmt9.cxx tmta.cxx tmtb.cxx tmtc.cxx
  24. #   tmtd.cxx tmte.cxx tmtf.cxx tmtg.cxx tmth.cxx tmti.cxx
  25. # Wrapped by robert@kea on Sun Jan 10 23:58:30 1993
  26. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  27. if test -f 'tmt.cxx' -a "${1}" != "-c" ; then 
  28.   echo shar: Will not clobber existing file \"'tmt.cxx'\"
  29. else
  30. echo shar: Extracting \"'tmt.cxx'\" \(4633 characters\)
  31. sed "s/^X//" >'tmt.cxx' <<'END_OF_FILE'
  32. X#define WANT_STREAM
  33. X
  34. X#include "include.h"
  35. X
  36. X#include "newmat.h"
  37. X
  38. X/**************************** test program ******************************/
  39. X
  40. Xclass PrintCounter
  41. X{
  42. X   int count;
  43. X   char* s;
  44. Xpublic:
  45. X   ~PrintCounter();
  46. X   PrintCounter(char * sx) : count(0), s(sx) {}
  47. X   void operator++() { count++; }
  48. X};
  49. X
  50. X   PrintCounter PCZ("Number of non-zero matrices = ");
  51. X   PrintCounter PCN("Number of matrices tested   = ");
  52. X
  53. XPrintCounter::~PrintCounter()
  54. X{ cout << s << count << "\n"; }
  55. X
  56. X
  57. Xvoid Print(const Matrix& X)
  58. X{
  59. X   ++PCN;
  60. X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
  61. X   cout << X.Nrows() << ", ";
  62. X   cout << X.Ncols() << ")\n\n";
  63. X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
  64. X   int nr=X.Nrows(); int nc=X.Ncols();
  65. X   for (int i=1; i<=nr; i++)
  66. X   {
  67. X      for (int j=1; j<=nc; j++)  cout << X(i,j) << "\t";
  68. X      cout << "\n";
  69. X   }
  70. X   cout << flush; ++PCZ;
  71. X}
  72. X
  73. Xvoid Print(const UpperTriangularMatrix& X)
  74. X{
  75. X   ++PCN;
  76. X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
  77. X   cout << X.Nrows() << ", ";
  78. X   cout << X.Ncols() << ")\n\n";
  79. X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
  80. X   int nr=X.Nrows(); int nc=X.Ncols();
  81. X   for (int i=1; i<=nr; i++)
  82. X   {
  83. X      for (int j=1; j<i; j++) cout << "\t";
  84. X      for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
  85. X      cout << "\n";
  86. X   }
  87. X   cout << flush; ++PCZ;
  88. X}
  89. X
  90. Xvoid Print(const DiagonalMatrix& X)
  91. X{
  92. X   ++PCN;
  93. X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
  94. X   cout << X.Nrows() << ", ";
  95. X   cout << X.Ncols() << ")\n\n";
  96. X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
  97. X   int nr=X.Nrows(); int nc=X.Ncols();
  98. X   for (int i=1; i<=nr; i++)
  99. X   {
  100. X      for (int j=1; j<i; j++) cout << "\t";
  101. X      if (i<=nc) cout << X(i,i) << "\t";
  102. X      cout << "\n";
  103. X   }
  104. X   cout << flush; ++PCZ;
  105. X}
  106. X
  107. Xvoid Print(const SymmetricMatrix& X)
  108. X{
  109. X   ++PCN;
  110. X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
  111. X   cout << X.Nrows() << ", ";
  112. X   cout << X.Ncols() << ")\n\n";
  113. X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
  114. X   int nr=X.Nrows(); int nc=X.Ncols();
  115. X   for (int i=1; i<=nr; i++)
  116. X   {
  117. X      for (int j=1; j<i; j++) cout << X(j,i) << "\t";
  118. X      for (j=i; j<=nc; j++)  cout << X(i,j) << "\t";
  119. X      cout << "\n";
  120. X   }
  121. X   cout << flush; ++PCZ;
  122. X}
  123. X
  124. Xvoid Print(const LowerTriangularMatrix& X)
  125. X{
  126. X   ++PCN;
  127. X   cout << "\nMatrix type: " << (char*)X.Type() << " (";
  128. X   cout << X.Nrows() << ", ";
  129. X   cout << X.Ncols() << ")\n\n";
  130. X   if (X.IsZero()) { cout << "All elements are zero\n" << flush; return; }
  131. X   int nr=X.Nrows();
  132. X   for (int i=1; i<=nr; i++)
  133. X   {
  134. X      for (int j=1; j<=i; j++) cout << X(i,j) << "\t";
  135. X      cout << "\n";
  136. X   }
  137. X   cout << flush; ++PCZ;
  138. X}
  139. X
  140. X
  141. Xvoid Clean(Matrix& A, Real c)
  142. X{
  143. X   int nr = A.Nrows(); int nc = A.Ncols();
  144. X   for (int i=1; i<=nr; i++)
  145. X   {
  146. X      for ( int j=1; j<=nc; j++)
  147. X      { Real a = A(i,j); if ((a < c) && (a > -c)) A(i,j) = 0.0; }
  148. X   }
  149. X}
  150. X
  151. Xvoid Clean(DiagonalMatrix& A, Real c)
  152. X{
  153. X   int nr = A.Nrows();
  154. X   for (int i=1; i<=nr; i++)
  155. X   { Real a = A(i,i); if ((a < c) && (a > -c)) A(i,i) = 0.0; }
  156. X}
  157. X
  158. X
  159. Xvoid trymat1(); void trymat2(); void trymat3();
  160. Xvoid trymat4(); void trymat5(); void trymat6();
  161. Xvoid trymat7(); void trymat8(); void trymat9();
  162. Xvoid trymata(); void trymatb(); void trymatc();
  163. Xvoid trymatd(); void trymate(); void trymatf();
  164. Xvoid trymatg(); void trymath(); void trymati();
  165. X
  166. Xmain()
  167. X{
  168. X   Real* s1; Real* s2; Real* s3; Real* s4;
  169. X   cout << "\nBegin test\n";   // Forces cout to allocate memory at beginning
  170. X   { Matrix A1(40,200); s1 = A1.Store(); }
  171. X   { Matrix A1(1,1); s3 = A1.Store(); }
  172. X   {
  173. X      Tracer et("Matrix test program");
  174. X
  175. X      Matrix A(25,150);
  176. X      {
  177. X     int i;
  178. X     RowVector A(8);
  179. X     for (i=1;i<=7;i++) A(i)=0.0; A(8)=1.0;
  180. X     Print(A);
  181. X      }
  182. X      cout << "\n";
  183. X
  184. X      TestTypeAdd(); TestTypeMult(); TestTypeOrder();
  185. X
  186. X      trymat1();
  187. X      trymat2();
  188. X      trymat3();
  189. X      trymat4();
  190. X      trymat5();
  191. X      trymat6();
  192. X      trymat7();
  193. X      trymat8();
  194. X      trymat9();
  195. X      trymata();
  196. X      trymatb();
  197. X      trymatc();
  198. X      trymatd();
  199. X      trymate();
  200. X      trymatf();
  201. X      trymatg();
  202. X      trymath();
  203. X      trymati();
  204. X   }
  205. X   { Matrix A1(40,200); s2 = A1.Store(); }
  206. X   cout << "\nChecking for lost memory: "
  207. X      << (unsigned long)s1 << " " << (unsigned long)s2 << " ";
  208. X   if (s1 != s2) cout << " - error\n"; else cout << " - ok\n\n";
  209. X   { Matrix A1(1,1); s4 = A1.Store(); }
  210. X   cout << "\nChecking for lost memory: "
  211. X      << (unsigned long)s3 << " " << (unsigned long)s4 << " ";
  212. X   if (s3 != s4) cout << " - error\n"; else cout << " - ok\n\n";
  213. X#ifdef DO_FREE_CHECK
  214. X   FreeCheck::Status();
  215. X#endif 
  216. X   return 0;
  217. X}
  218. END_OF_FILE
  219. if test 4633 -ne `wc -c <'tmt.cxx'`; then
  220.     echo shar: \"'tmt.cxx'\" unpacked with wrong size!
  221. fi
  222. # end of 'tmt.cxx'
  223. fi
  224. if test -f 'tmt1.cxx' -a "${1}" != "-c" ; then 
  225.   echo shar: Will not clobber existing file \"'tmt1.cxx'\"
  226. else
  227. echo shar: Extracting \"'tmt1.cxx'\" \(2304 characters\)
  228. sed "s/^X//" >'tmt1.cxx' <<'END_OF_FILE'
  229. X
  230. X//#define WANT_STREAM
  231. X
  232. X
  233. X
  234. X#include "include.h"
  235. X
  236. X#include "newmat.h"
  237. X
  238. X/**************************** test program ******************************/
  239. X
  240. Xvoid Print(const Matrix& X);
  241. Xvoid Print(const UpperTriangularMatrix& X);
  242. Xvoid Print(const DiagonalMatrix& X);
  243. Xvoid Print(const SymmetricMatrix& X);
  244. Xvoid Print(const LowerTriangularMatrix& X);
  245. X
  246. X
  247. Xvoid trymat1()
  248. X{
  249. X//   cout << "\nFirst test of Matrix package\n\n";
  250. X   Tracer et("First test of Matrix package");
  251. X   Exception::PrintTrace(TRUE);
  252. X   {
  253. X      Tracer et1("Stage 1");
  254. X      int i,j;
  255. X
  256. X      LowerTriangularMatrix L(10);
  257. X      for (i=1;i<=10;i++) for (j=1;j<=i;j++) L(i,j)=2.0+i*i+j;
  258. X      SymmetricMatrix S(10);
  259. X      for (i=1;i<=10;i++) for (j=1;j<=i;j++) S(i,j)=i*j+1.0;
  260. X      SymmetricMatrix S1 = S / 2.0;
  261. X      S = S1 * 2.0;
  262. X      UpperTriangularMatrix U=L.t()*2.0;
  263. X      Print(LowerTriangularMatrix(L-U.t()*0.5));
  264. X      DiagonalMatrix D(10);
  265. X      for (i=1;i<=10;i++) D(i,i)=(i-4)*(i-5)*(i-6);
  266. X      Matrix M=(S+U-D+L)*(L+U-D+S);
  267. X      DiagonalMatrix DD=D*D;
  268. X      LowerTriangularMatrix LD=L*D;
  269. X      // expressions split for Turbo C
  270. X      Matrix M1 = S*L + U*L - D*L + L*L + 10.0;
  271. X      { M1 = M1 + S*U + U*U - D*U + L*U - S*D; }
  272. X      { M1 = M1 - U*D + DD - LD + S*S; }
  273. X      { M1 = M1 + U*S - D*S + L*S - 10.0; }
  274. X      M=M1-M;
  275. X      Print(M);
  276. X   }
  277. X   {
  278. X      Tracer et1("Stage 2");
  279. X      int i,j;
  280. X
  281. X      LowerTriangularMatrix L(9);
  282. X      for (i=1;i<=9;i++) for (j=1;j<=i;j++) L(i,j)=1.0+j;
  283. X      UpperTriangularMatrix U1(9);
  284. X      for (j=1;j<=9;j++) for (i=1;i<=j;i++) U1(i,j)=1.0+i;
  285. X      LowerTriangularMatrix LX(9);
  286. X      for (i=1;i<=9;i++) for (j=1;j<=i;j++) LX(i,j)=1.0+i*i;
  287. X      UpperTriangularMatrix UX(9);
  288. X      for (j=1;j<=9;j++) for (i=1;i<=j;i++) UX(i,j)=1.0+j*j;
  289. X      {
  290. X         L=L+LX/0.5;   L=L-LX*3.0;   L=LX*2.0+L;
  291. X         U1=U1+UX*2.0; U1=U1-UX*3.0; U1=UX*2.0+U1;
  292. X      }
  293. X
  294. X
  295. X      SymmetricMatrix S(9);
  296. X      for (i=1;i<=9;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
  297. X      {
  298. X         SymmetricMatrix S1 = S;
  299. X         S=S1+5.0;
  300. X         S=S-3.0;
  301. X      }
  302. X
  303. X      DiagonalMatrix D(9);
  304. X      for (i=1;i<=9;i++) D(i,i)=S(i,i);
  305. X      UpperTriangularMatrix U=L.t()*2.0;
  306. X      {
  307. X         U1=U1*2.0 - U;  Print(U1);
  308. X         L=L*2.0-D; U=U-D;
  309. X      }
  310. X      Matrix M=U+L; S=S*2.0; M=S-M; Print(M);
  311. X   }
  312. X//   cout << "\nEnd of first test\n";
  313. X}
  314. X
  315. END_OF_FILE
  316. if test 2304 -ne `wc -c <'tmt1.cxx'`; then
  317.     echo shar: \"'tmt1.cxx'\" unpacked with wrong size!
  318. fi
  319. # end of 'tmt1.cxx'
  320. fi
  321. if test -f 'tmt2.cxx' -a "${1}" != "-c" ; then 
  322.   echo shar: Will not clobber existing file \"'tmt2.cxx'\"
  323. else
  324. echo shar: Extracting \"'tmt2.cxx'\" \(2264 characters\)
  325. sed "s/^X//" >'tmt2.cxx' <<'END_OF_FILE'
  326. X
  327. X//#define WANT_STREAM
  328. X
  329. X
  330. X#include "include.h"
  331. X
  332. X#include "newmat.h"
  333. X
  334. X/**************************** test program ******************************/
  335. X
  336. Xvoid Print(const Matrix& X);
  337. Xvoid Print(const UpperTriangularMatrix& X);
  338. Xvoid Print(const DiagonalMatrix& X);
  339. Xvoid Print(const SymmetricMatrix& X);
  340. Xvoid Print(const LowerTriangularMatrix& X);
  341. X
  342. X
  343. Xvoid trymat2()
  344. X{
  345. X//   cout << "\nSecond test of Matrix package\n\n";
  346. X   Tracer et("Second test of Matrix package");
  347. X   Exception::PrintTrace(TRUE);
  348. X
  349. X   int i,j;
  350. X   Matrix M(3,5);
  351. X   for (i=1; i<=3; i++) for (j=1; j<=5; j++) M(i,j) = 100*i + j;
  352. X   Matrix X(8,10);
  353. X   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
  354. X   Matrix Y = X; Matrix Z = X;
  355. X   { X.SubMatrix(2,4,3,7) << M; }
  356. X   for (i=1; i<=3; i++) for (j=1; j<=5; j++) Y(i+1,j+2) = 100*i + j;
  357. X   Print(Matrix(X-Y));
  358. X
  359. X
  360. X   Real a[15]; Real* r = a;
  361. X   for (i=1; i<=3; i++) for (j=1; j<=5; j++) *r++ = 100*i + j;
  362. X   { Z.SubMatrix(2,4,3,7) << a; }
  363. X   Print(Matrix(Z-Y));
  364. X
  365. X   { M=33; X.SubMatrix(2,4,3,7) << M; }
  366. X   { Z.SubMatrix(2,4,3,7) = 33; }
  367. X   Print(Matrix(Z-X));
  368. X
  369. X   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
  370. X   Y = X;
  371. X   UpperTriangularMatrix U(5);
  372. X   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
  373. X   { X.SubMatrix(3,7,5,9) << U; }
  374. X   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
  375. X   for (i=1; i<=5; i++) for (j=1; j<i; j++) Y(i+2,j+4) = 0.0;
  376. X   Print(Matrix(X-Y));
  377. X   for (i=1; i<=8; i++) for (j=1; j<=10; j++) X(i,j) = 1000*i + 10*j;
  378. X   Y = X;
  379. X   for (i=1; i<=5; i++) for (j=i; j<=5; j++) U(i,j) = 100*i + j;
  380. X   { X.SubMatrix(3,7,5,9).Inject(U); }
  381. X   for (i=1; i<=5; i++) for (j=i; j<=5; j++) Y(i+2,j+4) = 100*i + j;
  382. X   Print(Matrix(X-Y));
  383. X
  384. X
  385. X   // test growing and shrinking a vector
  386. X   ColumnVector V(100);
  387. X   for (i=1;i<=100;i++) V(i) = i*i+i;
  388. X   V = V.Rows(1,50);               // to get first 50 vlaues.
  389. X
  390. X   {
  391. X      V.Release(); ColumnVector VX=V;
  392. X      V.ReDimension(100); V = 0.0; V.Rows(1,50)=VX;
  393. X   }                               // V now length 100
  394. X
  395. X#if __ZTC__
  396. X   ColumnVector VM = V; VM = 100;
  397. X#else
  398. X   M=V; M=100;                     // to make sure V will hold its values
  399. X#endif
  400. X   for (i=1;i<=50;i++) V(i) -= i*i+i;
  401. X   Print(V);
  402. X
  403. X
  404. X//   cout << "\nEnd of second test\n";
  405. X}
  406. END_OF_FILE
  407. if test 2264 -ne `wc -c <'tmt2.cxx'`; then
  408.     echo shar: \"'tmt2.cxx'\" unpacked with wrong size!
  409. fi
  410. # end of 'tmt2.cxx'
  411. fi
  412. if test -f 'tmt3.cxx' -a "${1}" != "-c" ; then 
  413.   echo shar: Will not clobber existing file \"'tmt3.cxx'\"
  414. else
  415. echo shar: Extracting \"'tmt3.cxx'\" \(2192 characters\)
  416. sed "s/^X//" >'tmt3.cxx' <<'END_OF_FILE'
  417. X
  418. X//#define WANT_STREAM
  419. X
  420. X#include "include.h"
  421. X
  422. X#include "newmat.h"
  423. X
  424. X
  425. X/**************************** test program ******************************/
  426. X
  427. Xvoid Print(const Matrix& X);
  428. Xvoid Print(const UpperTriangularMatrix& X);
  429. Xvoid Print(const DiagonalMatrix& X);
  430. Xvoid Print(const SymmetricMatrix& X);
  431. Xvoid Print(const LowerTriangularMatrix& X);
  432. X
  433. Xvoid trymat3()
  434. X{
  435. X//   cout << "\nThird test of Matrix package\n";
  436. X   Tracer et("Third test of Matrix package");
  437. X   Exception::PrintTrace(TRUE);
  438. X
  439. X   int i,j;
  440. X
  441. X   {
  442. X      Tracer et1("Stage 1");
  443. X      SymmetricMatrix S(7);
  444. X      for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
  445. X      S=-S+2.0;
  446. X
  447. X      DiagonalMatrix D(7);
  448. X      for (i=1;i<=7;i++) D(i,i)=S(i,i);
  449. X
  450. X      Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0;  M4=M4-4.0; Print(M4); }
  451. X      SymmetricMatrix S2=D; Matrix M2=S2;  { M2=-D+M2; Print(M2); }
  452. X      UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
  453. X      LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
  454. X      M2=D; M2=M2-D; Print(M2);
  455. X      for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
  456. X      U2=L2.t(); D=D.t(); S=S.t();
  457. X      M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
  458. X      M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
  459. X   }
  460. X
  461. X   {
  462. X      Tracer et1("Stage 2");
  463. X      DiagonalMatrix D(6);
  464. X      for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
  465. X      UpperTriangularMatrix U2(7); LowerTriangularMatrix L2(7);
  466. X      for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
  467. X      { U2=L2.t(); }
  468. X      DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
  469. X      Matrix M2(6,7);
  470. X      for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
  471. X      Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
  472. X      Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
  473. X      MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
  474. X      {
  475. X         D.ReDimension(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
  476. X         LowerTriangularMatrix LD(1); LD=D;
  477. X         UpperTriangularMatrix UD(1); UD=D;
  478. X         M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
  479. X         M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
  480. X         M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
  481. X         M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
  482. X      }
  483. X   }
  484. X
  485. X//   cout << "\nEnd of third test\n";
  486. X}
  487. X
  488. END_OF_FILE
  489. if test 2192 -ne `wc -c <'tmt3.cxx'`; then
  490.     echo shar: \"'tmt3.cxx'\" unpacked with wrong size!
  491. fi
  492. # end of 'tmt3.cxx'
  493. fi
  494. if test -f 'tmt4.cxx' -a "${1}" != "-c" ; then 
  495.   echo shar: Will not clobber existing file \"'tmt4.cxx'\"
  496. else
  497. echo shar: Extracting \"'tmt4.cxx'\" \(2428 characters\)
  498. sed "s/^X//" >'tmt4.cxx' <<'END_OF_FILE'
  499. X
  500. X//#define WANT_STREAM
  501. X
  502. X
  503. X#include "include.h"
  504. X
  505. X#include "newmat.h"
  506. X
  507. X
  508. X/**************************** test program ******************************/
  509. X
  510. Xvoid Print(const Matrix& X);
  511. Xvoid Print(const UpperTriangularMatrix& X);
  512. Xvoid Print(const DiagonalMatrix& X);
  513. Xvoid Print(const SymmetricMatrix& X);
  514. Xvoid Print(const LowerTriangularMatrix& X);
  515. X
  516. Xvoid trymat4()
  517. X{
  518. X//   cout << "\nFourth test of Matrix package\n";
  519. X   Tracer et("Fourth test of Matrix package");
  520. X   Exception::PrintTrace(TRUE);
  521. X
  522. X   int i,j;
  523. X
  524. X   {
  525. X      Tracer et1("Stage 1");
  526. X      Matrix M(10,10);
  527. X      UpperTriangularMatrix U(10);
  528. X      for (i=1;i<=10;i++) for (j=1;j<=10;j++) M(i,j) = 100*i+j;
  529. X      U << -M;
  530. X      Matrix X1 = M.Rows(2,4);
  531. X      Matrix Y1 = U.t().Rows(2,4);
  532. X      Matrix X = U; { Print(Matrix(X.Columns(2,4).t()-Y1)); }
  533. X      RowVector RV = M.Row(5);
  534. X      {
  535. X         X.ReDimension(3,10);
  536. X         X.Row(1) << M.Row(2); X.Row(2) << M.Row(3); X.Row(3) << M.Row(4);
  537. X         Print(Matrix(X-X1));
  538. X      }
  539. X      {
  540. X         UpperTriangularMatrix V = U.SymSubMatrix(3,5);
  541. X         Matrix MV = U.SubMatrix(3,5,3,5); { Print(Matrix(MV-V)); }
  542. X         Matrix X2 = M.t().Columns(2,4); { Print(Matrix(X2-X1.t())); }
  543. X         Matrix Y2 = U.Columns(2,4); { Print(Matrix(Y2-Y1.t())); }
  544. X         ColumnVector CV = M.t().Column(5); { Print(ColumnVector(CV-RV.t())); }
  545. X         X.ReDimension(10,3); M = M.t();
  546. X         X.Column(1) << M.Column(2); X.Column(2) << M.Column(3);
  547. X         X.Column(3) << M.Column(4);
  548. X         Print(Matrix(X-X2));
  549. X      }
  550. X   }
  551. X
  552. X   {
  553. X      Tracer et1("Stage 2");
  554. X      Matrix M; Matrix X; M.ReDimension(5,8);
  555. X      for (i=1;i<=5;i++) for (j=1;j<=8;j++) M(i,j) = 100*i+j;
  556. X      {
  557. X         X = M.Columns(5,8); M.Columns(5,8) << M.Columns(1,4);
  558. X             M.Columns(1,4) << X;
  559. X         X = M.Columns(3,4); M.Columns(3,4) << M.Columns(1,2);
  560. X             M.Columns(1,2) << X;
  561. X         X = M.Columns(7,8); M.Columns(7,8) << M.Columns(5,6);
  562. X             M.Columns(5,6) << X;
  563. X      }
  564. X      {
  565. X         X = M.Column(2); M.Column(2) = M.Column(1); M.Column(1) = X;
  566. X         X = M.Column(4); M.Column(4) = M.Column(3); M.Column(3) = X;
  567. X         X = M.Column(6); M.Column(6) = M.Column(5); M.Column(5) = X;
  568. X         X = M.Column(8); M.Column(8) = M.Column(7); M.Column(7) = X;
  569. X         X.ReDimension(5,8);
  570. X      }
  571. X      for (i=1;i<=5;i++) for (j=1;j<=8;j++) X(i,9-j) = 100*i+j;
  572. X      Print(Matrix(X-M));
  573. X   }
  574. X
  575. X//   cout << "\nEnd of fourth test\n";
  576. X}
  577. X
  578. END_OF_FILE
  579. if test 2428 -ne `wc -c <'tmt4.cxx'`; then
  580.     echo shar: \"'tmt4.cxx'\" unpacked with wrong size!
  581. fi
  582. # end of 'tmt4.cxx'
  583. fi
  584. if test -f 'tmt5.cxx' -a "${1}" != "-c" ; then 
  585.   echo shar: Will not clobber existing file \"'tmt5.cxx'\"
  586. else
  587. echo shar: Extracting \"'tmt5.cxx'\" \(1725 characters\)
  588. sed "s/^X//" >'tmt5.cxx' <<'END_OF_FILE'
  589. X
  590. X//#define WANT_STREAM
  591. X
  592. X
  593. X#include "include.h"
  594. X
  595. X#include "newmat.h"
  596. X
  597. X
  598. X/**************************** test program ******************************/
  599. X
  600. Xvoid Print(const Matrix& X);
  601. Xvoid Print(const UpperTriangularMatrix& X);
  602. Xvoid Print(const DiagonalMatrix& X);
  603. Xvoid Print(const SymmetricMatrix& X);
  604. Xvoid Print(const LowerTriangularMatrix& X);
  605. X
  606. Xvoid trymat5()
  607. X{
  608. X//   cout << "\nFifth test of Matrix package\n";
  609. X   Tracer et("Fifth test of Matrix package");
  610. X   Exception::PrintTrace(TRUE);
  611. X
  612. X   int i,j;
  613. X
  614. X   Matrix A(5,6);
  615. X   for (i=1;i<=5;i++) for (j=1;j<=6;j++) A(i,j)=1+i*j+i*i+j*j;
  616. X   ColumnVector CV(6);
  617. X   for (i=1;i<=6;i++) CV(i)=i*i+3;
  618. X   ColumnVector CV2(5); for (i=1;i<=5;i++) CV2(i)=1.0;
  619. X   ColumnVector CV1=CV;
  620. X
  621. X   {
  622. X      CV=A*CV;
  623. X      RowVector RV=CV.t(); // RowVector RV; RV=CV.t();
  624. X      RV=RV-1.0;
  625. X      CV=(RV*A).t()+A.t()*CV2; CV1=(A.t()*A)*CV1 - CV;
  626. X      Print(CV1);
  627. X   }
  628. X
  629. X   CV1.ReDimension(6);
  630. X   CV2.ReDimension(6);
  631. X   CV.ReDimension(6);
  632. X   for (i=1;i<=6;i++) { CV1(i)=i*3+1; CV2(i)=10-i; CV(i)=11+i*2; }
  633. X   ColumnVector CX=CV2-CV; { CX=CX+CV1; Print(CX); }
  634. X   Print(ColumnVector(CV1+CV2-CV));
  635. X   RowVector RV=CV.t(); RowVector RV1=CV1.t();
  636. X   RowVector R=RV-RV1; Print(RowVector(R-CV2.t()));
  637. X
  638. X// test loading of list
  639. X
  640. X   RV.ReDimension(10);
  641. X   for (i=1;i<=10;i++) RV(i) = i*i;
  642. X   RV1.ReDimension(10);
  643. X   RV1 << 1 << 4 << 9 << 16 << 25 << 36 << 49 << 64 << 81 << 100; // << 121;
  644. X   Print(RowVector(RV-RV1));
  645. X   et.ReName("Fifth test of Matrix package - at end");
  646. X   Matrix X(2,3);
  647. X   X << 11 << 12 << 13
  648. X     << 21 << 22 << 23;
  649. X   X(1,1) -= 11; X(1,2) -= 12; X(1,3) -= 13;
  650. X   X(2,1) -= 21; X(2,2) -= 22; X(2,3) -= 23;
  651. X   Print(X);
  652. X//   BandMatrix BM(7,2,2); BM << 23;
  653. X
  654. X//   cout << "\nEnd of fifth test\n";
  655. X}
  656. END_OF_FILE
  657. if test 1725 -ne `wc -c <'tmt5.cxx'`; then
  658.     echo shar: \"'tmt5.cxx'\" unpacked with wrong size!
  659. fi
  660. # end of 'tmt5.cxx'
  661. fi
  662. if test -f 'tmt6.cxx' -a "${1}" != "-c" ; then 
  663.   echo shar: Will not clobber existing file \"'tmt6.cxx'\"
  664. else
  665. echo shar: Extracting \"'tmt6.cxx'\" \(1326 characters\)
  666. sed "s/^X//" >'tmt6.cxx' <<'END_OF_FILE'
  667. X
  668. X//#define WANT_STREAM
  669. X
  670. X#include "include.h"
  671. X
  672. X#include "newmat.h"
  673. X
  674. X
  675. X/**************************** test program ******************************/
  676. X
  677. Xvoid Print(const Matrix& X);
  678. Xvoid Print(const UpperTriangularMatrix& X);
  679. Xvoid Print(const DiagonalMatrix& X);
  680. Xvoid Print(const SymmetricMatrix& X);
  681. Xvoid Print(const LowerTriangularMatrix& X);
  682. X
  683. Xvoid trymat6()
  684. X{
  685. X//   cout << "\nSixth test of Matrix package\n";
  686. X   Tracer et("Sixth test of Matrix package");
  687. X   Exception::PrintTrace(TRUE);
  688. X
  689. X   int i,j;
  690. X
  691. X
  692. X   DiagonalMatrix D(6);
  693. X   UpperTriangularMatrix U(6);
  694. X   for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*i-50; D(i,i)=i*i+i-10; }
  695. X   LowerTriangularMatrix L=(U*3.0).t();
  696. X   SymmetricMatrix S(6);
  697. X   for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
  698. X   Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S;
  699. X   Matrix M(6,6);
  700. X   for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
  701. X   {
  702. X      Tracer et1("Stage 1");
  703. X      Print(Matrix(MS+(-MS)));
  704. X      Print(Matrix((S+M)-(MS+M)));
  705. X      Print(Matrix((M+U)-(M+MU)));
  706. X      Print(Matrix((M+L)-(M+ML)));
  707. X   }
  708. X   {
  709. X      Tracer et1("Stage 2");
  710. X      Print(Matrix((M+D)-(M+MD)));
  711. X      Print(Matrix((U+D)-(MU+MD)));
  712. X      Print(Matrix((D+L)-(ML+MD)));
  713. X      Print(Matrix((-U+D)+MU-MD));
  714. X      Print(Matrix((-L+D)+ML-MD));
  715. X   }
  716. X
  717. X
  718. X
  719. X//   cout << "\nEnd of sixth test\n";
  720. X}
  721. X
  722. END_OF_FILE
  723. if test 1326 -ne `wc -c <'tmt6.cxx'`; then
  724.     echo shar: \"'tmt6.cxx'\" unpacked with wrong size!
  725. fi
  726. # end of 'tmt6.cxx'
  727. fi
  728. if test -f 'tmt7.cxx' -a "${1}" != "-c" ; then 
  729.   echo shar: Will not clobber existing file \"'tmt7.cxx'\"
  730. else
  731. echo shar: Extracting \"'tmt7.cxx'\" \(1520 characters\)
  732. sed "s/^X//" >'tmt7.cxx' <<'END_OF_FILE'
  733. X
  734. X//#define WANT_STREAM
  735. X
  736. X#include "include.h"
  737. X
  738. X#include "newmat.h"
  739. X
  740. X
  741. X/**************************** test program ******************************/
  742. X
  743. Xvoid Print(const Matrix& X);
  744. Xvoid Print(const UpperTriangularMatrix& X);
  745. Xvoid Print(const DiagonalMatrix& X);
  746. Xvoid Print(const SymmetricMatrix& X);
  747. Xvoid Print(const LowerTriangularMatrix& X);
  748. X
  749. Xvoid Clean(Matrix&, Real);
  750. X
  751. Xvoid trymat7()
  752. X{
  753. X//   cout << "\nSeventh test of Matrix package\n";
  754. X   Tracer et("Seventh test of Matrix package");
  755. X   Exception::PrintTrace(TRUE);
  756. X
  757. X   int i,j;
  758. X
  759. X
  760. X   DiagonalMatrix D(6);
  761. X   UpperTriangularMatrix U(6);
  762. X   for (i=1;i<=6;i++) { for (j=i;j<=6;j++) U(i,j)=i*i*j-50; D(i,i)=i*i+i-10; }
  763. X   LowerTriangularMatrix L=(U*3.0).t();
  764. X   SymmetricMatrix S(6);
  765. X   for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j;
  766. X   Matrix MD=D; Matrix ML=L; Matrix MU=U;
  767. X   Matrix MS=S;
  768. X   Matrix M(6,6);
  769. X   for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0;  
  770. X   {
  771. X      Tracer et1("Stage 1");
  772. X      Print(Matrix((S-M)-(MS-M)));
  773. X      Print(Matrix((-M-S)+(MS+M)));
  774. X      Print(Matrix((U-M)-(MU-M)));
  775. X   }
  776. X   {
  777. X      Tracer et1("Stage 2");
  778. X      Print(Matrix((L-M)+(M-ML)));
  779. X      Print(Matrix((D-M)+(M-MD)));
  780. X      Print(Matrix((D-S)+(MS-MD)));
  781. X      Print(Matrix((D-L)+(ML-MD)));
  782. X   }
  783. X
  784. X   { M=MU.t(); }
  785. X   LowerTriangularMatrix LY=D.i()*U.t();
  786. X   {
  787. X      Tracer et1("Stage 3");
  788. X      MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
  789. X      L=U.t();
  790. X      LY=D.i()*L; MS=D*LY-M; Clean(MS,0.00000001); Print(MS);
  791. X   }
  792. X
  793. X//   cout << "\nEnd of seventh test\n";
  794. X}
  795. END_OF_FILE
  796. if test 1520 -ne `wc -c <'tmt7.cxx'`; then
  797.     echo shar: \"'tmt7.cxx'\" unpacked with wrong size!
  798. fi
  799. # end of 'tmt7.cxx'
  800. fi
  801. if test -f 'tmt8.cxx' -a "${1}" != "-c" ; then 
  802.   echo shar: Will not clobber existing file \"'tmt8.cxx'\"
  803. else
  804. echo shar: Extracting \"'tmt8.cxx'\" \(1419 characters\)
  805. sed "s/^X//" >'tmt8.cxx' <<'END_OF_FILE'
  806. X
  807. X//#define WANT_STREAM
  808. X
  809. X
  810. X#include "include.h"
  811. X
  812. X#include "newmat.h"
  813. X
  814. X
  815. X/**************************** test program ******************************/
  816. X
  817. Xvoid Print(const Matrix& X);
  818. Xvoid Print(const UpperTriangularMatrix& X);
  819. Xvoid Print(const DiagonalMatrix& X);
  820. Xvoid Print(const SymmetricMatrix& X);
  821. Xvoid Print(const LowerTriangularMatrix& X);
  822. X
  823. Xvoid trymat8()
  824. X{
  825. X//   cout << "\nEighth test of Matrix package\n";
  826. X   Tracer et("Eighth test of Matrix package");
  827. X   Exception::PrintTrace(TRUE);
  828. X
  829. X   int i;
  830. X
  831. X
  832. X   DiagonalMatrix D(6);
  833. X   for (i=1;i<=6;i++)  D(i,i)=i*i+i-10;
  834. X   DiagonalMatrix D2=D;
  835. X   Matrix MD=D;
  836. X
  837. X   DiagonalMatrix D1(6); for (i=1;i<=6;i++) D1(i,i)=-100+i*i*i;
  838. X   Matrix MD1=D1;
  839. X   Print(Matrix(D*D1-MD*MD1));
  840. X   Print(Matrix((-D)*D1+MD*MD1));
  841. X   Print(Matrix(D*(-D1)+MD*MD1));
  842. X   DiagonalMatrix DX=D;
  843. X   {
  844. X      Tracer et1("Stage 1");
  845. X      DX=(DX+D1)*DX; Print(Matrix(DX-(MD+MD1)*MD));
  846. X      DX=D;
  847. X      DX=-DX*DX+(DX-(-D1))*((-D1)+DX);
  848. X#ifdef __GNUG__
  849. X      Matrix MX = Matrix(MD1);
  850. X      MD1=DX+(MX.t())*(MX.t()); Print(MD1);
  851. X#else
  852. X      MD1=DX+(Matrix(MD1).t())*(Matrix(MD1).t()); Print(MD1);
  853. X#endif
  854. X      DX=D; DX=DX; DX=D2-DX; Print(DiagonalMatrix(DX));
  855. X      DX=D;
  856. X   }
  857. X   {
  858. X      Tracer et1("Stage 2");
  859. X      D.Release(2);
  860. X      D1=D; D2=D;
  861. X      Print(DiagonalMatrix(D1-DX));
  862. X      Print(DiagonalMatrix(D2-DX));
  863. X      MD1=1.0;
  864. X      Print(Matrix(MD1-1.0));
  865. X   }
  866. X
  867. X//   cout << "\nEnd of eighth test\n";
  868. X}
  869. END_OF_FILE
  870. if test 1419 -ne `wc -c <'tmt8.cxx'`; then
  871.     echo shar: \"'tmt8.cxx'\" unpacked with wrong size!
  872. fi
  873. # end of 'tmt8.cxx'
  874. fi
  875. if test -f 'tmt9.cxx' -a "${1}" != "-c" ; then 
  876.   echo shar: Will not clobber existing file \"'tmt9.cxx'\"
  877. else
  878. echo shar: Extracting \"'tmt9.cxx'\" \(1577 characters\)
  879. sed "s/^X//" >'tmt9.cxx' <<'END_OF_FILE'
  880. X
  881. X//#define WANT_STREAM
  882. X
  883. X
  884. X#include "include.h"
  885. X#include "newmat.h"
  886. X
  887. X
  888. X/**************************** test program ******************************/
  889. X
  890. Xvoid Print(const Matrix& X);
  891. Xvoid Print(const UpperTriangularMatrix& X);
  892. Xvoid Print(const DiagonalMatrix& X);
  893. Xvoid Print(const SymmetricMatrix& X);
  894. Xvoid Print(const LowerTriangularMatrix& X);
  895. X
  896. Xvoid Clean(Matrix&, Real);
  897. Xvoid Clean(DiagonalMatrix&, Real);
  898. X
  899. Xvoid trymat9()
  900. X{
  901. X//   cout << "\nNinth test of Matrix package\n";
  902. X   Tracer et("Ninth test of Matrix package");
  903. X   Exception::PrintTrace(TRUE);
  904. X
  905. X
  906. X   int i; int j;
  907. X   Matrix A(7,7); Matrix X(7,3);
  908. X   for (i=1;i<=7;i++) for (j=1;j<=7;j++) A(i,j)=i*i+j+((i==j) ? 1 : 0);
  909. X   for (i=1;i<=7;i++) for (j=1;j<=3;j++) X(i,j)=i-j;
  910. X   Matrix B = A.i(); DiagonalMatrix D(7); D=1.0;
  911. X   {
  912. X      Tracer et1("Stage 1");
  913. X      Matrix Q = B*A-D; Clean(Q, 0.000000001); Print(Q);
  914. X      Q=A; Q = Q.i() * X; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
  915. X      Q=X; Q = A.i() * Q; Q = A*Q - X; Clean(Q, 0.000000001); Print(Q);
  916. X   }
  917. X   for (i=1;i<=7;i++) D(i,i)=i*i+1;
  918. X   DiagonalMatrix E(3); for (i=1;i<=3;i++) E(i,i)=i+23;
  919. X   {
  920. X      Tracer et1("Stage 2");
  921. X      Matrix DXE = D.i() * X * E;
  922. X      DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE); 
  923. X      E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1;
  924. X   }
  925. X   DiagonalMatrix F=D;
  926. X   {
  927. X      Tracer et1("Stage 3");
  928. X      F=E.i()*F; F=F*E-D; Clean(F,0.00000001); Print(F);
  929. X      F=E.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F);
  930. X   }
  931. X   { F=E; F=F.i()*D; F=F*E-D; Clean(F,0.00000001); Print(F); }
  932. X//   cout << "\nEnd of ninth test\n";
  933. X}
  934. X
  935. END_OF_FILE
  936. if test 1577 -ne `wc -c <'tmt9.cxx'`; then
  937.     echo shar: \"'tmt9.cxx'\" unpacked with wrong size!
  938. fi
  939. # end of 'tmt9.cxx'
  940. fi
  941. if test -f 'tmta.cxx' -a "${1}" != "-c" ; then 
  942.   echo shar: Will not clobber existing file \"'tmta.cxx'\"
  943. else
  944. echo shar: Extracting \"'tmta.cxx'\" \(2837 characters\)
  945. sed "s/^X//" >'tmta.cxx' <<'END_OF_FILE'
  946. X
  947. X//#define WANT_STREAM
  948. X
  949. X
  950. X#include "include.h"
  951. X
  952. X#include "newmatap.h"
  953. X
  954. X
  955. X/**************************** test program ******************************/
  956. X
  957. Xvoid Print(const Matrix& X);
  958. Xvoid Print(const UpperTriangularMatrix& X);
  959. Xvoid Print(const DiagonalMatrix& X);
  960. Xvoid Print(const SymmetricMatrix& X);
  961. Xvoid Print(const LowerTriangularMatrix& X);
  962. X
  963. Xvoid Clean(Matrix&, Real);
  964. X
  965. X
  966. Xstatic void process(const GeneralMatrix& A,
  967. X   const ColumnVector& X1, const ColumnVector& X2)
  968. X{
  969. X      Matrix B = A;
  970. X      LinearEquationSolver L(A);
  971. X      Matrix Y(4,2);
  972. X      Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
  973. X      Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
  974. X      Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
  975. X}
  976. X
  977. X
  978. X
  979. Xvoid trymata()
  980. X{
  981. X//   cout << "\nTenth test of Matrix package\n";
  982. X   Tracer et("Tenth test of Matrix package");
  983. X   Exception::PrintTrace(TRUE);
  984. X   int i; int j;
  985. X   UpperTriangularMatrix U(8);
  986. X   for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
  987. X   Matrix X(8,6);
  988. X   for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
  989. X   Matrix Y = U.i()*X; Matrix MU=U;
  990. X   Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
  991. X   Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
  992. X   UpperTriangularMatrix UX(8);
  993. X   for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
  994. X   UX(4,4)=0; UX(4,5)=0;
  995. X   UpperTriangularMatrix UY = U.i() * UX;
  996. X   { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
  997. X   LowerTriangularMatrix LY = U.t().i() * UX.t();
  998. X   { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
  999. X   DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
  1000. X   { X=D.i()*MU; }
  1001. X   { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
  1002. X   { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
  1003. X//   X=MU.t();
  1004. X//   LY=D.i()*U.t(); Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
  1005. X//   LowerTriangularMatrix L=U.t();
  1006. X//   LY=D.i()*L; Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
  1007. X   U.ReDimension(8);
  1008. X   for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
  1009. X   MU = U;
  1010. X   MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
  1011. X   MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
  1012. X
  1013. X   // test LINEQ
  1014. X   {
  1015. X      ColumnVector X1(4), X2(4);
  1016. X      X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
  1017. X      X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
  1018. X
  1019. X
  1020. X      Matrix A(4,4);
  1021. X      A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
  1022. X      A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
  1023. X      A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
  1024. X      A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
  1025. X      process(A,X1,X2);
  1026. X
  1027. X      BandMatrix B(4,1,1);  B.Inject(A);
  1028. X      process(B,X1,X2);
  1029. X
  1030. X      UpperTriangularMatrix U(4);
  1031. X      U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
  1032. X        U(2,2)=8; U(2,3)=7; U(2,4)=6;
  1033. X              U(3,3)=2; U(3,4)=7;
  1034. X                    U(4,4)=1;
  1035. X      process(U,X1,X2);
  1036. X
  1037. X      LowerTriangularMatrix L = U.t();
  1038. X      process(L,X1,X2);
  1039. X   }
  1040. X
  1041. X//   cout << "\nEnd of tenth test\n";
  1042. X}
  1043. END_OF_FILE
  1044. if test 2837 -ne `wc -c <'tmta.cxx'`; then
  1045.     echo shar: \"'tmta.cxx'\" unpacked with wrong size!
  1046. fi
  1047. # end of 'tmta.cxx'
  1048. fi
  1049. if test -f 'tmtb.cxx' -a "${1}" != "-c" ; then 
  1050.   echo shar: Will not clobber existing file \"'tmtb.cxx'\"
  1051. else
  1052. echo shar: Extracting \"'tmtb.cxx'\" \(4020 characters\)
  1053. sed "s/^X//" >'tmtb.cxx' <<'END_OF_FILE'
  1054. X
  1055. X//#define WANT_STREAM
  1056. X
  1057. X#include "include.h"
  1058. X
  1059. X#include "newmat.h"
  1060. X
  1061. X
  1062. X/**************************** test program ******************************/
  1063. X
  1064. Xvoid Print(const Matrix& X);
  1065. Xvoid Print(const UpperTriangularMatrix& X);
  1066. Xvoid Print(const DiagonalMatrix& X);
  1067. Xvoid Print(const SymmetricMatrix& X);
  1068. Xvoid Print(const LowerTriangularMatrix& X);
  1069. X
  1070. Xvoid Clean(Matrix&, Real);
  1071. X
  1072. Xvoid trymatb()
  1073. X{
  1074. X//   cout << "\nEleventh test of Matrix package\n";
  1075. X   Tracer et("Eleventh test of Matrix package");
  1076. X   Exception::PrintTrace(TRUE);
  1077. X   int i; int j;
  1078. X   RowVector RV; RV.ReDimension(10);
  1079. X   {
  1080. X      Tracer et1("Stage 1");
  1081. X      for (i=1;i<=10;i++) RV(i)=i*i-3;
  1082. X      Matrix X(1,1); X(1,1) = .25;
  1083. X      Print(RowVector(X.i() * RV - RV / .25));
  1084. X//      Print(RowVector(X.i() * Matrix(RV) - RV / .25)); // != zortech, AT&T
  1085. X      Print(RowVector(X.i() * RV - RV / .25));
  1086. X   }
  1087. X   LowerTriangularMatrix L(5); UpperTriangularMatrix U(5);
  1088. X   for (i=1; i<=5; i++) for (j=1; j<=i; j++)
  1089. X   { L(i,j) = i*i + j -2.0; U(j,i) = i*i*j+3; }
  1090. X   DiagonalMatrix D(5);
  1091. X   for (i=1; i<=5; i++) D(i,i) = i*i + i + 2;
  1092. X   Matrix M1 = -L; Matrix M2 = L-U; Matrix M3 = U*3; Matrix M4 = U-L;
  1093. X   Matrix M5 = M1 - D; M1 = D * (-3) - M3;
  1094. X   {
  1095. X      Tracer et1("Stage 2");
  1096. X      Print(Matrix((M2-M4*2)+M5*3-M1));
  1097. X      M1 = L.t(); Print(Matrix(M1.t()-L));
  1098. X      M1 = U.t(); Print(Matrix(M1.t()-U));
  1099. X   }
  1100. X   {
  1101. X      Tracer et1("Stage 3");
  1102. X      SymmetricMatrix S(5);
  1103. X      for (i=1; i<=5; i++) for (j=1; j<=i; j++) S(i,j) = i*j+i-j+5;
  1104. X      M2 = S.i() * M4; M1 = S; M3=M1*M2-M4; Clean(M3,0.00000001); Print(M3);
  1105. X      SymmetricMatrix T(5);
  1106. X      for (i=1; i<=5; i++) for (j=1; j<=i; j++) T(i,j) = i*i*j*j+i-j+5-i*j;
  1107. X      M1 = S.i() * T; M1 = S * M1; M1 = M1-T; Clean(M1,0.00000001); Print(M1);
  1108. X      ColumnVector CV(5); for (i=1; i<=5; i++) CV(i) = i*i*i+10;
  1109. X      M1 = CV * RV;
  1110. X   }
  1111. X   {
  1112. X      Tracer et1("Stage 4");
  1113. X      M4.ReDimension(5,10);
  1114. X      for (i=1; i<=5; i++) for (j=1; j<=10; j++) M4(i,j) = (i*i*i+10)*(j*j-3);
  1115. X      Print(Matrix(M1-M4));
  1116. X      M1 = L.t(); M2 = U.t(); M3 = L+U; Print(Matrix(M1-M3.t()+M2));
  1117. X   }
  1118. X//   UpperTriangularMatrix U2((const UpperTriangularMatrix&)U); // != zortech
  1119. X   UpperTriangularMatrix U2((UpperTriangularMatrix&)U);
  1120. X   {
  1121. X      Tracer et1("Stage 5");
  1122. X      Print(Matrix(U2-U));
  1123. X      M2.ReDimension(10,10);
  1124. X      for (i=1; i<=10; i++) for (j=1; j<=10; j++) M2(i,j) = (i*i*i+10)*(j*j-3);
  1125. X      D << M2; L << M2; U << M2;               // check copy into
  1126. X      Print( Matrix( (D+M2)-(L+U) ) );
  1127. X   }
  1128. X   {
  1129. X      Tracer et1("Stage 6");
  1130. X      M1.ReDimension(6,10);
  1131. X      for (i=1; i<=6; i++) for (j=1; j<=10; j++)  M1(i,j) = 100*i + j;
  1132. X      M2 = M1.SubMatrix(3,5,4,7);  M3.ReDimension(3,4);
  1133. X      for (i=3; i<=5; i++) for (j=4; j<=7; j++)   M3(i-2,j-3) = 100*i + j;
  1134. X      Print(Matrix(M2-M3));
  1135. X   }
  1136. X   int a1,a2,a3,a4;
  1137. X   {
  1138. X      Tracer et1("Stage 7");
  1139. X      int a1,a2,a3,a4;
  1140. X      a1=4; a2=9; a3=4; a4=7;
  1141. X      U.ReDimension(10);
  1142. X      for (i=1; i<=10; i++) for (j=i; j<=10; j++)  U(i,j) = 100*i + j;
  1143. X      M2 = U.SubMatrix(a1,a2,a3,a4);
  1144. X      M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
  1145. X      for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
  1146. X         M3(i-a1+1,j-a3+1) = 100*i + j;
  1147. X      Print(Matrix(M2-M3));
  1148. X   }
  1149. X   {
  1150. X      Tracer et1("Stage 8");
  1151. X      a1=3; a2=9; a3=2; a4=7;
  1152. X      U.ReDimension(10);
  1153. X      for (i=1; i<=10; i++) for (j=i; j<=10; j++)  U(i,j) = 100*i + j;
  1154. X      M2 = U.SubMatrix(a1,a2,a3,a4);
  1155. X      M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
  1156. X      for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
  1157. X         M3(i-a1+1,j-a3+1) = 100*i + j;
  1158. X      Print(Matrix(M2-M3));
  1159. X   }
  1160. X   {
  1161. X      Tracer et1("Stage 9");
  1162. X      a1=4; a2=6; a3=2; a4=5;
  1163. X      U.ReDimension(10);
  1164. X      for (i=1; i<=10; i++) for (j=i; j<=10; j++)  U(i,j) = 100*i + j;
  1165. X      M2 = U.SubMatrix(a1,a2,a3,a4);
  1166. X      M3.ReDimension(a2-a1+1,a4-a3+1); M3=0.0;
  1167. X      for (i=a1; i<=a2; i++) for (j=(i>a3) ? i : a3; j<=a4; j++)
  1168. X         M3(i-a1+1,j-a3+1) = 100*i + j;
  1169. X      Print(Matrix(M2-M3));
  1170. X   }
  1171. X
  1172. X//   cout << "\nEnd of eleventh test\n";
  1173. X}
  1174. END_OF_FILE
  1175. if test 4020 -ne `wc -c <'tmtb.cxx'`; then
  1176.     echo shar: \"'tmtb.cxx'\" unpacked with wrong size!
  1177. fi
  1178. # end of 'tmtb.cxx'
  1179. fi
  1180. if test -f 'tmtc.cxx' -a "${1}" != "-c" ; then 
  1181.   echo shar: Will not clobber existing file \"'tmtc.cxx'\"
  1182. else
  1183. echo shar: Extracting \"'tmtc.cxx'\" \(2900 characters\)
  1184. sed "s/^X//" >'tmtc.cxx' <<'END_OF_FILE'
  1185. X
  1186. X//#define WANT_STREAM
  1187. X
  1188. X
  1189. X#include "include.h"
  1190. X#include "newmat.h"
  1191. X
  1192. Xvoid Print(const Matrix& X);
  1193. Xvoid Print(const UpperTriangularMatrix& X);
  1194. Xvoid Print(const DiagonalMatrix& X);
  1195. Xvoid Print(const SymmetricMatrix& X);
  1196. Xvoid Print(const LowerTriangularMatrix& X);
  1197. X
  1198. Xvoid Clean(Matrix&, Real);
  1199. X
  1200. X
  1201. X
  1202. Xvoid trymatc()
  1203. X{
  1204. X//   cout << "\nTwelfth test of Matrix package\n";
  1205. X   Tracer et("Twelfth test of Matrix package");
  1206. X   Exception::PrintTrace(TRUE);
  1207. X   DiagonalMatrix D(15); D=1.5;
  1208. X   Matrix A(15,15);
  1209. X   int i,j;
  1210. X   for (i=1;i<=15;i++) for (j=1;j<=15;j++) A(i,j)=i*i+j-150;
  1211. X   { A = A + D; }
  1212. X   ColumnVector B(15);
  1213. X   for (i=1;i<=15;i++) B(i)=i+i*i-150.0;
  1214. X   {
  1215. X      Tracer et1("Stage 1");
  1216. X      ColumnVector B1=B;
  1217. X      B=(A*2.0).i() * B1;
  1218. X      Matrix X = A*B-B1/2.0;
  1219. X      Clean(X, 0.000000001); Print(X);
  1220. X      A.ReDimension(3,5);
  1221. X      for (i=1; i<=3; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
  1222. X
  1223. X      B = A.AsColumn()+10000;
  1224. X      RowVector R = (A+10000).AsColumn().t();
  1225. X      Print( RowVector(R-B.t()) );
  1226. X   }
  1227. X
  1228. X   {
  1229. X      Tracer et1("Stage 2");
  1230. X      B = A.AsColumn()+10000;
  1231. X      Matrix XR = (A+10000).AsMatrix(15,1).t();
  1232. X      Print( RowVector(XR-B.t()) );
  1233. X   }
  1234. X
  1235. X   {
  1236. X      Tracer et1("Stage 3");
  1237. X      B = (A.AsMatrix(15,1)+A.AsColumn())/2.0+10000;
  1238. X      Matrix MR = (A+10000).AsColumn().t();
  1239. X      Print( RowVector(MR-B.t()) );
  1240. X
  1241. X      B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
  1242. X      MR = A.AsColumn().t();
  1243. X      Print( RowVector(MR-B.t()) );
  1244. X   }
  1245. X
  1246. X   {
  1247. X      Tracer et1("Stage 4");
  1248. X      B = (A.AsMatrix(15,1)+A.AsColumn())/2.0;
  1249. X      RowVector R = A.AsColumn().t();
  1250. X      Print( RowVector(R-B.t()) );
  1251. X   }
  1252. X
  1253. X   {
  1254. X      Tracer et1("Stage 5");
  1255. X      RowVector R = (A.AsColumn()-5000).t();
  1256. X      B = ((R.t()+10000) - A.AsColumn())-5000;
  1257. X      Print( RowVector(B.t()) );
  1258. X   }
  1259. X
  1260. X   {
  1261. X      Tracer et1("Stage 6");
  1262. X      B = A.AsColumn(); ColumnVector B1 = (A+10000).AsColumn() - 10000;
  1263. X      Print(ColumnVector(B1-B));
  1264. X   }
  1265. X
  1266. X   {
  1267. X      Tracer et1("Stage 7");
  1268. X      Matrix X = B.AsMatrix(3,5); Print(Matrix(X-A));
  1269. X      for (i=1; i<=3; i++) for (j=1; j<=5; j++) B(5*(i-1)+j) -= i+100*j;
  1270. X      Print(B);
  1271. X   }
  1272. X
  1273. X   {
  1274. X      Tracer et1("Stage 8");
  1275. X      A.ReDimension(7,7); D.ReDimension(7);
  1276. X      for (i=1; i<=7; i++) for (j=1; j<=7; j++) A(i,j) = i*j*j;
  1277. X      for (i=1; i<=7; i++) D(i,i) = i;
  1278. X      UpperTriangularMatrix U; U << A;
  1279. X      Matrix X = A; for (i=1; i<=7; i++) X(i,i) = i;
  1280. X      A.Inject(D); Print(Matrix(X-A));
  1281. X      X = U; U.Inject(D); A = U; for (i=1; i<=7; i++) X(i,i) = i;
  1282. X      Print(Matrix(X-A));
  1283. X   }
  1284. X
  1285. X   {
  1286. X      Tracer et1("Stage 9");
  1287. X      A.ReDimension(7,5);
  1288. X      for (i=1; i<=7; i++) for (j=1; j<=5; j++) A(i,j) = i+100*j;
  1289. X      Matrix Y = A; Y = Y - ((const Matrix&)A); Print(Y);
  1290. X      Matrix X = A; // X.Release();
  1291. X      Y = A; Y = ((const Matrix&)X) - A; Print(Y); Y = 0.0;
  1292. X      Y = ((const Matrix&)X) - ((const Matrix&)A); Print(Y);
  1293. X   }
  1294. X
  1295. X//   cout << "\nEnd of twelfth test\n";
  1296. X}
  1297. END_OF_FILE
  1298. if test 2900 -ne `wc -c <'tmtc.cxx'`; then
  1299.     echo shar: \"'tmtc.cxx'\" unpacked with wrong size!
  1300. fi
  1301. # end of 'tmtc.cxx'
  1302. fi
  1303. if test -f 'tmtd.cxx' -a "${1}" != "-c" ; then 
  1304.   echo shar: Will not clobber existing file \"'tmtd.cxx'\"
  1305. else
  1306. echo shar: Extracting \"'tmtd.cxx'\" \(3317 characters\)
  1307. sed "s/^X//" >'tmtd.cxx' <<'END_OF_FILE'
  1308. X
  1309. X//#define WANT_STREAM
  1310. X
  1311. X#include "include.h"
  1312. X#include "newmatap.h"
  1313. X
  1314. Xvoid Print(const Matrix& X);
  1315. Xvoid Print(const UpperTriangularMatrix& X);
  1316. Xvoid Print(const DiagonalMatrix& X);
  1317. Xvoid Print(const SymmetricMatrix& X);
  1318. Xvoid Print(const LowerTriangularMatrix& X);
  1319. X
  1320. Xvoid Clean(Matrix&, Real);
  1321. X
  1322. X
  1323. X
  1324. X
  1325. Xvoid trymatd()
  1326. X{
  1327. X//   cout << "\nThirteenth test of Matrix package\n";
  1328. X   Tracer et("Thirteenth test of Matrix package");
  1329. X   Exception::PrintTrace(TRUE);
  1330. X   Matrix X(5,20);
  1331. X   int i,j;
  1332. X   for (j=1;j<=20;j++) X(1,j) = j+1;
  1333. X   for (i=2;i<=5;i++) for (j=1;j<=20; j++) X(i,j) = (long)X(i-1,j) * j % 1001;
  1334. X   SymmetricMatrix S; S << X * X.t();
  1335. X   Matrix SM = X * X.t() - S;
  1336. X   Print(SM);
  1337. X   LowerTriangularMatrix L = Cholesky(S);
  1338. X   Matrix Diff = L*L.t()-S; Clean(Diff, 0.000000001);
  1339. X   Print(Diff);
  1340. X   LowerTriangularMatrix L1(5);
  1341. X   HHDecompose(X,L1);
  1342. X   Diff = L - L1; Clean(Diff,0.000000001); Print(Diff);
  1343. X
  1344. X   ColumnVector C1(5);
  1345. X   {
  1346. X      Tracer et1("Stage 1");
  1347. X      X.ReDimension(5,5);
  1348. X      for (j=1;j<=5;j++) X(1,j) = j+1;
  1349. X      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  1350. X         X(i,j) = (long)X(i-1,j) * j % 1001;
  1351. X      for (i=1;i<=5;i++) C1(i) = i*i;
  1352. X      CroutMatrix A = X;
  1353. X      ColumnVector C2 = A.i() * C1; C1 = X.i()  * C1;
  1354. X      X = C1 - C2; Clean(X,0.000000001); Print(X); 
  1355. X   }
  1356. X
  1357. X   {
  1358. X      Tracer et1("Stage 2");
  1359. X      X.ReDimension(7,7);
  1360. X      for (j=1;j<=7;j++) X(1,j) = j+1;
  1361. X      for (i=2;i<=7;i++) for (j=1;j<=7; j++)
  1362. X         X(i,j) = (long)X(i-1,j) * j % 1001;
  1363. X      C1.ReDimension(7);
  1364. X      for (i=1;i<=7;i++) C1(i) = i*i;
  1365. X      RowVector R1 = C1.t();
  1366. X      Diff = R1 * X.i() - ( X.t().i() * R1.t() ).t(); Clean(Diff,0.000000001);
  1367. X      Print(Diff);
  1368. X   }
  1369. X
  1370. X   {
  1371. X      Tracer et1("Stage 3");
  1372. X      X.ReDimension(5,5);
  1373. X      for (j=1;j<=5;j++) X(1,j) = j+1;
  1374. X      for (i=2;i<=5;i++) for (j=1;j<=5; j++)
  1375. X         X(i,j) = (long)X(i-1,j) * j % 1001;
  1376. X      C1.ReDimension(5);
  1377. X      for (i=1;i<=5;i++) C1(i) = i*i;
  1378. X      CroutMatrix A1 = X*X;
  1379. X      ColumnVector C2 = A1.i() * C1; C1 = X.i()  * C1; C1 = X.i()  * C1;
  1380. X      X = C1 - C2; Clean(X,0.000000001); Print(X); 
  1381. X   }
  1382. X
  1383. X
  1384. X   {
  1385. X      Tracer et1("Stage 4");
  1386. X      int n = 40;
  1387. X      SymmetricBandMatrix B(n,2); B = 0.0;
  1388. X      for (i=1; i<=n; i++)
  1389. X      {
  1390. X     B(i,i) = 6;
  1391. X     if (i<=n-1) B(i,i+1) = -4;
  1392. X     if (i<=n-2) B(i,i+2) = 1;
  1393. X      }
  1394. X      B(1,1) = 5; B(n,n) = 5;
  1395. X      SymmetricMatrix A = B;
  1396. X      ColumnVector X(n); X = 0; X(1) = 1;
  1397. X      ColumnVector Y1 = A.i() * X;
  1398. X      LowerTriangularMatrix C1 = Cholesky(A);
  1399. X      ColumnVector Y2 = C1.t().i() * (C1.i() * X) - Y1;
  1400. X      Clean(Y2, 0.000000001); Print(Y2);
  1401. X  
  1402. X      LowerBandMatrix C2 = Cholesky(B);
  1403. X      Matrix M = C2 - C1; Clean(M, 0.000000001); Print(M);
  1404. X      ColumnVector Y3 = C2.t().i() * (C2.i() * X) - Y1;
  1405. X      Clean(Y3, 0.000000001); Print(Y3);
  1406. X   }
  1407. X
  1408. X   {
  1409. X      Tracer et1("Stage 5");
  1410. X      SymmetricMatrix A(4), B(4);
  1411. X      A << 5
  1412. X        << 1 << 4
  1413. X        << 2 << 1 << 6
  1414. X        << 1 << 0 << 1 << 7;
  1415. X      B << 8
  1416. X        << 1 << 5
  1417. X        << 1 << 0 << 9
  1418. X        << 2 << 1 << 0 << 6;
  1419. X      LowerTriangularMatrix AB = Cholesky(A) * Cholesky(B);
  1420. X      Matrix M = Cholesky(A) * B * Cholesky(A).t() - AB*AB.t();
  1421. X      Clean(M, 0.000000001); Print(M);
  1422. X      M = A * Cholesky(B); M = M * M.t() - A * B * A;
  1423. X      Clean(M, 0.000000001); Print(M);
  1424. X   }
  1425. X
  1426. X
  1427. X//   cout << "\nEnd of Thirteenth test\n";
  1428. X}
  1429. END_OF_FILE
  1430. if test 3317 -ne `wc -c <'tmtd.cxx'`; then
  1431.     echo shar: \"'tmtd.cxx'\" unpacked with wrong size!
  1432. fi
  1433. # end of 'tmtd.cxx'
  1434. fi
  1435. if test -f 'tmte.cxx' -a "${1}" != "-c" ; then 
  1436.   echo shar: Will not clobber existing file \"'tmte.cxx'\"
  1437. else
  1438. echo shar: Extracting \"'tmte.cxx'\" \(5726 characters\)
  1439. sed "s/^X//" >'tmte.cxx' <<'END_OF_FILE'
  1440. X
  1441. X//#define WANT_STREAM
  1442. X#define WANT_MATH
  1443. X
  1444. X#include "include.h"
  1445. X
  1446. X#include "newmatap.h"
  1447. X
  1448. Xvoid Print(const Matrix& X);
  1449. Xvoid Print(const UpperTriangularMatrix& X);
  1450. Xvoid Print(const DiagonalMatrix& X);
  1451. Xvoid Print(const SymmetricMatrix& X);
  1452. Xvoid Print(const LowerTriangularMatrix& X);
  1453. X
  1454. Xvoid Clean(Matrix&, Real);
  1455. Xvoid Clean(DiagonalMatrix&, Real);
  1456. X
  1457. X
  1458. X
  1459. X
  1460. Xvoid trymate()
  1461. X{
  1462. X//   cout << "\nFourteenth test of Matrix package\n";
  1463. X   Tracer et("Fourteenth test of Matrix package");
  1464. X   Exception::PrintTrace(TRUE);
  1465. X
  1466. X   {
  1467. X      Tracer et1("Stage 1");
  1468. X      Matrix A(8,5);
  1469. X#ifndef ATandT
  1470. X      Real   a[] =   { 22, 10,  2,  3,  7,
  1471. X               14,  7, 10,  0,  8,
  1472. X               -1, 13, -1,-11,  3,
  1473. X               -3, -2, 13, -2,  4,
  1474. X            9,  8,  1, -2,  4,
  1475. X            9,  1, -7,  5, -1,
  1476. X            2, -6,  6,  5,  1,
  1477. X            4,  5,  0, -2,  2 };
  1478. X#else
  1479. X      Real a[40];
  1480. X      a[ 0]=22; a[ 1]=10; a[ 2]= 2; a[ 3]= 3; a[ 4]= 7;
  1481. X      a[ 5]=14; a[ 6]= 7; a[ 7]=10; a[ 8]= 0; a[ 9]= 8;
  1482. X      a[10]=-1; a[11]=13; a[12]=-1; a[13]=-11;a[14]= 3;
  1483. X      a[15]=-3; a[16]=-2; a[17]=13; a[18]=-2; a[19]= 4;
  1484. X      a[20]= 9; a[21]= 8; a[22]= 1; a[23]=-2; a[24]= 4;
  1485. X      a[25]= 9; a[26]= 1; a[27]=-7; a[28]= 5; a[29]=-1;
  1486. X      a[30]= 2; a[31]=-6; a[32]= 6; a[33]= 5; a[34]= 1;
  1487. X      a[35]= 4; a[36]= 5; a[37]= 0; a[38]=-2; a[39]= 2;
  1488. X#endif
  1489. X      A << a;
  1490. X      DiagonalMatrix D; Matrix U; Matrix V;
  1491. X#ifdef ATandT
  1492. X      int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  1493. X#else
  1494. X      DiagonalMatrix I(A.Ncols());
  1495. X#endif
  1496. X      I=1.0;
  1497. X      SymmetricMatrix S1; S1 << A.t() * A;
  1498. X      SymmetricMatrix S2; S2 << A * A.t();
  1499. X      Real zero = 0.0; SVD(A+zero,D,U,V);
  1500. X      DiagonalMatrix D1; SVD(A,D1); Print(DiagonalMatrix(D-D1));
  1501. X      Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  1502. X      Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  1503. X      Matrix B = U * D * V.t() - A; Clean(B,0.000000001);Print(B);
  1504. X      D1=0.0;  SVD(A,D1,A); Print(Matrix(A-U));
  1505. X      SortDescending(D);
  1506. X      D(1) -= sqrt(1248); D(2) -= 20; D(3) -= sqrt(384);
  1507. X      Clean(D,0.000000001); Print(D);
  1508. X      Jacobi(S1, D, V);
  1509. X      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1510. X      SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
  1511. X      Clean(D,0.000000001); Print(D);
  1512. X      Jacobi(S2, D, V);
  1513. X      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1514. X      SortDescending(D); D(1)-=1248; D(2)-=400; D(3)-=384;
  1515. X      Clean(D,0.000000001); Print(D);
  1516. X
  1517. X      EigenValues(S1, D, V);
  1518. X      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1519. X      D(5)-=1248; D(4)-=400; D(3)-=384;
  1520. X      Clean(D,0.000000001); Print(D);
  1521. X      EigenValues(S2, D, V);
  1522. X      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1523. X      D(8)-=1248; D(7)-=400; D(6)-=384;
  1524. X      Clean(D,0.000000001); Print(D);
  1525. X
  1526. X      EigenValues(S1, D);
  1527. X      D(5)-=1248; D(4)-=400; D(3)-=384;
  1528. X      Clean(D,0.000000001); Print(D);
  1529. X      EigenValues(S2, D);
  1530. X      D(8)-=1248; D(7)-=400; D(6)-=384;
  1531. X      Clean(D,0.000000001); Print(D);
  1532. X   }
  1533. X   {
  1534. X      Tracer et1("Stage 2");
  1535. X      Matrix A(20,21);
  1536. X      for (int i=1; i<=20; i++) for (int j=1; j<=21; j++)
  1537. X      { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 21-i; else A(i,j) = -1; }
  1538. X      A = A.t();
  1539. X      SymmetricMatrix S1; S1 << A.t() * A;
  1540. X      SymmetricMatrix S2; S2 << A * A.t();
  1541. X      DiagonalMatrix D; Matrix U; Matrix V;
  1542. X#ifdef ATandT
  1543. X      int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  1544. X#else
  1545. X      DiagonalMatrix I(A.Ncols());
  1546. X#endif
  1547. X      I=1.0;
  1548. X      SVD(A,D,U,V);
  1549. X      Matrix SU = U.t() * U - I;    Clean(SU,0.000000001); Print(SU);
  1550. X      Matrix SV = V.t() * V - I;    Clean(SV,0.000000001); Print(SV);
  1551. X      Matrix B = U * D * V.t() - A; Clean(B,0.000000001);  Print(B);
  1552. X      for (i=1; i<=20; i++)  D(i) -= sqrt((22-i)*(21-i));
  1553. X      Clean(D,0.000000001); Print(D);
  1554. X      Jacobi(S1, D, V);
  1555. X      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1556. X      SortDescending(D);
  1557. X      for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
  1558. X      Clean(D,0.000000001); Print(D);
  1559. X      Jacobi(S2, D, V);
  1560. X      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1561. X      SortDescending(D);
  1562. X      for (i=1; i<=20; i++)  D(i) -= (22-i)*(21-i);
  1563. X      Clean(D,0.000000001); Print(D);
  1564. X
  1565. X      EigenValues(S1, D, V);
  1566. X      V = S1 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1567. X      for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
  1568. X      Clean(D,0.000000001); Print(D);
  1569. X      EigenValues(S2, D, V);
  1570. X      V = S2 - V * D * V.t(); Clean(V,0.000000001); Print(V);
  1571. X      for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
  1572. X      Clean(D,0.000000001); Print(D);
  1573. X
  1574. X      EigenValues(S1, D);
  1575. X      for (i=1; i<=20; i++)  D(i) -= (i+1)*i;
  1576. X      Clean(D,0.000000001); Print(D);
  1577. X      EigenValues(S2, D);
  1578. X      for (i=2; i<=21; i++)  D(i) -= (i-1)*i;
  1579. X      Clean(D,0.000000001); Print(D);
  1580. X   }
  1581. X
  1582. X   {
  1583. X      Tracer et1("Stage 3");
  1584. X      Matrix A(30,30);
  1585. X      for (int i=1; i<=30; i++) for (int j=1; j<=30; j++)
  1586. X      { if (i>j) A(i,j) = 0; else if (i==j) A(i,j) = 1; else A(i,j) = -1; }
  1587. X      Real d1 = A.LogDeterminant().Value();
  1588. X      DiagonalMatrix D; Matrix U; Matrix V;
  1589. X#ifdef ATandT
  1590. X      int anc = A.Ncols(); DiagonalMatrix I(anc);     // AT&T 2.1 bug
  1591. X#else
  1592. X      DiagonalMatrix I(A.Ncols());
  1593. X#endif
  1594. X      I=1.0;
  1595. X      SVD(A,D,U,V);
  1596. X      Matrix SU = U.t() * U - I; Clean(SU,0.000000001); Print(SU);
  1597. X      Matrix SV = V.t() * V - I; Clean(SV,0.000000001); Print(SV);
  1598. X      Real d2 = D.LogDeterminant().Value();
  1599. X      Matrix B = U * D * V.t() - A; Clean(B,0.000000001); Print(B);
  1600. X      SortDescending(D);  // Print(D);
  1601. X      Real d3 = D.LogDeterminant().Value();
  1602. X      ColumnVector Test(3);
  1603. X      Test(1) = d1 - 1; Test(2) = d2 - 1; Test(3) = d3 - 1;
  1604. X      Clean(Test,0.00000001); Print(Test); // only 8 decimal figures
  1605. X   }
  1606. X
  1607. X//   cout << "\nEnd of Fourteenth test\n";
  1608. X}
  1609. END_OF_FILE
  1610. if test 5726 -ne `wc -c <'tmte.cxx'`; then
  1611.     echo shar: \"'tmte.cxx'\" unpacked with wrong size!
  1612. fi
  1613. # end of 'tmte.cxx'
  1614. fi
  1615. if test -f 'tmtf.cxx' -a "${1}" != "-c" ; then 
  1616.   echo shar: Will not clobber existing file \"'tmtf.cxx'\"
  1617. else
  1618. echo shar: Extracting \"'tmtf.cxx'\" \(3322 characters\)
  1619. sed "s/^X//" >'tmtf.cxx' <<'END_OF_FILE'
  1620. X
  1621. X//#define WANT_STREAM
  1622. X#define WANT_MATH
  1623. X
  1624. X#include "include.h"
  1625. X
  1626. X#include "newmatap.h"
  1627. X
  1628. Xvoid Print(const Matrix& X);
  1629. Xvoid Print(const UpperTriangularMatrix& X);
  1630. Xvoid Print(const DiagonalMatrix& X);
  1631. Xvoid Print(const SymmetricMatrix& X);
  1632. Xvoid Print(const LowerTriangularMatrix& X);
  1633. X
  1634. Xvoid Clean(Matrix&, Real);
  1635. X
  1636. X
  1637. Xstatic void SlowFT(const ColumnVector& a, const ColumnVector&b,
  1638. X   ColumnVector& x, ColumnVector& y)
  1639. X{
  1640. X   int n = a.Nrows();
  1641. X   x.ReDimension(n); y.ReDimension(n);
  1642. X   Real f = 6.2831853071795864769/n;
  1643. X   for (int j=1; j<=n; j++)
  1644. X   {
  1645. X      Real sumx = 0.0; Real sumy = 0.0;
  1646. X      for (int k=1; k<=n; k++)
  1647. X      {
  1648. X     Real theta = - (j-1) * (k-1) * f;
  1649. X     Real c = cos(theta);  Real s = sin(theta);
  1650. X     sumx += c * a(k) - s * b(k);
  1651. X     sumy += s * a(k) + c * b(k);
  1652. X      }
  1653. X      x(j) = sumx; y(j) = sumy;
  1654. X   }
  1655. X}
  1656. X
  1657. Xstatic void test(int n)
  1658. X{
  1659. X   Tracer et("Test FFT");
  1660. X
  1661. X   ColumnVector A(n), B(n), X, Y;
  1662. X   for (int i=0; i<n; i++)
  1663. X   {
  1664. X      Real f = 6.2831853071795864769*i/n;
  1665. X      A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  1666. X      B.element(i) = fabs(0.25 * cos(31.0 * f)) + (Real)i/(Real)n;
  1667. X   }
  1668. X   FFT(A, B, X, Y);
  1669. X   FFTI(X, Y, X, Y);
  1670. X   X = X - A; Y = Y - B;
  1671. X   Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
  1672. X}
  1673. X
  1674. Xstatic void test1(int n)
  1675. X{
  1676. X   Tracer et("Test RealFFT");
  1677. X
  1678. X   ColumnVector A(n), B(n), X, Y;
  1679. X   for (int i=0; i<n; i++)
  1680. X   {
  1681. X      Real f = 6.2831853071795864769*i/n;
  1682. X      A.element(i) = fabs(sin(7.0*f) + 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  1683. X   }
  1684. X   B = 0.0;
  1685. X   FFT(A, B, X, Y);
  1686. X   int n2 = n/2+1;
  1687. X   ColumnVector X1,Y1,X2,Y2;
  1688. X   RealFFT(A, X1, Y1);
  1689. X   X2 = X1 - X.Rows(1,n2); Y2 = Y1 - Y.Rows(1,n2);
  1690. X   Clean(X2,0.000000001); Clean(Y2,0.000000001); Print(X2); Print(Y2);
  1691. X   RealFFTI(X1,Y1,B);
  1692. X   B = A - B;
  1693. X   Clean(B,0.000000001); Print(B);
  1694. X}
  1695. X
  1696. Xstatic void test2(int n)
  1697. X{
  1698. X   Tracer et("cf FFT and slow FT");
  1699. X
  1700. X   ColumnVector A(n), B(n), X, Y, X1, Y1;
  1701. X   for (int i=0; i<n; i++)
  1702. X   {
  1703. X      Real f = 6.2831853071795864769*i/n;
  1704. X      A.element(i) = fabs(sin(7.0*f) - 0.5 * cos(19.0 * f)) + (Real)i/(Real)n;
  1705. X      B.element(i) = fabs(0.25 * cos(31.0 * f)) - (Real)i/(Real)n;
  1706. X   }
  1707. X   FFT(A, B, X, Y);
  1708. X   SlowFT(A, B, X1, Y1);
  1709. X   X = X - X1; Y = Y - Y1;
  1710. X   Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y);
  1711. X}
  1712. X
  1713. Xvoid trymatf()
  1714. X{
  1715. X   Tracer et("Fifteenth test of Matrix package");
  1716. X   Exception::PrintTrace(TRUE);
  1717. X//   cout << "\nFifteenth test of Matrix package\n";
  1718. X//   cout << "\n";
  1719. X
  1720. X   ColumnVector A(12), B(12);
  1721. X   for (int i = 1; i <=12; i++)
  1722. X   {
  1723. X      Real i1 = i - 1;
  1724. X      A(i) = .7
  1725. X       + .2 * cos(6.2831853071795864769 * 4.0 * i1 / 12)
  1726. X       + .3 * sin(6.2831853071795864769 * 3.0 * i1 / 12);
  1727. X      B(i) = .9
  1728. X       + .5 * sin(6.2831853071795864769 * 2.0 * i1 / 12)
  1729. X       + .4 * cos(6.2831853071795864769 * 1.0 * i1 / 12);
  1730. X   }
  1731. X   FFT(A, B, A, B);
  1732. X   A(1) -= 8.4; A(3) -= 3.0; A(5) -= 1.2; A(9) -= 1.2; A(11) += 3.0;
  1733. X   B(1) -= 10.8; B(2) -= 2.4; B(4) += 1.8; B(10) -= 1.8; B(12) -= 2.4;
  1734. X   Clean(A,0.000000001); Clean(B,0.000000001); Print(A); Print(B);
  1735. X
  1736. X
  1737. X
  1738. X   test(2048);
  1739. X   test(2000);
  1740. X   test(27*81);
  1741. X   test(2310);
  1742. X   test(49*49);
  1743. X   test(13*13*13);
  1744. X   test(43*47);
  1745. X   test1(98);
  1746. X   test1(100);
  1747. X   test1(2048);
  1748. X   test1(2000);
  1749. X   test1(35*35*2);
  1750. X   test2(13);
  1751. X   test2(12);
  1752. X   test2(9);
  1753. X   test2(16);
  1754. X
  1755. X//   cout << "\nEnd of Fifteenth test\n";
  1756. X}
  1757. END_OF_FILE
  1758. if test 3322 -ne `wc -c <'tmtf.cxx'`; then
  1759.     echo shar: \"'tmtf.cxx'\" unpacked with wrong size!
  1760. fi
  1761. # end of 'tmtf.cxx'
  1762. fi
  1763. if test -f 'tmtg.cxx' -a "${1}" != "-c" ; then 
  1764.   echo shar: Will not clobber existing file \"'tmtg.cxx'\"
  1765. else
  1766. echo shar: Extracting \"'tmtg.cxx'\" \(3855 characters\)
  1767. sed "s/^X//" >'tmtg.cxx' <<'END_OF_FILE'
  1768. X
  1769. X//#define WANT_STREAM
  1770. X
  1771. X#include "include.h"
  1772. X
  1773. X#include "newmatap.h"
  1774. X
  1775. Xvoid Print(const Matrix& X);
  1776. Xvoid Print(const UpperTriangularMatrix& X);
  1777. Xvoid Print(const DiagonalMatrix& X);
  1778. Xvoid Print(const SymmetricMatrix& X);
  1779. Xvoid Print(const LowerTriangularMatrix& X);
  1780. X
  1781. Xvoid Clean(Matrix&, Real);
  1782. X
  1783. X
  1784. X
  1785. X
  1786. Xvoid trymatg()
  1787. X{
  1788. X//   cout << "\nSixteenth test of Matrix package\n";
  1789. X//   cout << "\n";
  1790. X   Tracer et("Sixteenth test of Matrix package");
  1791. X   Exception::PrintTrace(TRUE);
  1792. X
  1793. X   int i,j;
  1794. X   Matrix M(4,7);
  1795. X   for (i=1; i<=4; i++) for (j=1; j<=7; j++) M(i,j) = 100 * i + j;
  1796. X   ColumnVector CV = M.AsColumn();
  1797. X   {
  1798. X      Tracer et1("Stage 1");
  1799. X      RowVector test(7);
  1800. X      test(1) = SumSquare(M);
  1801. X      test(2) = SumSquare(CV);
  1802. X      test(3) = SumSquare(CV.t());
  1803. X      test(4) = SumSquare(CV.AsDiagonal());
  1804. X      test(5) = SumSquare(M.AsColumn());
  1805. X      test(6) = Matrix(CV.t()*CV)(1,1);
  1806. X      test(7) = (CV.t()*CV).AsScalar();
  1807. X      test = test - 2156560; Print(test);
  1808. X   }
  1809. X
  1810. X   UpperTriangularMatrix U(6);
  1811. X   for (i=1; i<=6; i++) for (j=i; j<=6; j++) U(i,j) = i + (i-j) * (i-j);
  1812. X   M = U; DiagonalMatrix D; D << U;
  1813. X   LowerTriangularMatrix L = U.t(); SymmetricMatrix S; S << (L+U)/2.0;
  1814. X   {
  1815. X      Tracer et1("Stage 2");
  1816. X      RowVector test(6);
  1817. X      test(1) = U.Trace();
  1818. X      test(2) = L.Trace();
  1819. X      test(3) = D.Trace();
  1820. X      test(4) = S.Trace();
  1821. X      test(5) = M.Trace();
  1822. X      test(6) = ((L+U)/2.0).Trace();
  1823. X      test = test - 21; Print(test);
  1824. X      test.ReDimension(5);
  1825. X      test(1) = LogDeterminant(U).Value();
  1826. X      test(2) = LogDeterminant(L).Value();
  1827. X      test(3) = LogDeterminant(D).Value();
  1828. X      test(4) = LogDeterminant(D).Value();
  1829. X      test(5) =  LogDeterminant((L+D)/2.0).Value(); // (!=Glockenspiel)
  1830. X      test = test - 720; Clean(test,0.000000001); Print(test);
  1831. X   }
  1832. X
  1833. X   {
  1834. X      Tracer et1("Stage 3");
  1835. X      S << L*U; M = S;
  1836. X      RowVector test(4);
  1837. X      test(1) = LogDeterminant(S).Value();
  1838. X      test(2) = LogDeterminant(M).Value();
  1839. X      test(3) = LogDeterminant(L*U).Value();           // (!=Glockenspiel)
  1840. X      test(4) = LogDeterminant(Matrix(L*L)).Value();   // (!=Glockenspiel)
  1841. X      test = test - 720.0 * 720.0; Clean(test,0.000000001); Print(test);
  1842. X   }
  1843. X
  1844. X   {
  1845. X      Tracer et1("Stage 4");
  1846. X      M = S * S;
  1847. X      Matrix SX = S;
  1848. X      RowVector test(3);
  1849. X      test(1) = SumSquare(S);
  1850. X      test(2) = SumSquare(SX);
  1851. X      test(3) = Trace(M);
  1852. X      test = test - 3925961; Print(test);
  1853. X   }
  1854. X   {
  1855. X      Tracer et1("Stage 5");
  1856. X      SymmetricMatrix SM(10), SN(10);
  1857. X      for (i=1; i<=10; i++) for (j=i; j<=10; j++)
  1858. X      {
  1859. X         SM(i,j) = 1.5 * i - j; SN(i,j) = SM(i,j) * SM(i,j);
  1860. X         if (SM(i,j) < 0.0)   SN(i,j) = - SN(i,j);
  1861. X      }
  1862. X      Matrix M = SM, N = SN; RowVector test(4);
  1863. X      test(1) = SumAbsoluteValue(SN);
  1864. X      test(2) = SumAbsoluteValue(-SN);
  1865. X      test(3) = SumAbsoluteValue(N);
  1866. X      test(4) = SumAbsoluteValue(-N);
  1867. X      test = test - 1168.75; Print(test);
  1868. X      test(1) = MaximumAbsoluteValue(SM);
  1869. X      test(2) = MaximumAbsoluteValue(-SM);
  1870. X      test(3) = MaximumAbsoluteValue(M);
  1871. X      test(4) = MaximumAbsoluteValue(-M);
  1872. X      test = test - 8.5; Print(test);
  1873. X   }
  1874. X   {
  1875. X      Tracer et1("Stage 6");
  1876. X      Matrix M(15,20); Real value = 0.0;
  1877. X      for (i=1; i<=15; i++) { for (j=1; j<=20; j++) M(i,j) = 1.5 * i - j; }
  1878. X      for (i=1; i<=20; i++)
  1879. X      { Real v = SumAbsoluteValue(M.Column(i)); if (value<v) value = v; }
  1880. X      RowVector test(3);
  1881. X      test(1) = value;
  1882. X      test(2) = Norm1(M);
  1883. X      test(3) = NormInfinity(M.t());
  1884. X      test = test - 165; Print(test);
  1885. X      test.ReDimension(5);
  1886. X      ColumnVector CV = M.AsColumn(); M = CV.t();
  1887. X      test(1) = Norm1(CV.t());
  1888. X      test(2) = MaximumAbsoluteValue(M);
  1889. X      test(3) = NormInfinity(CV);
  1890. X      test(4) = Norm1(M);
  1891. X      test(5) = NormInfinity(M.t());
  1892. X      test = test - 21.5; Print(test);
  1893. X   }   
  1894. X
  1895. X//   cout << "\nEnd of Sixteenth test\n";
  1896. X}
  1897. END_OF_FILE
  1898. if test 3855 -ne `wc -c <'tmtg.cxx'`; then
  1899.     echo shar: \"'tmtg.cxx'\" unpacked with wrong size!
  1900. fi
  1901. # end of 'tmtg.cxx'
  1902. fi
  1903. if test -f 'tmth.cxx' -a "${1}" != "-c" ; then 
  1904.   echo shar: Will not clobber existing file \"'tmth.cxx'\"
  1905. else
  1906. echo shar: Extracting \"'tmth.cxx'\" \(6032 characters\)
  1907. sed "s/^X//" >'tmth.cxx' <<'END_OF_FILE'
  1908. X
  1909. X//#define WANT_STREAM
  1910. X
  1911. X#include "include.h"
  1912. X
  1913. X#include "newmatap.h"
  1914. X//#include "newmatio.h"
  1915. X
  1916. Xvoid Print(const Matrix& X);
  1917. Xvoid Print(const UpperTriangularMatrix& X);
  1918. Xvoid Print(const DiagonalMatrix& X);
  1919. Xvoid Print(const SymmetricMatrix& X);
  1920. Xvoid Print(const LowerTriangularMatrix& X);
  1921. X
  1922. Xvoid Clean(Matrix&, Real);
  1923. X
  1924. X
  1925. X
  1926. X
  1927. Xvoid trymath()
  1928. X{
  1929. X//   cout << "\nSeventeenth test of Matrix package\n";
  1930. X//   cout << "\n";
  1931. X   Tracer et("Seventeenth test of Matrix package");
  1932. X   Exception::PrintTrace(TRUE);
  1933. X
  1934. X
  1935. X   {
  1936. X      Tracer et1("Stage 1");
  1937. X      int i, j;
  1938. X      BandMatrix B(8,3,1);
  1939. X      for (i=1; i<=8; i++) for (j=-3; j<=1; j++)
  1940. X     { int k = i+j; if (k>0 && k<=8) B(i,k) = i + k/64.0; }
  1941. X
  1942. X      DiagonalMatrix I(8); I = 1; BandMatrix B1; B1 = I;
  1943. X      UpperTriangularMatrix UT = I; Print(Matrix(B1-UT));
  1944. X      Print(Matrix(B * B - B * 2 + I - (B - I) * (B - I)));
  1945. X      Matrix A = B; BandMatrix C; C = B;
  1946. X      Print(Matrix(B * A - C * 2 + I - (A - I) * (B - I)));
  1947. X
  1948. X      C.ReDimension(8,2,2); C = 0.0; C.Inject(B);
  1949. X      Matrix X = A.t();
  1950. X      B1.ReDimension(8,2,2); B1.Inject(X); Print(Matrix(C.t()-B1));
  1951. X
  1952. X      Matrix BI = B.i(); A = A.i()-BI; Clean(A,0.000000001); Print(A);
  1953. X
  1954. X
  1955. X      UpperBandMatrix U; U << B; LowerBandMatrix L; L << B;
  1956. X      DiagonalMatrix D; D << B;
  1957. X      Print( Matrix(L + (U - D - B)) );
  1958. X
  1959. X
  1960. X
  1961. X      for (i=1; i<=8; i++)  A.Column(i) << B.Column(i);
  1962. X      Print(Matrix(A-B));
  1963. X   }
  1964. X   {
  1965. X      Tracer et1("Stage 2");
  1966. X      BandMatrix A(7,2,2);
  1967. X      int i,j;
  1968. X      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  1969. X      {
  1970. X     int k=i-j; if (k<0) k = -k;
  1971. X     if (k==0) A(i,j)=6;
  1972. X     else if (k==1) A(i,j) = -4;
  1973. X     else if (k==2) A(i,j) = 1;
  1974. X     A(1,1) = A(7,7) = 5;
  1975. X      }
  1976. X      DiagonalMatrix D(7); D = 0.0; A = A - D;
  1977. X      BandLUMatrix B(A); Matrix M = A;
  1978. X      ColumnVector V(3);
  1979. X      V(1) = LogDeterminant(B).Value();
  1980. X      V(2) = LogDeterminant(A).Value();
  1981. X      V(3) = LogDeterminant(M).Value();
  1982. X      V = V / 64 - 1; Clean(V,0.000000001); Print(V);
  1983. X      ColumnVector X(7);
  1984. X#ifdef ATandT
  1985. X      Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  1986. X#else
  1987. X      Real a[] = {1,2,3,4,5,6,7};
  1988. X#endif
  1989. X      X << a;
  1990. X      M = (M.i()*X).t() - (B.i()*X).t() * 2.0 + (A.i()*X).t();
  1991. X      Clean(M,0.000000001); Print(M);
  1992. X
  1993. X
  1994. X      BandMatrix P(80,2,5); ColumnVector CX(80);
  1995. X      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  1996. X      { int d = i-j; if (d<=2 && d>=-5) P(i,j) = i + j/100.0; }
  1997. X      for (i=1; i<=80; i++)  CX(i) = i*100.0;
  1998. X      Matrix MP = P;
  1999. X      ColumnVector V1 = P.i() * CX; ColumnVector V2 = MP.i() * CX;
  2000. X      V = V1 - V2; Clean(V,0.000000001); Print(V);
  2001. X
  2002. X      V1 = P * V1; V2 = MP * V2; V = V1 - V2; Clean(V,0.000000001); Print(V);
  2003. X      RowVector XX(1);
  2004. X      XX = LogDeterminant(P).Value() / LogDeterminant(MP).Value() - 1.0;
  2005. X      Clean(XX,0.000000001); Print(XX);
  2006. X
  2007. X      LowerBandMatrix LP(80,5);
  2008. X      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  2009. X      { int d = i-j; if (d<=5 && d>=0) LP(i,j) = i + j/100.0; }
  2010. X      MP = LP;
  2011. X      XX.ReDimension(4);
  2012. X      XX(1) = LogDeterminant(LP).Value();
  2013. X      XX(2) = LogDeterminant(MP).Value();
  2014. X      V1 = LP.i() * CX; V2 = MP.i() * CX;
  2015. X      V = V1 - V2; Clean(V,0.000000001); Print(V);
  2016. X
  2017. X      UpperBandMatrix UP(80,4);
  2018. X      for (i=1; i<=80; i++) for (j=1; j<=80; j++)
  2019. X      { int d = i-j; if (d<=0 && d>=-4) UP(i,j) = i + j/100.0; }
  2020. X      MP = UP;
  2021. X      XX(3) = LogDeterminant(UP).Value();
  2022. X      XX(4) = LogDeterminant(MP).Value();
  2023. X      V1 = UP.i() * CX; V2 = MP.i() * CX;
  2024. X      V = V1 - V2; Clean(V,0.000000001); Print(V);
  2025. X      XX = XX / SumAbsoluteValue(XX) - .25; Clean(XX,0.000000001); Print(XX);
  2026. X   }
  2027. X   {
  2028. X      Tracer et1("Stage 3");
  2029. X      SymmetricBandMatrix SA(8,5);
  2030. X      int i,j;
  2031. X      for (i=1; i<=8; i++) for (j=1; j<=8; j++)
  2032. X     if (i-j<=5 && 0<=i-j) SA(i,j) =i + j/128.0;
  2033. X      DiagonalMatrix D(8); D = 10; SA = SA + D;
  2034. X
  2035. X      Matrix MA1(8,8); Matrix MA2(8,8);
  2036. X      for (i=1; i<=8; i++)
  2037. X     { MA1.Column(i) << SA.Column(i); MA2.Row(i) << SA.Row(i); }
  2038. X      Print(Matrix(MA1-MA2));
  2039. X
  2040. X      D = 10; SA = SA.t() + D; MA1 = MA1 + D;
  2041. X      Print(Matrix(MA1-SA));
  2042. X
  2043. X      UpperBandMatrix UB(8,3); LowerBandMatrix LB(8,4);
  2044. X      D << SA; UB << SA; LB << SA;
  2045. X      SA = SA * 5.0; D = D * 5.0; LB = LB * 5.0; UB = UB * 5.0;
  2046. X      BandMatrix B = LB - D + UB - SA; Print(Matrix(B));
  2047. X
  2048. X      SymmetricBandMatrix A(7,2); A = 100.0;
  2049. X      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  2050. X      {
  2051. X     int k=i-j;
  2052. X     if (k==0) A(i,j)=6;
  2053. X     else if (k==1) A(i,j) = -4;
  2054. X     else if (k==2) A(i,j) = 1;
  2055. X     A(1,1) = A(7,7) = 5;
  2056. X      }
  2057. X      BandLUMatrix C(A); Matrix M = A;
  2058. X      ColumnVector X(7);
  2059. X      X(1) = LogDeterminant(C).Value() - 64;
  2060. X      X(2) = LogDeterminant(A).Value() - 64;
  2061. X      X(3) = LogDeterminant(M).Value() - 64;
  2062. X      X(4) = SumSquare(M) - SumSquare(A);
  2063. X      X(5) = SumAbsoluteValue(M) - SumAbsoluteValue(A);
  2064. X      X(6) = MaximumAbsoluteValue(M) - MaximumAbsoluteValue(A);
  2065. X      X(7) = Trace(M) - Trace(A);
  2066. X      Clean(X,0.000000001); Print(X);
  2067. X
  2068. X#ifdef ATandT
  2069. X      Real a[7]; a[0]=1; a[1]=2; a[2]=3; a[3]=4; a[4]=5; a[5]=6; a[6]=7;
  2070. X#else
  2071. X      Real a[] = {1,2,3,4,5,6,7};
  2072. X#endif
  2073. X      X << a;
  2074. X      X = M.i()*X - C.i()*X * 2 + A.i()*X;
  2075. X      Clean(X,0.000000001); Print(X);
  2076. X
  2077. X
  2078. X      LB << A; UB << A; D << A;
  2079. X      BandMatrix XA = LB + (UB - D);
  2080. X      Print(Matrix(XA - A));
  2081. X
  2082. X      for (i=1; i<=7; i++) for (j=1; j<=7; j++)
  2083. X      {
  2084. X     int k=i-j;
  2085. X     if (k==0) A(i,j)=6;
  2086. X     else if (k==1) A(i,j) = -4;
  2087. X     else if (k==2) A(i,j) = 1;
  2088. X     A(1,1) = A(7,7) = 5;
  2089. X      }
  2090. X      D = 1;
  2091. X
  2092. X      M = LB.i() * LB - D; Clean(M,0.000000001); Print(M);
  2093. X      M = UB.i() * UB - D; Clean(M,0.000000001); Print(M);
  2094. X      M = XA.i() * XA - D; Clean(M,0.000000001); Print(M);
  2095. X#ifdef __GNUG__
  2096. X      Matrix MUB = UB; Matrix MLB = LB;
  2097. X      M = LB.i() * UB - LB.i() * MUB; Clean(M,0.000000001); Print(M);
  2098. X      M = UB.i() * LB - UB.i() * MLB; Clean(M,0.000000001); Print(M);
  2099. X#else
  2100. X      M = LB.i() * UB - LB.i() * Matrix(UB); Clean(M,0.000000001); Print(M);
  2101. X      M = UB.i() * LB - UB.i() * Matrix(LB); Clean(M,0.000000001); Print(M);
  2102. X#endif
  2103. X   }
  2104. X
  2105. X
  2106. X//   cout << "\nEnd of Seventeenth test\n";
  2107. X}
  2108. END_OF_FILE
  2109. if test 6032 -ne `wc -c <'tmth.cxx'`; then
  2110.     echo shar: \"'tmth.cxx'\" unpacked with wrong size!
  2111. fi
  2112. # end of 'tmth.cxx'
  2113. fi
  2114. if test -f 'tmti.cxx' -a "${1}" != "-c" ; then 
  2115.   echo shar: Will not clobber existing file \"'tmti.cxx'\"
  2116. else
  2117. echo shar: Extracting \"'tmti.cxx'\" \(5259 characters\)
  2118. sed "s/^X//" >'tmti.cxx' <<'END_OF_FILE'
  2119. X
  2120. X//#define WANT_STREAM
  2121. X
  2122. X#include "include.h"
  2123. X
  2124. X#include "newmatap.h"
  2125. X//#include "newmatio.h"
  2126. X
  2127. Xvoid Print(const Matrix& X);
  2128. Xvoid Print(const UpperTriangularMatrix& X);
  2129. Xvoid Print(const DiagonalMatrix& X);
  2130. Xvoid Print(const SymmetricMatrix& X);
  2131. Xvoid Print(const LowerTriangularMatrix& X);
  2132. X
  2133. Xvoid Clean(Matrix&, Real);
  2134. X
  2135. Xvoid WillNotConverge()
  2136. X{
  2137. X   Matrix A(10,10);
  2138. X   Throw(ConvergenceException(A));
  2139. X}
  2140. X
  2141. Xvoid trymati()
  2142. X{
  2143. X   Tracer et("Eighteenth test of Matrix package");
  2144. X   Exception::PrintTrace(TRUE);
  2145. X   ProgramException::SetAction(0);           // turn off error messages
  2146. X   DataException::SetAction(0);
  2147. X   ConvergenceException::SetAction(0);
  2148. X   ColumnVector checks(14); checks = 1.0; checks(1) = 0.0;
  2149. X
  2150. X   Try { WillNotConverge(); }
  2151. X   Catch(ConvergenceException) { checks(2) = 0; }
  2152. X   Catch(InternalException) { checks(1) = 1; }
  2153. X   Catch(ProgramException) { checks(1) = 1; }
  2154. X   Catch(DataException) { checks(1) = 1; }
  2155. X   Catch(SpaceException) { checks(1) = 1; }
  2156. X   CatchAndThrow;
  2157. X
  2158. X   Try { Matrix M(10,10); SymmetricMatrix S = M; } 
  2159. X   Catch(ConvergenceException) { checks(1) = 1; }
  2160. X   Catch(InternalException) { checks(1) = 1; }
  2161. X   Catch(ProgramException) { checks(3) = 0; }
  2162. X   Catch(DataException) { checks(1) = 1; }
  2163. X   Catch(SpaceException) { checks(1) = 1; }
  2164. X   CatchAndThrow;
  2165. X
  2166. X   Try { Matrix M(10,10); M(10,11) = 2.0; } 
  2167. X   Catch(ConvergenceException) { checks(1) = 1; }
  2168. X   Catch(InternalException) { checks(1) = 1; }
  2169. X   Catch(ProgramException) { checks(4) = 0; }
  2170. X   Catch(DataException) { checks(1) = 1; }
  2171. X   Catch(SpaceException) { checks(1) = 1; }
  2172. X   CatchAndThrow;
  2173. X
  2174. X   Try { Matrix M(10,10); M = 0.0; M = M.i(); } 
  2175. X   Catch(ConvergenceException) { checks(1) = 1; }
  2176. X   Catch(InternalException) { checks(1) = 1; }
  2177. X   Catch(ProgramException) { checks(1) = 1; }
  2178. X   Catch(DataException) { checks(5) = 0; }
  2179. X   Catch(SpaceException) { checks(1) = 1; }
  2180. X   CatchAndThrow;
  2181. X
  2182. X   Try { ColumnVector A(30), B(50);  A = 5; B = 3; FFT(A,B,A,B); } 
  2183. X   Catch(ConvergenceException) { checks(1) = 1; }
  2184. X   Catch(InternalException) { checks(1) = 1; }
  2185. X   Catch(ProgramException) { checks(6) = 0; }
  2186. X   Catch(DataException) { checks(1) = 1; }
  2187. X   Catch(SpaceException) { checks(1) = 1; }
  2188. X   CatchAndThrow;
  2189. X
  2190. X   Try
  2191. X   {
  2192. X      ColumnVector A(30); A = 5; Matrix At = A.t();
  2193. X      DiagonalMatrix D;
  2194. X      SVD(At,D);
  2195. X   } 
  2196. X   Catch(ConvergenceException) { checks(1) = 1; }
  2197. X   Catch(InternalException) { checks(1) = 1; }
  2198. X   Catch(ProgramException) { checks(6) = 0; }
  2199. X   Catch(DataException) { checks(1) = 1; }
  2200. X   Catch(SpaceException) { checks(1) = 1; }
  2201. X   CatchAndThrow;
  2202. X
  2203. X   Try { BandMatrix X(10,3,4); X(1,10) = 4.0; } 
  2204. X   Catch(ConvergenceException) { checks(1) = 1; }
  2205. X   Catch(InternalException) { checks(1) = 1; }
  2206. X   Catch(ProgramException) { checks(7) = 0; }
  2207. X   Catch(DataException) { checks(1) = 1; }
  2208. X   Catch(SpaceException) { checks(1) = 1; }
  2209. X   CatchAndThrow;
  2210. X
  2211. X   Try
  2212. X   {
  2213. X      SymmetricMatrix S(10); S = 5;
  2214. X      LowerTriangularMatrix L = Cholesky(S);
  2215. X   } 
  2216. X   Catch(ConvergenceException) { checks(1) = 1; }
  2217. X   Catch(InternalException) { checks(1) = 1; }
  2218. X   Catch(ProgramException) { checks(1) = 1; }
  2219. X   Catch(DataException) { checks(8) = 0; }
  2220. X   Catch(SpaceException) { checks(1) = 1; }
  2221. X   CatchAndThrow;
  2222. X
  2223. X   Try { BandMatrix M(10,3,5); M = 0.0; Matrix XM = M.i(); } 
  2224. X   Catch(ConvergenceException) { checks(1) = 1; }
  2225. X   Catch(InternalException) { checks(1) = 1; }
  2226. X   Catch(ProgramException) { checks(1) = 1; }
  2227. X   Catch(DataException) { checks(9) = 0; }
  2228. X   Catch(SpaceException) { checks(1) = 1; }
  2229. X   CatchAndThrow;
  2230. X
  2231. X
  2232. X   Try { ColumnVector X(10); ColumnVector Y; X = 5; X = X - Y; }
  2233. X   Catch(ConvergenceException) { checks(1) = 1; }
  2234. X   Catch(InternalException) { checks(1) = 1; }
  2235. X   Catch(ProgramException) { checks(10) = 0; }
  2236. X   Catch(DataException) { checks(1) = 1; }
  2237. X   Catch(SpaceException) { checks(1) = 1; }
  2238. X   CatchAndThrow;
  2239. X
  2240. X   Try
  2241. X   {
  2242. X      UpperTriangularMatrix U(3); RowVector RV(3); RV = 10;
  2243. X      U.Row(2) = RV;
  2244. X   }
  2245. X   Catch(ConvergenceException) { checks(1) = 1; }
  2246. X   Catch(InternalException) { checks(1) = 1; }
  2247. X   Catch(ProgramException) { checks(11) = 0; }
  2248. X   Catch(DataException) { checks(1) = 1; }
  2249. X   Catch(SpaceException) { checks(1) = 1; }
  2250. X   CatchAndThrow;
  2251. X
  2252. X   Try { DiagonalMatrix D(3); D << 12 << 13 << 14 << 15; }
  2253. X   Catch(ConvergenceException) { checks(1) = 1; }
  2254. X   Catch(InternalException) { checks(1) = 1; }
  2255. X   Catch(ProgramException) { checks(12) = 0; }
  2256. X   Catch(DataException) { checks(1) = 1; }
  2257. X   Catch(SpaceException) { checks(1) = 1; }
  2258. X   CatchAndThrow;
  2259. X
  2260. X   Try { ColumnVector D(3); D << 12 << 13; D << 1 << 2 << 3; }
  2261. X   Catch(ConvergenceException) { checks(1) = 1; }
  2262. X   Catch(InternalException) { checks(1) = 1; }
  2263. X   Catch(ProgramException) { checks(13) = 0; }
  2264. X   Catch(DataException) { checks(1) = 1; }
  2265. X   Catch(SpaceException) { checks(1) = 1; }
  2266. X   CatchAndThrow;
  2267. X
  2268. X
  2269. X   Try {  { ColumnVector D(3); D << 12 << 13; }  }
  2270. X   Catch(ConvergenceException) { checks(1) = 1; }
  2271. X   Catch(InternalException) { checks(1) = 1; }
  2272. X   Catch(ProgramException) { checks(14) = 0; }
  2273. X   Catch(DataException) { checks(1) = 1; }
  2274. X   Catch(SpaceException) { checks(1) = 1; }
  2275. X   CatchAndThrow;
  2276. X
  2277. X
  2278. X   ProgramException::SetAction(1);           // restore error messages
  2279. X   DataException::SetAction(1);
  2280. X   ConvergenceException::SetAction(1);
  2281. X
  2282. X   Print(checks);
  2283. X
  2284. X}
  2285. END_OF_FILE
  2286. if test 5259 -ne `wc -c <'tmti.cxx'`; then
  2287.     echo shar: \"'tmti.cxx'\" unpacked with wrong size!
  2288. fi
  2289. # end of 'tmti.cxx'
  2290. fi
  2291. echo shar: End of archive 7 \(of 8\).
  2292. cp /dev/null ark7isdone
  2293. MISSING=""
  2294. for I in 1 2 3 4 5 6 7 8 ; do
  2295.     if test ! -f ark${I}isdone ; then
  2296.     MISSING="${MISSING} ${I}"
  2297.     fi
  2298. done
  2299. if test "${MISSING}" = "" ; then
  2300.     echo You have unpacked all 8 archives.
  2301.     rm -f ark[1-9]isdone
  2302. else
  2303.     echo You still need to unpack the following archives:
  2304.     echo "        " ${MISSING}
  2305. fi
  2306. ##  End of shell archive.
  2307. exit 0
  2308.  
  2309. exit 0 # Just in case...
  2310.