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