home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1994 March / Source_Code_CD-ROM_Walnut_Creek_March_1994.iso / compsrcs / misc / volume34 / newmat07 / part05 < prev    next >
Encoding:
Text File  |  1993-01-10  |  58.0 KB  |  2,173 lines

  1. Newsgroups: comp.sources.misc
  2. From: robertd@kauri.vuw.ac.nz (Robert Davies)
  3. Subject: v34i111:  newmat07 - A matrix package in C++, Part05/08
  4. Message-ID: <1993Jan11.153215.2460@sparky.imd.sterling.com>
  5. X-Md4-Signature: 95c3bfe3e634dc03ab1b546ccee8137a
  6. Date: Mon, 11 Jan 1993 15:32:15 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 111
  11. Archive-name: newmat07/part05
  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 5 (of 8)."
  22. # Contents:  newmat5.cxx newmat6.cxx newmat7.cxx newmat8.cxx
  23. #   newmat9.cxx newmatrm.cxx
  24. # Wrapped by robert@kea on Sun Jan 10 23:58:07 1993
  25. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  26. if test -f 'newmat5.cxx' -a "${1}" != "-c" ; then 
  27.   echo shar: Will not clobber existing file \"'newmat5.cxx'\"
  28. else
  29. echo shar: Extracting \"'newmat5.cxx'\" \(11567 characters\)
  30. sed "s/^X//" >'newmat5.cxx' <<'END_OF_FILE'
  31. X//$$ newmat5.cxx         Transpose, evaluate etc
  32. X
  33. X// Copyright (C) 1991,2,3: R B Davies
  34. X
  35. X#include "include.h"
  36. X
  37. X#include "newmat.h"
  38. X#include "newmatrc.h"
  39. X
  40. X//#define REPORT { static ExeCounter ExeCount(__LINE__,5); ExeCount++; }
  41. X
  42. X#define REPORT {}
  43. X
  44. X
  45. X/************************ carry out operations ******************************/
  46. X
  47. X
  48. XGeneralMatrix* GeneralMatrix::Transpose(TransposedMatrix* tm, MatrixType mt)
  49. X{
  50. X   GeneralMatrix* gm1;
  51. X
  52. X   if (Compare(Type().t(),mt))
  53. X   {
  54. X      REPORT
  55. X      gm1 = mt.New(ncols,nrows,tm);
  56. X      for (int i=0; i<ncols; i++)
  57. X      {
  58. X     MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
  59. X         MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
  60. X      }
  61. X   }
  62. X   else
  63. X   {
  64. X      REPORT
  65. X      gm1 = mt.New(ncols,nrows,tm);
  66. X      MatrixRow mr(this, LoadOnEntry);
  67. X      MatrixCol mc(gm1, StoreOnExit+DirectPart);
  68. X      int i = nrows;
  69. X      while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
  70. X   }
  71. X   tDelete(); gm1->ReleaseAndDelete(); return gm1;
  72. X}
  73. X
  74. XGeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  75. X{ REPORT  return Evaluate(mt); }
  76. X
  77. X
  78. XGeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
  79. X{ REPORT return Evaluate(mt); }
  80. X
  81. XBoolean GeneralMatrix::IsZero() const
  82. X{
  83. X   REPORT
  84. X   Real* s=store; int i=storage;
  85. X   while (i--) { if (*s++) return FALSE; }
  86. X   return TRUE;
  87. X}
  88. X
  89. XGeneralMatrix* ColumnVector::Transpose(TransposedMatrix*, MatrixType mt)
  90. X{
  91. X   REPORT
  92. X   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  93. X   gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
  94. X   return BorrowStore(gmx,mt);
  95. X}
  96. X
  97. XGeneralMatrix* RowVector::Transpose(TransposedMatrix*, MatrixType mt)
  98. X{
  99. X   REPORT
  100. X   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  101. X   gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
  102. X   return BorrowStore(gmx,mt);
  103. X}
  104. X
  105. XGeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
  106. X{
  107. X   if (Compare(this->Type(),mt)) { REPORT return this; }
  108. X   REPORT
  109. X   GeneralMatrix* gmx = mt.New(nrows,ncols,this);
  110. X   MatrixRow mr(this, LoadOnEntry);
  111. X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  112. X   int i=nrows;
  113. X   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  114. X   tDelete(); gmx->ReleaseAndDelete(); return gmx;
  115. X}
  116. X
  117. XGeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
  118. X{
  119. X   if (Compare(cgm->Type(),mt))
  120. X   {
  121. X      REPORT
  122. X#ifdef TEMPS_DESTROYED_QUICKLY
  123. X      GeneralMatrix* gmx = (GeneralMatrix*)cgm; delete this; return gmx;
  124. X#else
  125. X      return (GeneralMatrix*)cgm;
  126. X#endif
  127. X   }
  128. X   REPORT
  129. X   GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols(),this);
  130. X   MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
  131. X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  132. X   int i=cgm->Nrows();
  133. X   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
  134. X   gmx->ReleaseAndDelete();
  135. X#ifdef TEMPS_DESTROYED_QUICKLY
  136. X   delete this;
  137. X#endif
  138. X   return gmx; // no tDelete
  139. X}
  140. X
  141. XGeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
  142. X{
  143. X   gm=((BaseMatrix*&)bm)->Evaluate();
  144. X   int nr=gm->Nrows(); int nc=gm->Ncols();
  145. X   Compare(gm->Type().AddEqualEl(),mt);
  146. X   if (!(mt==gm->Type()))
  147. X   {
  148. X      REPORT
  149. X      GeneralMatrix* gmx = mt.New(nr,nc,this);
  150. X      MatrixRow mr(gm, LoadOnEntry); 
  151. X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  152. X      while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
  153. X      gmx->ReleaseAndDelete(); gm->tDelete();
  154. X#ifdef TEMPS_DESTROYED_QUICKLY
  155. X      delete this;
  156. X#endif
  157. X      return gmx;
  158. X   }
  159. X   else if (gm->reuse())
  160. X   {
  161. X      REPORT gm->Add(f);
  162. X#ifdef TEMPS_DESTROYED_QUICKLY
  163. X      GeneralMatrix* gmx = gm; delete this; return gmx;
  164. X#else
  165. X      return gm;
  166. X#endif
  167. X   }
  168. X   else
  169. X   {
  170. X      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
  171. X      gmy->ReleaseAndDelete(); gmy->Add(gm,f);
  172. X#ifdef TEMPS_DESTROYED_QUICKLY
  173. X      delete this;
  174. X#endif
  175. X      return gmy;
  176. X   }
  177. X}
  178. X
  179. XGeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
  180. X{
  181. X   gm=((BaseMatrix*&)bm)->Evaluate();
  182. X   int nr=gm->Nrows(); int nc=gm->Ncols();
  183. X   if (Compare(gm->Type(),mt))
  184. X   {
  185. X      if (gm->reuse())
  186. X      {
  187. X         REPORT gm->Multiply(f);
  188. X#ifdef TEMPS_DESTROYED_QUICKLY
  189. X         GeneralMatrix* gmx = gm; delete this; return gmx;
  190. X#else
  191. X         return gm;
  192. X#endif
  193. X      }
  194. X      else
  195. X      {
  196. X         REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  197. X         gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
  198. X#ifdef TEMPS_DESTROYED_QUICKLY
  199. X         delete this;
  200. X#endif
  201. X         return gmx;
  202. X      }
  203. X   }
  204. X   else
  205. X   {
  206. X      REPORT
  207. X      GeneralMatrix* gmx = mt.New(nr,nc,this);
  208. X      MatrixRow mr(gm, LoadOnEntry); 
  209. X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  210. X      while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
  211. X      gmx->ReleaseAndDelete(); gm->tDelete();
  212. X#ifdef TEMPS_DESTROYED_QUICKLY
  213. X      delete this;
  214. X#endif
  215. X      return gmx;
  216. X   }
  217. X}
  218. X
  219. XGeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
  220. X{
  221. X   gm=((BaseMatrix*&)bm)->Evaluate();
  222. X   int nr=gm->Nrows(); int nc=gm->Ncols();
  223. X   if (Compare(gm->Type(),mt))
  224. X   {
  225. X      if (gm->reuse())
  226. X      {
  227. X         REPORT gm->Negate();
  228. X#ifdef TEMPS_DESTROYED_QUICKLY
  229. X         GeneralMatrix* gmx = gm; delete this; return gmx;
  230. X#else
  231. X         return gm;
  232. X#endif
  233. X      }
  234. X      else
  235. X      {
  236. X         REPORT
  237. X         GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
  238. X         gmx->ReleaseAndDelete(); gmx->Negate(gm);
  239. X#ifdef TEMPS_DESTROYED_QUICKLY
  240. X         delete this;
  241. X#endif
  242. X         return gmx;
  243. X      }
  244. X   }
  245. X   else
  246. X   {
  247. X      REPORT
  248. X      GeneralMatrix* gmx = mt.New(nr,nc,this);
  249. X      MatrixRow mr(gm, LoadOnEntry); 
  250. X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  251. X      while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
  252. X      gmx->ReleaseAndDelete(); gm->tDelete();
  253. X#ifdef TEMPS_DESTROYED_QUICKLY
  254. X      delete this;
  255. X#endif
  256. X      return gmx;
  257. X   }
  258. X}   
  259. X
  260. XGeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
  261. X{
  262. X   REPORT
  263. X   gm=((BaseMatrix*&)bm)->Evaluate();
  264. X   Compare(gm->Type().t(),mt);
  265. X   GeneralMatrix* gmx=gm->Transpose(this, mt);
  266. X#ifdef TEMPS_DESTROYED_QUICKLY
  267. X   delete this;
  268. X#endif
  269. X   return gmx;
  270. X}
  271. X   
  272. XGeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
  273. X{
  274. X   gm = ((BaseMatrix*&)bm)->Evaluate();
  275. X   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
  276. X   gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
  277. X#ifdef TEMPS_DESTROYED_QUICKLY
  278. X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  279. X#else
  280. X   return gm->BorrowStore(gmx,mt);
  281. X#endif
  282. X}
  283. X
  284. XGeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
  285. X{
  286. X   gm = ((BaseMatrix*&)bm)->Evaluate();
  287. X   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
  288. X   gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
  289. X#ifdef TEMPS_DESTROYED_QUICKLY
  290. X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  291. X#else
  292. X   return gm->BorrowStore(gmx,mt);
  293. X#endif
  294. X}
  295. X
  296. XGeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
  297. X{
  298. X   gm = ((BaseMatrix*&)bm)->Evaluate();
  299. X   GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
  300. X   gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
  301. X#ifdef TEMPS_DESTROYED_QUICKLY
  302. X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  303. X#else
  304. X   return gm->BorrowStore(gmx,mt);
  305. X#endif
  306. X}
  307. X
  308. XGeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
  309. X{
  310. X   Tracer tr("MatedMatrix::Evaluate");
  311. X   gm = ((BaseMatrix*&)bm)->Evaluate();
  312. X   GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
  313. X   gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
  314. X   if (nr*nc != gmx->storage)
  315. X      Throw(IncompatibleDimensionsException());
  316. X#ifdef TEMPS_DESTROYED_QUICKLY
  317. X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
  318. X#else
  319. X   return gm->BorrowStore(gmx,mt);
  320. X#endif
  321. X}
  322. X
  323. XGeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
  324. X{
  325. X   REPORT
  326. X   Tracer tr("SubMatrix(evaluate)");
  327. X   gm = ((BaseMatrix*&)bm)->Evaluate();
  328. X   if (row_number < 0) row_number = gm->Nrows();
  329. X   if (col_number < 0) col_number = gm->Ncols();
  330. X   if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  331. X      Throw(SubMatrixDimensionException());
  332. X   if (IsSym) Compare(gm->Type().ssub(), mt); 
  333. X   else Compare(gm->Type().sub(), mt);
  334. X   GeneralMatrix* gmx = mt.New(row_number, col_number, this);
  335. X   int i = row_number;
  336. X   MatrixRow mr(gm, LoadOnEntry, row_skip); 
  337. X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
  338. X   MatrixRowCol sub;
  339. X   while (i--)
  340. X   {
  341. X      mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  342. X      mrx.Copy(sub); mrx.Next(); mr.Next();
  343. X   }
  344. X   gmx->ReleaseAndDelete(); gm->tDelete();
  345. X#ifdef TEMPS_DESTROYED_QUICKLY
  346. X   delete this;
  347. X#endif
  348. X   return gmx;
  349. X}   
  350. X
  351. X
  352. XGeneralMatrix* ReturnMatrixX::Evaluate(MatrixType mt)
  353. X{
  354. X#ifdef TEMPS_DESTROYED_QUICKLY
  355. X   GeneralMatrix* gmx = gm; delete this; return gmx->Evaluate(mt);
  356. X#else
  357. X   return gm->Evaluate(mt);
  358. X#endif
  359. X}
  360. X
  361. X
  362. Xvoid GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
  363. X{
  364. X   REPORT
  365. X   Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  366. X   while (i--)
  367. X   { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
  368. X   i = storage & 3; while (i--) *s++ = *s1++ + f;
  369. X}
  370. X   
  371. Xvoid GeneralMatrix::Add(Real f)
  372. X{
  373. X   REPORT
  374. X   Real* s=store; int i=(storage >> 2);
  375. X   while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
  376. X   i = storage & 3; while (i--) *s++ += f;
  377. X}
  378. X   
  379. Xvoid GeneralMatrix::Negate(GeneralMatrix* gm1)
  380. X{
  381. X   // change sign of elements
  382. X   REPORT
  383. X   Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
  384. X   while (i--)
  385. X   { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
  386. X   i = storage & 3; while(i--) *s++ = -(*s1++);
  387. X}
  388. X   
  389. Xvoid GeneralMatrix::Negate()
  390. X{
  391. X   REPORT
  392. X   Real* s=store; int i=(storage >> 2);
  393. X   while (i--)
  394. X   { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
  395. X   i = storage & 3; while(i--) { *s = -(*s); s++; }
  396. X}
  397. X   
  398. Xvoid GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
  399. X{
  400. X   REPORT
  401. X   Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
  402. X   while (i--)
  403. X   { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
  404. X   i = storage & 3; while (i--) *s++ = *s1++ * f;
  405. X}
  406. X   
  407. Xvoid GeneralMatrix::Multiply(Real f)
  408. X{
  409. X   REPORT
  410. X   Real* s=store; int i=(storage >> 2);
  411. X   while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
  412. X   i = storage & 3; while (i--) *s++ *= f;
  413. X}
  414. X   
  415. X
  416. X/************************ MatrixInput routines ****************************/
  417. X
  418. Xint MatrixInput::n = 0;       // number values still to be read
  419. XReal* MatrixInput::r;         // pointer to last thing read
  420. Xint MatrixInput::depth = 0;   // number of objects of this class in existence
  421. X
  422. XMatrixInput MatrixInput::operator<<(Real f)
  423. X{
  424. X   if (!(n--))
  425. X   { depth=0;  n=0; Throw(ProgramException("List of values too long")); }
  426. X   *(++r) = f;
  427. X   return MatrixInput();
  428. X}
  429. X
  430. X
  431. XMatrixInput BandMatrix::operator<<(Real)
  432. X{
  433. X   Throw(ProgramException("Cannot use list read with a BandMatrix"));
  434. X   return MatrixInput(); 
  435. X}
  436. X
  437. Xvoid BandMatrix::operator<<(const Real*)
  438. X{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
  439. X
  440. XMatrixInput GeneralMatrix::operator<<(Real f)
  441. X{
  442. X   if (MatrixInput::n)
  443. X   {
  444. X      MatrixInput::depth=0;            // so we can recover
  445. X      MatrixInput::n=0;                // so we can recover
  446. X      Throw(ProgramException("A list of values was too short"));
  447. X   }
  448. X   MatrixInput::n = Storage();
  449. X   if (MatrixInput::n<=0)
  450. X      Throw(ProgramException("Loading data to zero dimension matrix"));
  451. X   MatrixInput::r = Store(); *(MatrixInput::r) = f; MatrixInput::n--;
  452. X   return MatrixInput(); 
  453. X}
  454. X
  455. XMatrixInput::~MatrixInput()
  456. X{
  457. X   if (depth && --depth == 0 && n != 0)
  458. X   {
  459. X      depth = 0; n = 0;                // so we can recover
  460. X      Throw(ProgramException("A list of values was too short"));
  461. X   }
  462. X}
  463. X
  464. END_OF_FILE
  465. if test 11567 -ne `wc -c <'newmat5.cxx'`; then
  466.     echo shar: \"'newmat5.cxx'\" unpacked with wrong size!
  467. fi
  468. # end of 'newmat5.cxx'
  469. fi
  470. if test -f 'newmat6.cxx' -a "${1}" != "-c" ; then 
  471.   echo shar: Will not clobber existing file \"'newmat6.cxx'\"
  472. else
  473. echo shar: Extracting \"'newmat6.cxx'\" \(13261 characters\)
  474. sed "s/^X//" >'newmat6.cxx' <<'END_OF_FILE'
  475. X//$$ newmat6.cxx            Operators, element access, submatrices
  476. X
  477. X// Copyright (C) 1991,2,3: R B Davies
  478. X
  479. X#include "include.h"
  480. X
  481. X#include "newmat.h"
  482. X#include "newmatrc.h"
  483. X
  484. X
  485. X//#define REPORT { static ExeCounter ExeCount(__LINE__,6); ExeCount++; }
  486. X
  487. X#define REPORT {}
  488. X
  489. X/*************************** general utilities *************************/
  490. X
  491. Xstatic int tristore(int n)                      // els in triangular matrix
  492. X{ return (n*(n+1))/2; }
  493. X
  494. X
  495. X/****************************** operators *******************************/
  496. X
  497. XReal& Matrix::operator()(int m, int n)
  498. X{
  499. X   if (m<=0 || m>nrows || n<=0 || n>ncols)
  500. X      Throw(IndexException(m,n,*this));
  501. X   return store[(m-1)*ncols+n-1];
  502. X}
  503. X
  504. XReal& SymmetricMatrix::operator()(int m, int n)
  505. X{
  506. X   if (m<=0 || n<=0 || m>nrows || n>ncols)
  507. X      Throw(IndexException(m,n,*this));
  508. X   if (m>=n) return store[tristore(m-1)+n-1];
  509. X   else return store[tristore(n-1)+m-1];
  510. X}
  511. X
  512. XReal& UpperTriangularMatrix::operator()(int m, int n)
  513. X{
  514. X   if (m<=0 || n<m || n>ncols)
  515. X      Throw(IndexException(m,n,*this));
  516. X   return store[(m-1)*ncols+n-1-tristore(m-1)];
  517. X}
  518. X
  519. XReal& LowerTriangularMatrix::operator()(int m, int n)
  520. X{
  521. X   if (n<=0 || m<n || m>nrows)
  522. X      Throw(IndexException(m,n,*this));
  523. X   return store[tristore(m-1)+n-1];
  524. X}
  525. X
  526. XReal& DiagonalMatrix::operator()(int m, int n)
  527. X{
  528. X   if (n<=0 || m!=n || m>nrows || n>ncols)
  529. X      Throw(IndexException(m,n,*this));
  530. X   return store[n-1];
  531. X}
  532. X
  533. XReal& DiagonalMatrix::operator()(int m)
  534. X{
  535. X   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  536. X   return store[m-1];
  537. X}
  538. X
  539. XReal& ColumnVector::operator()(int m)
  540. X{
  541. X   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  542. X   return store[m-1];
  543. X}
  544. X
  545. XReal& RowVector::operator()(int n)
  546. X{
  547. X   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  548. X   return store[n-1];
  549. X}
  550. X
  551. XReal& BandMatrix::operator()(int m, int n)
  552. X{
  553. X   int w = upper+lower+1; int i = lower+n-m;
  554. X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  555. X      Throw(IndexException(m,n,*this));
  556. X   return store[w*(m-1)+i];
  557. X}
  558. X
  559. XReal& UpperBandMatrix::operator()(int m, int n)
  560. X{
  561. X   int w = upper+1; int i = n-m;
  562. X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  563. X      Throw(IndexException(m,n,*this));
  564. X   return store[w*(m-1)+i];
  565. X}
  566. X
  567. XReal& LowerBandMatrix::operator()(int m, int n)
  568. X{
  569. X   int w = lower+1; int i = lower+n-m;
  570. X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  571. X      Throw(IndexException(m,n,*this));
  572. X   return store[w*(m-1)+i];
  573. X}
  574. X
  575. XReal& SymmetricBandMatrix::operator()(int m, int n)
  576. X{
  577. X   int w = lower+1;
  578. X   if (m>=n)
  579. X   {
  580. X      int i = lower+n-m;
  581. X      if ( m>nrows || n<=0 || i<0 )
  582. X         Throw(IndexException(m,n,*this));
  583. X      return store[w*(m-1)+i];
  584. X   }
  585. X   else
  586. X   {
  587. X      int i = lower+m-n;
  588. X      if ( n>nrows || m<=0 || i<0 )
  589. X         Throw(IndexException(m,n,*this));
  590. X      return store[w*(n-1)+i];
  591. X   }
  592. X}
  593. X
  594. X#ifndef __ZTC__
  595. X
  596. XReal Matrix::operator()(int m, int n) const
  597. X{
  598. X   if (m<=0 || m>nrows || n<=0 || n>ncols)
  599. X      Throw(IndexException(m,n,*this));
  600. X   return store[(m-1)*ncols+n-1];
  601. X}
  602. X
  603. XReal SymmetricMatrix::operator()(int m, int n) const
  604. X{
  605. X   if (m<=0 || n<=0 || m>nrows || n>ncols)
  606. X      Throw(IndexException(m,n,*this));
  607. X   if (m>=n) return store[tristore(m-1)+n-1];
  608. X   else return store[tristore(n-1)+m-1];
  609. X}
  610. X
  611. XReal UpperTriangularMatrix::operator()(int m, int n) const
  612. X{
  613. X   if (m<=0 || n<m || n>ncols)
  614. X      Throw(IndexException(m,n,*this));
  615. X   return store[(m-1)*ncols+n-1-tristore(m-1)];
  616. X}
  617. X
  618. XReal LowerTriangularMatrix::operator()(int m, int n) const
  619. X{
  620. X   if (n<=0 || m<n || m>nrows)
  621. X      Throw(IndexException(m,n,*this));
  622. X   return store[tristore(m-1)+n-1];
  623. X}
  624. X
  625. XReal DiagonalMatrix::operator()(int m, int n) const
  626. X{
  627. X   if (n<=0 || m!=n || m>nrows || n>ncols)
  628. X      Throw(IndexException(m,n,*this));
  629. X   return store[n-1];
  630. X}
  631. X
  632. XReal DiagonalMatrix::operator()(int m) const
  633. X{
  634. X   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
  635. X   return store[m-1];
  636. X}
  637. X
  638. XReal ColumnVector::operator()(int m) const
  639. X{
  640. X   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
  641. X   return store[m-1];
  642. X}
  643. X
  644. XReal RowVector::operator()(int n) const
  645. X{
  646. X   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
  647. X   return store[n-1];
  648. X}
  649. X
  650. XReal BandMatrix::operator()(int m, int n) const
  651. X{
  652. X   int w = upper+lower+1; int i = lower+n-m;
  653. X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
  654. X      Throw(IndexException(m,n,*this));
  655. X   return store[w*(m-1)+i];
  656. X}
  657. X
  658. XReal SymmetricBandMatrix::operator()(int m, int n) const
  659. X{
  660. X   int w = lower+1;
  661. X   if (m>=n)
  662. X   {
  663. X      int i = lower+n-m;
  664. X      if ( m>nrows || n<=0 || i<0 )
  665. X         Throw(IndexException(m,n,*this));
  666. X      return store[w*(m-1)+i];
  667. X   }
  668. X   else
  669. X   {
  670. X      int i = lower+m-n;
  671. X      if ( n>nrows || m<=0 || i<0 )
  672. X         Throw(IndexException(m,n,*this));
  673. X      return store[w*(n-1)+i];
  674. X   }
  675. X}
  676. X
  677. X#endif
  678. X
  679. XReal BaseMatrix::AsScalar() const
  680. X{
  681. X   REPORT
  682. X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  683. X   if (gm->nrows!=1 || gm->ncols!=1)
  684. X   {
  685. X      Try
  686. X      {
  687. X         Tracer tr("AsScalar");
  688. X     Throw(ProgramException("Cannot convert to scalar", *gm));
  689. X      }
  690. X      CatchAll { gm->tDelete(); Throw(); }
  691. X   }
  692. X   Real x = *(gm->store); gm->tDelete(); return x;
  693. X}
  694. X
  695. X#ifdef TEMPS_DESTROYED_QUICKLY
  696. X
  697. XAddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const
  698. X{
  699. X   REPORT
  700. X   AddedMatrix* x = new AddedMatrix(this, &bm);
  701. X   MatrixErrorNoSpace(x);
  702. X   return *x;
  703. X}
  704. X
  705. XMultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const
  706. X{
  707. X   REPORT
  708. X   MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
  709. X   MatrixErrorNoSpace(x);
  710. X   return *x;
  711. X}
  712. X
  713. X//SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
  714. XSolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
  715. X{
  716. X   REPORT
  717. X   SolvedMatrix* x;
  718. X   Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
  719. X   CatchAll { delete this; Throw(); }
  720. X   delete this;                // since we are using bm rather than this
  721. X   return *x;
  722. X}
  723. X
  724. XSubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const
  725. X{
  726. X   REPORT
  727. X   SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
  728. X   MatrixErrorNoSpace(x);
  729. X   return *x;
  730. X}
  731. X
  732. XShiftedMatrix& BaseMatrix::operator+(Real f) const
  733. X{
  734. X   REPORT
  735. X   ShiftedMatrix* x = new ShiftedMatrix(this, f);
  736. X   MatrixErrorNoSpace(x);
  737. X   return *x;
  738. X}
  739. X
  740. XScaledMatrix& BaseMatrix::operator*(Real f) const
  741. X{
  742. X   REPORT
  743. X   ScaledMatrix* x = new ScaledMatrix(this, f);
  744. X   MatrixErrorNoSpace(x);
  745. X   return *x;
  746. X}
  747. X
  748. XScaledMatrix& BaseMatrix::operator/(Real f) const
  749. X{
  750. X   REPORT
  751. X   ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
  752. X   MatrixErrorNoSpace(x);
  753. X   return *x;
  754. X}
  755. X
  756. XShiftedMatrix& BaseMatrix::operator-(Real f) const
  757. X{
  758. X   REPORT
  759. X   ShiftedMatrix* x = new ShiftedMatrix(this, -f);
  760. X   MatrixErrorNoSpace(x);
  761. X   return *x;
  762. X}
  763. X
  764. XTransposedMatrix& BaseMatrix::t() const
  765. X{
  766. X   REPORT
  767. X   TransposedMatrix* x = new TransposedMatrix(this);
  768. X   MatrixErrorNoSpace(x);
  769. X   return *x;
  770. X}
  771. X
  772. XNegatedMatrix& BaseMatrix::operator-() const
  773. X{
  774. X   REPORT
  775. X   NegatedMatrix* x = new NegatedMatrix(this);
  776. X   MatrixErrorNoSpace(x);
  777. X   return *x;
  778. X}
  779. X
  780. XInvertedMatrix& BaseMatrix::i() const
  781. X{
  782. X   REPORT
  783. X   InvertedMatrix* x = new InvertedMatrix(this);
  784. X   MatrixErrorNoSpace(x);
  785. X   return *x;
  786. X}
  787. X
  788. XConstMatrix& GeneralMatrix::c() const
  789. X{
  790. X   if (tag != -1)
  791. X      Throw(ProgramException(".c() applied to temporary matrix"));
  792. X   REPORT
  793. X   ConstMatrix* x = new ConstMatrix(this);
  794. X   MatrixErrorNoSpace(x);
  795. X   return *x;
  796. X}
  797. X
  798. XRowedMatrix& BaseMatrix::AsRow() const
  799. X{
  800. X   REPORT
  801. X   RowedMatrix* x = new RowedMatrix(this);
  802. X   MatrixErrorNoSpace(x);
  803. X   return *x;
  804. X}
  805. X
  806. XColedMatrix& BaseMatrix::AsColumn() const
  807. X{
  808. X   REPORT
  809. X   ColedMatrix* x = new ColedMatrix(this);
  810. X   MatrixErrorNoSpace(x);
  811. X   return *x;
  812. X}
  813. X
  814. XDiagedMatrix& BaseMatrix::AsDiagonal() const
  815. X{
  816. X   REPORT
  817. X   DiagedMatrix* x = new DiagedMatrix(this);
  818. X   MatrixErrorNoSpace(x);
  819. X   return *x;
  820. X}
  821. X
  822. XMatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const
  823. X{
  824. X   REPORT
  825. X   MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
  826. X   MatrixErrorNoSpace(x);
  827. X   return *x;
  828. X}
  829. X
  830. X#else
  831. X
  832. XAddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
  833. X{ REPORT return AddedMatrix(this, &bm); }
  834. X
  835. XMultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
  836. X{ REPORT return MultipliedMatrix(this, &bm); }
  837. X
  838. XSolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
  839. X{ REPORT return SolvedMatrix(bm, &bmx); }
  840. X
  841. XSubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
  842. X{ REPORT return SubtractedMatrix(this, &bm); }
  843. X
  844. XShiftedMatrix BaseMatrix::operator+(Real f) const
  845. X{ REPORT return ShiftedMatrix(this, f); }
  846. X
  847. XScaledMatrix BaseMatrix::operator*(Real f) const
  848. X{ REPORT return ScaledMatrix(this, f); }
  849. X
  850. XScaledMatrix BaseMatrix::operator/(Real f) const
  851. X{ REPORT return ScaledMatrix(this, 1.0/f); }
  852. X
  853. XShiftedMatrix BaseMatrix::operator-(Real f) const
  854. X{ REPORT return ShiftedMatrix(this, -f); }
  855. X
  856. XTransposedMatrix BaseMatrix::t() const
  857. X{ REPORT return TransposedMatrix(this); }
  858. X
  859. XNegatedMatrix BaseMatrix::operator-() const
  860. X{ REPORT return NegatedMatrix(this); }
  861. X
  862. XInvertedMatrix BaseMatrix::i() const
  863. X{ REPORT return InvertedMatrix(this); }
  864. X
  865. XConstMatrix GeneralMatrix::c() const
  866. X{
  867. X   if (tag != -1)
  868. X      Throw(ProgramException(".c() applied to temporary matrix"));
  869. X   REPORT return ConstMatrix(this);
  870. X}
  871. X
  872. XRowedMatrix BaseMatrix::AsRow() const
  873. X{ REPORT return RowedMatrix(this); }
  874. X
  875. XColedMatrix BaseMatrix::AsColumn() const
  876. X{ REPORT return ColedMatrix(this); }
  877. X
  878. XDiagedMatrix BaseMatrix::AsDiagonal() const
  879. X{ REPORT return DiagedMatrix(this); }
  880. X
  881. XMatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
  882. X{ REPORT return MatedMatrix(this,nrx,ncx); }
  883. X
  884. X#endif
  885. X
  886. Xvoid GeneralMatrix::operator=(Real f)
  887. X{ REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
  888. X
  889. Xvoid Matrix::operator=(const BaseMatrix& X)
  890. X{
  891. X   REPORT //CheckConversion(X);
  892. X   MatrixConversionCheck mcc;
  893. X   Eq(X,MatrixType::Rt);
  894. X} 
  895. X
  896. Xvoid RowVector::operator=(const BaseMatrix& X)
  897. X{
  898. X   REPORT  // CheckConversion(X);
  899. X   MatrixConversionCheck mcc;
  900. X   Eq(X,MatrixType::RV);
  901. X   if (nrows!=1)
  902. X   {
  903. X      Tracer tr("RowVector(=)");
  904. X      Throw(VectorException(*this));
  905. X   }
  906. X}
  907. X
  908. Xvoid ColumnVector::operator=(const BaseMatrix& X)
  909. X{
  910. X   REPORT //CheckConversion(X);
  911. X   MatrixConversionCheck mcc;
  912. X   Eq(X,MatrixType::CV);
  913. X   if (ncols!=1)
  914. X   {
  915. X      Tracer tr("ColumnVector(=)");
  916. X      Throw(VectorException(*this));
  917. X   }
  918. X}
  919. X
  920. Xvoid SymmetricMatrix::operator=(const BaseMatrix& X)
  921. X{
  922. X   REPORT // CheckConversion(X);
  923. X   MatrixConversionCheck mcc;
  924. X   Eq(X,MatrixType::Sm);
  925. X}
  926. Xvoid UpperTriangularMatrix::operator=(const BaseMatrix& X)
  927. X{
  928. X   REPORT //CheckConversion(X);
  929. X   MatrixConversionCheck mcc;
  930. X   Eq(X,MatrixType::UT);
  931. X}
  932. X
  933. Xvoid LowerTriangularMatrix::operator=(const BaseMatrix& X)
  934. X{
  935. X   REPORT //CheckConversion(X);
  936. X   MatrixConversionCheck mcc;
  937. X   Eq(X,MatrixType::LT);
  938. X}
  939. X
  940. Xvoid DiagonalMatrix::operator=(const BaseMatrix& X)
  941. X{
  942. X   REPORT // CheckConversion(X);
  943. X   MatrixConversionCheck mcc;
  944. X   Eq(X,MatrixType::Dg);
  945. X}
  946. X
  947. Xvoid GeneralMatrix::operator<<(const Real* r)
  948. X{
  949. X   REPORT
  950. X   int i = storage; Real* s=store;
  951. X   while(i--) *s++ = *r++;
  952. X}
  953. X
  954. X
  955. X/************************* element access *********************************/
  956. X
  957. XReal& Matrix::element(int m, int n)
  958. X{
  959. X   if (m<0 || m>= nrows || n<0 || n>= ncols)
  960. X      Throw(IndexException(m,n,*this,TRUE));
  961. X   return store[m*ncols+n];
  962. X}
  963. X
  964. XReal& SymmetricMatrix::element(int m, int n)
  965. X{
  966. X   if (m<0 || n<0 || m >= nrows || n>=ncols)
  967. X      Throw(IndexException(m,n,*this,TRUE));
  968. X   if (m>=n) return store[tristore(m)+n];
  969. X   else return store[tristore(n)+m];
  970. X}
  971. X
  972. XReal& UpperTriangularMatrix::element(int m, int n)
  973. X{
  974. X   if (m<0 || n<m || n>=ncols)
  975. X      Throw(IndexException(m,n,*this,TRUE));
  976. X   return store[m*ncols+n-tristore(m)];
  977. X}
  978. X
  979. XReal& LowerTriangularMatrix::element(int m, int n)
  980. X{
  981. X   if (n<0 || m<n || m>=nrows)
  982. X      Throw(IndexException(m,n,*this,TRUE));
  983. X   return store[tristore(m)+n];
  984. X}
  985. X
  986. XReal& DiagonalMatrix::element(int m, int n)
  987. X{
  988. X   if (n<0 || m!=n || m>=nrows || n>=ncols)
  989. X      Throw(IndexException(m,n,*this,TRUE));
  990. X   return store[n];
  991. X}
  992. X
  993. XReal& DiagonalMatrix::element(int m)
  994. X{
  995. X   if (m<0 || m>=nrows) Throw(IndexException(m,*this,TRUE));
  996. X   return store[m];
  997. X}
  998. X
  999. XReal& ColumnVector::element(int m)
  1000. X{
  1001. X   if (m<0 || m>= nrows) Throw(IndexException(m,*this,TRUE));
  1002. X   return store[m];
  1003. X}
  1004. X
  1005. XReal& RowVector::element(int n)
  1006. X{
  1007. X   if (n<0 || n>= ncols)  Throw(IndexException(n,*this,TRUE));
  1008. X   return store[n];
  1009. X}
  1010. X
  1011. XReal& BandMatrix::element(int m, int n)
  1012. X{
  1013. X   int w = upper+lower+1; int i = lower+n-m;
  1014. X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  1015. X      Throw(IndexException(m,n,*this,TRUE));
  1016. X   return store[w*m+i];
  1017. X}
  1018. X
  1019. XReal& UpperBandMatrix::element(int m, int n)
  1020. X{
  1021. X   int w = upper+1; int i = n-m;
  1022. X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  1023. X      Throw(IndexException(m,n,*this,TRUE));
  1024. X   return store[w*m+i];
  1025. X}
  1026. X
  1027. XReal& LowerBandMatrix::element(int m, int n)
  1028. X{
  1029. X   int w = lower+1; int i = lower+n-m;
  1030. X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
  1031. X      Throw(IndexException(m,n,*this,TRUE));
  1032. X   return store[w*m+i];
  1033. X}
  1034. X
  1035. XReal& SymmetricBandMatrix::element(int m, int n)
  1036. X{
  1037. X   int w = lower+1;
  1038. X   if (m>=n)
  1039. X   {
  1040. X      int i = lower+n-m;
  1041. X      if ( m>=nrows || n<0 || i<0 )
  1042. X         Throw(IndexException(m,n,*this,TRUE));
  1043. X      return store[w*m+i];
  1044. X   }
  1045. X   else
  1046. X   {
  1047. X      int i = lower+m-n;
  1048. X      if ( n>=nrows || m<0 || i<0 )
  1049. X         Throw(IndexException(m,n,*this,TRUE));
  1050. X      return store[w*n+i];
  1051. X   }
  1052. X}
  1053. X
  1054. END_OF_FILE
  1055. if test 13261 -ne `wc -c <'newmat6.cxx'`; then
  1056.     echo shar: \"'newmat6.cxx'\" unpacked with wrong size!
  1057. fi
  1058. # end of 'newmat6.cxx'
  1059. fi
  1060. if test -f 'newmat7.cxx' -a "${1}" != "-c" ; then 
  1061.   echo shar: Will not clobber existing file \"'newmat7.cxx'\"
  1062. else
  1063. echo shar: Extracting \"'newmat7.cxx'\" \(16089 characters\)
  1064. sed "s/^X//" >'newmat7.cxx' <<'END_OF_FILE'
  1065. X//$$ newmat7.cxx     Invert, solve, binary operations
  1066. X
  1067. X// Copyright (C) 1991,2,3: R B Davies
  1068. X
  1069. X#include "include.h"
  1070. X
  1071. X#include "newmat.h"
  1072. X#include "newmatrc.h"
  1073. X
  1074. X//#define REPORT { static ExeCounter ExeCount(__LINE__,7); ExeCount++; }
  1075. X
  1076. X#define REPORT {}
  1077. X
  1078. X
  1079. X/***************************** solve routines ******************************/
  1080. X
  1081. XGeneralMatrix* GeneralMatrix::MakeSolver()
  1082. X{
  1083. X   REPORT
  1084. X   GeneralMatrix* gm = new CroutMatrix(*this);
  1085. X   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  1086. X}
  1087. X
  1088. XGeneralMatrix* Matrix::MakeSolver()
  1089. X{
  1090. X   REPORT
  1091. X   GeneralMatrix* gm = new CroutMatrix(*this);
  1092. X   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  1093. X}
  1094. X
  1095. Xvoid CroutMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  1096. X{
  1097. X   REPORT
  1098. X   Real* el = mcin.store; int i = mcin.skip;
  1099. X   while (i--) *el++ = 0.0;
  1100. X   el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  1101. X   while (i--) *el++ = 0.0;
  1102. X   lubksb(mcin.store, mcout.skip);
  1103. X}
  1104. X
  1105. X
  1106. X// Do we need check for entirely zero output?
  1107. X
  1108. Xvoid UpperTriangularMatrix::Solver(MatrixRowCol& mcout,
  1109. X   const MatrixRowCol& mcin)
  1110. X{
  1111. X   REPORT
  1112. X   Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  1113. X   while (i-- > 0) *elx++ = 0.0;
  1114. X   int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; Real* el = elx;
  1115. X   int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
  1116. X   while (j-- > 0) *elx++ = 0.0;
  1117. X   Real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
  1118. X   while (i-- > 0)
  1119. X   {
  1120. X      elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
  1121. X      while (jx--) sum += *(--Ael) * *(--elx);
  1122. X      elx--; *elx = (*elx - sum) / *(--Ael);
  1123. X   }
  1124. X}
  1125. X
  1126. Xvoid LowerTriangularMatrix::Solver(MatrixRowCol& mcout,
  1127. X   const MatrixRowCol& mcin)
  1128. X{
  1129. X   REPORT
  1130. X   Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  1131. X   while (i-- > 0) *elx++ = 0.0;
  1132. X   int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  1133. X   int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  1134. X   while (j-- > 0) *elx++ = 0.0;
  1135. X   Real* el = mcin.store+nc; Real* Ael = store + (nc*(nc+1))/2; j = 0;
  1136. X   while (i-- > 0)
  1137. X   {
  1138. X      elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
  1139. X      while (jx--) sum += *Ael++ * *elx++;
  1140. X      *elx = (*elx - sum) / *Ael++;
  1141. X   }
  1142. X}
  1143. X
  1144. X/******************* carry out binary operations *************************/
  1145. X
  1146. Xstatic GeneralMatrix*
  1147. X   GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
  1148. Xstatic GeneralMatrix*
  1149. X   GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
  1150. Xstatic GeneralMatrix*
  1151. X   GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
  1152. Xstatic GeneralMatrix*
  1153. X   GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
  1154. X
  1155. XGeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
  1156. X{
  1157. X   REPORT
  1158. X   gm1=((BaseMatrix*&)bm1)->Evaluate();
  1159. X   gm2=((BaseMatrix*&)bm2)->Evaluate();
  1160. X#ifdef TEMPS_DESTROYED_QUICKLY
  1161. X   GeneralMatrix* gmx;
  1162. X   Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
  1163. X   CatchAll { delete this; Throw(); }
  1164. X   delete this; return gmx;
  1165. X#else
  1166. X   return GeneralAdd(gm1,gm2,this,mt);
  1167. X#endif   
  1168. X}
  1169. X
  1170. XGeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  1171. X{
  1172. X   REPORT
  1173. X   gm1=((BaseMatrix*&)bm1)->Evaluate();
  1174. X   gm2=((BaseMatrix*&)bm2)->Evaluate();
  1175. X#ifdef TEMPS_DESTROYED_QUICKLY
  1176. X   GeneralMatrix* gmx;
  1177. X   Try { gmx = GeneralSub(gm1,gm2,this,mt); }
  1178. X   CatchAll { delete this; Throw(); }
  1179. X   delete this; return gmx;
  1180. X#else
  1181. X   return GeneralSub(gm1,gm2,this,mt);
  1182. X#endif   
  1183. X}
  1184. X
  1185. XGeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  1186. X{
  1187. X   REPORT
  1188. X   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  1189. X   gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
  1190. X   gm1=((BaseMatrix*&)bm1)->Evaluate();
  1191. X#ifdef TEMPS_DESTROYED_QUICKLY
  1192. X   GeneralMatrix* gmx;
  1193. X   Try { gmx = GeneralMult(gm1, gm2, this, mt); }
  1194. X   CatchAll { delete this; Throw(); }
  1195. X   delete this; return gmx;
  1196. X#else
  1197. X   return GeneralMult(gm1, gm2, this, mt);
  1198. X#endif   
  1199. X}
  1200. X
  1201. XGeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  1202. X{
  1203. X   REPORT
  1204. X   gm1=((BaseMatrix*&)bm1)->Evaluate();
  1205. X   gm2=((BaseMatrix*&)bm2)->Evaluate();
  1206. X#ifdef TEMPS_DESTROYED_QUICKLY
  1207. X   GeneralMatrix* gmx;
  1208. X   Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
  1209. X   CatchAll { delete this; Throw(); }
  1210. X   delete this; return gmx;
  1211. X#else
  1212. X   return GeneralSolv(gm1,gm2,this,mt);
  1213. X#endif   
  1214. X}
  1215. X
  1216. X// routines for adding or subtracting matrices of identical storage structure
  1217. X
  1218. Xstatic void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  1219. X{
  1220. X   REPORT
  1221. X   Real* s1=gm1->Store(); Real* s2=gm2->Store();
  1222. X   Real* s=gm->Store(); int i=gm->Storage() >> 2;
  1223. X   while (i--)
  1224. X   {
  1225. X       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  1226. X       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  1227. X   }
  1228. X   i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
  1229. X}
  1230. X   
  1231. Xstatic void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  1232. X{
  1233. X   REPORT
  1234. X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  1235. X   while (i--)
  1236. X   { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
  1237. X   i=gm->Storage() & 3; while (i--) *s++ += *s2++;
  1238. X}
  1239. X
  1240. Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  1241. X{
  1242. X   REPORT
  1243. X   Real* s1=gm1->Store(); Real* s2=gm2->Store();
  1244. X   Real* s=gm->Store(); int i=gm->Storage() >> 2;
  1245. X   while (i--)
  1246. X   {
  1247. X       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  1248. X       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  1249. X   }
  1250. X   i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
  1251. X}
  1252. X
  1253. Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  1254. X{
  1255. X   REPORT
  1256. X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  1257. X   while (i--)
  1258. X   { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
  1259. X   i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
  1260. X}
  1261. X
  1262. Xstatic void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  1263. X{
  1264. X   REPORT
  1265. X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  1266. X   while (i--)
  1267. X   {
  1268. X      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  1269. X      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  1270. X   }
  1271. X   i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
  1272. X}
  1273. X
  1274. X// routines for adding or subtracting matrices of different storage structure
  1275. X
  1276. Xstatic void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  1277. X{
  1278. X   int nr = gm->Nrows();
  1279. X   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  1280. X   MatrixRow mr(gm, StoreOnExit+DirectPart);
  1281. X   while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  1282. X}
  1283. X   
  1284. Xstatic void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  1285. X// Add into first argument
  1286. X{
  1287. X   int nr = gm->Nrows();
  1288. X   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
  1289. X   MatrixRow mr2(gm2, LoadOnEntry);
  1290. X   while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
  1291. X}
  1292. X
  1293. Xstatic void SubtractDS
  1294. X   (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  1295. X{
  1296. X   int nr = gm->Nrows();
  1297. X   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  1298. X   MatrixRow mr(gm, StoreOnExit+DirectPart);
  1299. X   while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  1300. X}
  1301. X
  1302. Xstatic void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  1303. X{
  1304. X   int nr = gm->Nrows();
  1305. X   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  1306. X   MatrixRow mr2(gm2, LoadOnEntry);
  1307. X   while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
  1308. X}
  1309. X
  1310. Xstatic void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  1311. X{
  1312. X   int nr = gm->Nrows();
  1313. X   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  1314. X   MatrixRow mr2(gm2, LoadOnEntry);
  1315. X   while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
  1316. X}
  1317. X
  1318. X#ifdef __GNUG__
  1319. Xvoid AddedMatrix::SelectVersion
  1320. X   (MatrixType mtx, int& c1, int& c2) const
  1321. X#else
  1322. Xvoid AddedMatrix::SelectVersion
  1323. X   (MatrixType mtx, Boolean& c1, Boolean& c2) const
  1324. X#endif
  1325. X// for determining version of add and subtract routines
  1326. X// will need to modify if further matrix structures are introduced
  1327. X{
  1328. X   MatrixBandWidth bm1 = gm1->BandWidth();
  1329. X   MatrixBandWidth bm2 = gm2->BandWidth();
  1330. X   MatrixBandWidth bmx = bm1 + bm2;
  1331. X   c1 = (mtx == gm1->Type()) && (bmx == bm1);
  1332. X   c2 = (mtx == gm2->Type()) && (bmx == bm2);
  1333. X}
  1334. X
  1335. Xstatic GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1336. X   AddedMatrix* am, MatrixType mtx)
  1337. X{
  1338. X   Tracer tr("GeneralAdd");
  1339. X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
  1340. X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
  1341. X      Throw(IncompatibleDimensionsException());
  1342. X   Compare(gm1->Type() + gm2->Type(),mtx);
  1343. X#ifdef __GNUG__
  1344. X   int c1,c2; am->SelectVersion(mtx,c1,c2);
  1345. X#else
  1346. X   Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
  1347. X#endif
  1348. X   if (c1 && c2)
  1349. X   {
  1350. X      if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  1351. X      else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  1352. X      else
  1353. X      {
  1354. X         REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
  1355. X         gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  1356. X      }
  1357. X   }
  1358. X   else
  1359. X   {
  1360. X      if (c1 && gm1->reuse() )               // must have type test first
  1361. X      { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
  1362. X      else if (c2 && gm2->reuse() )
  1363. X      { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
  1364. X      else
  1365. X      {
  1366. X         REPORT
  1367. X     GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
  1368. X     if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  1369. X         gmx->ReleaseAndDelete(); return gmx;
  1370. X      }
  1371. X   }
  1372. X}
  1373. X
  1374. X
  1375. Xstatic GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1376. X   SubtractedMatrix* sm, MatrixType mtx)
  1377. X{
  1378. X   Tracer tr("GeneralSub");
  1379. X   Compare(gm1->Type() + gm2->Type(),mtx);
  1380. X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
  1381. X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  1382. X      Throw(IncompatibleDimensionsException());
  1383. X#ifdef __GNUG__
  1384. X   int c1,c2; sm->SelectVersion(mtx,c1,c2);
  1385. X#else
  1386. X   Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
  1387. X#endif
  1388. X   if (c1 && c2)
  1389. X   {
  1390. X      if (gm1->reuse())
  1391. X      { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  1392. X      else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  1393. X      else
  1394. X      {
  1395. X         REPORT
  1396. X     GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
  1397. X         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  1398. X      }
  1399. X   }
  1400. X   else
  1401. X   {
  1402. X      if ( c1 && gm1->reuse() )
  1403. X      { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
  1404. X      else if ( c2 && gm2->reuse() )
  1405. X      {
  1406. X         REPORT
  1407. X         ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
  1408. X      }
  1409. X      else
  1410. X      {
  1411. X         REPORT
  1412. X     GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
  1413. X     if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  1414. X     gmx->ReleaseAndDelete(); return gmx;
  1415. X      }
  1416. X   }
  1417. X}
  1418. X
  1419. Xstatic GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1420. X   MultipliedMatrix* mm, MatrixType mtx)
  1421. X{
  1422. X   REPORT
  1423. X   Tracer tr("GeneralMult1");
  1424. X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
  1425. X   if (gm1->Ncols() !=gm2->Nrows())
  1426. X      Throw(IncompatibleDimensionsException());
  1427. X   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  1428. X
  1429. X   MatrixCol mcx(gmx, StoreOnExit+DirectPart);
  1430. X   MatrixCol mc2(gm2, LoadOnEntry);
  1431. X   while (nc--)
  1432. X   {
  1433. X      MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
  1434. X      Real* el = mcx();                              // pointer to an element
  1435. X      int n = mcx.Storage();
  1436. X      while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
  1437. X      mc2.Next(); mcx.Next();
  1438. X   }
  1439. X   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  1440. X}
  1441. X
  1442. Xstatic GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1443. X   MultipliedMatrix* mm, MatrixType mtx)
  1444. X{
  1445. X   // version that accesses by row only - not good for thin matrices
  1446. X   // or column vectors in right hand term. Needs fixing
  1447. X   REPORT
  1448. X   Tracer tr("GeneralMult2");
  1449. X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
  1450. X   if (gm1->Ncols() !=gm2->Nrows())
  1451. X      Throw(IncompatibleDimensionsException());
  1452. X   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  1453. X
  1454. X   Real* el = gmx->Store(); int n = gmx->Storage();
  1455. X   while (n--) *el++ = 0.0;
  1456. X   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
  1457. X   MatrixRow mr1(gm1, LoadOnEntry);
  1458. X   while (nr--)
  1459. X   {
  1460. X      MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
  1461. X      el = mr1();                              // pointer to an element
  1462. X      n = mr1.Storage();
  1463. X      while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
  1464. X      mr1.Next(); mrx.Next();
  1465. X   }
  1466. X   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  1467. X}
  1468. X
  1469. Xstatic GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
  1470. X{
  1471. X   // matrix multiplication for type Matrix only
  1472. X   REPORT
  1473. X   Tracer tr("MatrixMult");
  1474. X
  1475. X   int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
  1476. X   if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
  1477. X
  1478. X   Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
  1479. X
  1480. X   Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
  1481. X   
  1482. X   if (ncr)
  1483. X   {
  1484. X      while (nr--)
  1485. X      {
  1486. X         Real* s2x = s2; int j = ncr;
  1487. X         Real* sx = s; Real f = *s1++; int k = nc;
  1488. X         while (k--) *sx++ = f * *s2x++;
  1489. X         while (--j)
  1490. X            { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
  1491. X         s = sx;
  1492. X      }
  1493. X   }
  1494. X   else *gm = 0.0;
  1495. X
  1496. X   gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
  1497. X}
  1498. X
  1499. Xstatic GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1500. X   MultipliedMatrix* mm, MatrixType mtx)
  1501. X{
  1502. X   if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
  1503. X   else
  1504. X   {
  1505. X      Compare(gm1->Type() * gm2->Type(),mtx);
  1506. X      int nr = gm2->Nrows(); int nc = gm2->Ncols();
  1507. X      if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
  1508. X      else return GeneralMult2(gm1, gm2, mm, mtx);
  1509. X   }
  1510. X}
  1511. X
  1512. Xstatic GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
  1513. X   BaseMatrix* sm, MatrixType mtx)
  1514. X{
  1515. X   REPORT
  1516. X   Tracer tr("GeneralSolv");
  1517. X   Compare(gm1->Type().i() * gm2->Type(),mtx);
  1518. X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
  1519. X   if (gm1->Ncols() !=gm2->Nrows())
  1520. X      Throw(IncompatibleDimensionsException());
  1521. X   GeneralMatrix* gmx = mtx.New(nr,nc,sm);
  1522. X
  1523. X   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
  1524. X#ifndef ATandT
  1525. X   MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
  1526. X                           // deleted for ATandT, to balance deletion below
  1527. X#endif
  1528. X   MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
  1529. X   MatrixCol mc2(gm2, r, LoadOnEntry);
  1530. X   GeneralMatrix* gms = gm1->MakeSolver();
  1531. X   Try
  1532. X   {
  1533. X      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
  1534. X   }
  1535. X   CatchAll
  1536. X   {
  1537. X      gms->tDelete(); delete gmx; gm2->tDelete();
  1538. X#ifndef ATandT
  1539. X      MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  1540. X                          // ATandT version 2.1 gives an internal error
  1541. X#endif
  1542. X#ifdef Version21
  1543. X      delete [] r;
  1544. X#else
  1545. X      delete [nr] r;
  1546. X#endif
  1547. X      Throw();
  1548. X   }
  1549. X   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
  1550. X#ifndef ATandT
  1551. X   MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  1552. X                          // ATandT version 2.1 gives an internal error
  1553. X#endif
  1554. X#ifdef Version21
  1555. X   delete [] r;
  1556. X#else
  1557. X   delete [nr] r;
  1558. X#endif
  1559. X   return gmx;
  1560. X}
  1561. X
  1562. XGeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
  1563. X{
  1564. X   // Matrix Inversion - use solve routines
  1565. X   REPORT
  1566. X   gm=((BaseMatrix*&)bm)->Evaluate();
  1567. X   int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
  1568. X#ifdef TEMPS_DESTROYED_QUICKLY
  1569. X   GeneralMatrix* gmx;
  1570. X   Try { gmx = GeneralSolv(gm,&I,this,mtx); }
  1571. X   CatchAll { delete this; Throw(); }
  1572. X   delete this; return gmx;
  1573. X#else
  1574. X   return GeneralSolv(gm,&I,this,mtx);
  1575. X#endif   
  1576. X}
  1577. X
  1578. X
  1579. X/*************************** norm functions ****************************/
  1580. X
  1581. XReal BaseMatrix::Norm1() const
  1582. X{
  1583. X   // maximum of sum of absolute values of a column
  1584. X   REPORT
  1585. X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1586. X   int nc = gm->Ncols(); Real value = 0.0;
  1587. X   MatrixCol mc(gm, LoadOnEntry);
  1588. X   while (nc--)
  1589. X      { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
  1590. X   gm->tDelete(); return value;
  1591. X}
  1592. X
  1593. XReal BaseMatrix::NormInfinity() const
  1594. X{
  1595. X   // maximum of sum of absolute values of a row
  1596. X   REPORT
  1597. X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1598. X   int nr = gm->Nrows(); Real value = 0.0;
  1599. X   MatrixRow mr(gm, LoadOnEntry);
  1600. X   while (nr--)
  1601. X      { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
  1602. X   gm->tDelete(); return value;
  1603. X}
  1604. END_OF_FILE
  1605. if test 16089 -ne `wc -c <'newmat7.cxx'`; then
  1606.     echo shar: \"'newmat7.cxx'\" unpacked with wrong size!
  1607. fi
  1608. # end of 'newmat7.cxx'
  1609. fi
  1610. if test -f 'newmat8.cxx' -a "${1}" != "-c" ; then 
  1611.   echo shar: Will not clobber existing file \"'newmat8.cxx'\"
  1612. else
  1613. echo shar: Extracting \"'newmat8.cxx'\" \(8288 characters\)
  1614. sed "s/^X//" >'newmat8.cxx' <<'END_OF_FILE'
  1615. X//$$ newmat8.cxx         Advanced LU transform, scalar functions
  1616. X
  1617. X// Copyright (C) 1991,2,3: R B Davies
  1618. X
  1619. X#define WANT_MATH
  1620. X
  1621. X#include "include.h"
  1622. X
  1623. X#include "newmatap.h"
  1624. X
  1625. X//#define REPORT { static ExeCounter ExeCount(__LINE__,8); ExeCount++; }
  1626. X
  1627. X#define REPORT {}
  1628. X
  1629. X
  1630. X/************************** LU transformation ****************************/
  1631. X
  1632. Xvoid CroutMatrix::ludcmp()
  1633. X// LU decomposition - from numerical recipes in C
  1634. X{
  1635. X   REPORT
  1636. X   Tracer trace("Crout(ludcmp)");
  1637. X   int i,j;
  1638. X
  1639. X   Real* vv=new Real [nrows]; MatrixErrorNoSpace(vv);
  1640. X   MONITOR_REAL_NEW("Make  (CroutMat)",nrows,vv)
  1641. X   Real* a;
  1642. X
  1643. X   a=store;
  1644. X   for (i=0;i<nrows;i++)
  1645. X   {
  1646. X      Real big=0.0;
  1647. X      j=nrows; while (j--) { Real temp=fabs(*a++); if (temp > big) big=temp; }
  1648. X      if (big == 0.0)
  1649. X      {
  1650. X         MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  1651. X#ifdef Version21
  1652. X         sing = TRUE; delete [] vv; return;
  1653. X#else
  1654. X         sing = TRUE; delete [nrows] vv; return;
  1655. X#endif
  1656. X      }
  1657. X      vv[i]=1.0/big;
  1658. X   }
  1659. X
  1660. X   Real* aj=store;
  1661. X   for (j=0;j<nrows;j++)
  1662. X   {
  1663. X      Real* ai=store;
  1664. X      for (i=0;i<j;i++)
  1665. X      {
  1666. X         Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  1667. X         int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  1668. X         *ajx = sum; ai += nrows;
  1669. X      }
  1670. X
  1671. X      Real big=0.0; int imax;
  1672. X      for (i=j;i<nrows;i++)
  1673. X      {
  1674. X         Real sum=*(ai+j); Real* aix=ai; Real* ajx=aj;
  1675. X         int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
  1676. X         *aix = sum; ai += nrows;
  1677. X         Real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
  1678. X      }
  1679. X
  1680. X      if (j != imax)
  1681. X      {
  1682. X         Real* amax=store+imax*nrows; Real* ajx=store+j*nrows;
  1683. X         int k=nrows; while(k--) { Real dum=*amax; *amax++=*ajx; *ajx++=dum; }
  1684. X         d=!d; vv[imax]=vv[j];
  1685. X      }
  1686. X
  1687. X      indx[j]=imax; ai=aj+j*nrows;
  1688. X      if (*ai == 0.0)
  1689. X      {
  1690. X         MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  1691. X#ifdef Version21
  1692. X         sing = TRUE; delete [] vv; return;
  1693. X#else
  1694. X         sing = TRUE; delete [nrows] vv; return;
  1695. X#endif
  1696. X      }
  1697. X      Real dum=1.0/(*ai);
  1698. X      i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
  1699. X
  1700. X      aj++;
  1701. X   }
  1702. X   MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
  1703. X#ifdef Version21
  1704. X   delete [] vv;
  1705. X#else
  1706. X   delete [nrows] vv;
  1707. X#endif
  1708. X}
  1709. X
  1710. Xvoid CroutMatrix::lubksb(Real* B, int mini)
  1711. X{
  1712. X   REPORT
  1713. X   Tracer trace("Crout(lubksb)");
  1714. X   if (sing) Throw(SingularException(*this));   
  1715. X   int i,j; int ii=-1; Real* ai=store;
  1716. X
  1717. X   for (i=0;i<nrows;i++)
  1718. X   {
  1719. X      int ip=indx[i]; Real sum=B[ip]; B[ip]=B[i];
  1720. X      if (ii>=0)
  1721. X      {
  1722. X         Real* aij=ai+ii; Real* bj=B+ii; j=i-ii;
  1723. X         while (j--) { sum -= *aij++ * *bj; bj++; }
  1724. X      }
  1725. X      else if (sum) ii=i;
  1726. X      B[i]=sum; ai += nrows;
  1727. X   }
  1728. X
  1729. X   for (i=nrows-1;i>=mini;i--)
  1730. X   {
  1731. X      Real* bj=B+i; ai -= nrows; Real* ajx=ai+i; Real sum=*bj; bj++;
  1732. X      j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
  1733. X      B[i]=sum/(*(ai+i));
  1734. X   }
  1735. X}
  1736. X
  1737. X
  1738. X/****************************** scalar functions ****************************/
  1739. X
  1740. Xinline Real square(Real x) { return x*x; }
  1741. X
  1742. XReal GeneralMatrix::SumSquare() const
  1743. X{
  1744. X   REPORT
  1745. X   Real sum = 0.0; int i = storage; Real* s = store;
  1746. X   while (i--) sum += square(*s++);
  1747. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1748. X}
  1749. X
  1750. XReal GeneralMatrix::SumAbsoluteValue() const
  1751. X{
  1752. X   REPORT
  1753. X   Real sum = 0.0; int i = storage; Real* s = store;
  1754. X   while (i--) sum += fabs(*s++);
  1755. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1756. X}
  1757. X
  1758. XReal GeneralMatrix::MaximumAbsoluteValue() const
  1759. X{
  1760. X   REPORT
  1761. X   Real maxval = 0.0; int i = storage; Real* s = store;
  1762. X   while (i--) { Real a = fabs(*s++); if (maxval < a) maxval = a; }
  1763. X   ((GeneralMatrix&)*this).tDelete(); return maxval;
  1764. X}
  1765. X
  1766. XReal SymmetricMatrix::SumSquare() const
  1767. X{
  1768. X   REPORT
  1769. X   Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  1770. X   for (int i = 0; i<nr; i++)
  1771. X   {
  1772. X      int j = i;
  1773. X      while (j--) sum2 += square(*s++);
  1774. X      sum1 += square(*s++);
  1775. X   }
  1776. X   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  1777. X}
  1778. X
  1779. XReal SymmetricMatrix::SumAbsoluteValue() const
  1780. X{
  1781. X   REPORT
  1782. X   Real sum1 = 0.0; Real sum2 = 0.0; Real* s = store; int nr = nrows;
  1783. X   for (int i = 0; i<nr; i++)
  1784. X   {
  1785. X      int j = i;
  1786. X      while (j--) sum2 += fabs(*s++);
  1787. X      sum1 += fabs(*s++);
  1788. X   }
  1789. X   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
  1790. X}
  1791. X
  1792. XReal BaseMatrix::SumSquare() const
  1793. X{
  1794. X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1795. X   Real s = gm->SumSquare(); return s;
  1796. X}
  1797. X
  1798. XReal BaseMatrix::SumAbsoluteValue() const
  1799. X{
  1800. X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1801. X   Real s = gm->SumAbsoluteValue(); return s;
  1802. X}
  1803. X
  1804. XReal BaseMatrix::MaximumAbsoluteValue() const
  1805. X{
  1806. X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1807. X   Real s = gm->MaximumAbsoluteValue(); return s;
  1808. X}
  1809. X
  1810. XReal Matrix::Trace() const
  1811. X{
  1812. X   REPORT
  1813. X   Tracer trace("Trace");
  1814. X   int i = nrows; int d = i+1;
  1815. X   if (i != ncols) Throw(NotSquareException(*this));
  1816. X   Real sum = 0.0; Real* s = store;
  1817. X   while (i--) { sum += *s; s += d; }
  1818. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1819. X}
  1820. X
  1821. XReal DiagonalMatrix::Trace() const
  1822. X{
  1823. X   REPORT
  1824. X   int i = nrows; Real sum = 0.0; Real* s = store;
  1825. X   while (i--) sum += *s++;
  1826. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1827. X}
  1828. X
  1829. XReal SymmetricMatrix::Trace() const
  1830. X{
  1831. X   REPORT
  1832. X   int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  1833. X   while (i--) { sum += *s; s += j++; }
  1834. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1835. X}
  1836. X
  1837. XReal LowerTriangularMatrix::Trace() const
  1838. X{
  1839. X   REPORT
  1840. X   int i = nrows; Real sum = 0.0; Real* s = store; int j = 2;
  1841. X   while (i--) { sum += *s; s += j++; }
  1842. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1843. X}
  1844. X
  1845. XReal UpperTriangularMatrix::Trace() const
  1846. X{
  1847. X   REPORT
  1848. X   int i = nrows; Real sum = 0.0; Real* s = store;
  1849. X   while (i) { sum += *s; s += i--; }
  1850. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1851. X}
  1852. X
  1853. XReal BandMatrix::Trace() const
  1854. X{
  1855. X   REPORT
  1856. X   int i = nrows; int w = lower+upper+1;
  1857. X   Real sum = 0.0; Real* s = store+lower;
  1858. X   while (i--) { sum += *s; s += w; }
  1859. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1860. X}
  1861. X
  1862. XReal SymmetricBandMatrix::Trace() const
  1863. X{
  1864. X   REPORT
  1865. X   int i = nrows; int w = lower+1;
  1866. X   Real sum = 0.0; Real* s = store+lower;
  1867. X   while (i--) { sum += *s; s += w; }
  1868. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1869. X}
  1870. X
  1871. XReal BaseMatrix::Trace() const
  1872. X{
  1873. X   REPORT
  1874. X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(MatrixType::Dg);
  1875. X   Real sum = gm->Trace(); return sum;
  1876. X}
  1877. X
  1878. Xvoid LogAndSign::operator*=(Real x)
  1879. X{
  1880. X   if (x > 0.0) { log_value += log(x); }
  1881. X   else if (x < 0.0) { log_value += log(-x); sign = -sign; }
  1882. X   else sign = 0;
  1883. X}
  1884. X
  1885. XReal LogAndSign::Value() const { return sign * exp(log_value); }
  1886. X
  1887. XLogAndSign::LogAndSign(Real f)
  1888. X{
  1889. X   if (f == 0.0) { log_value = 0.0; sign = 0; return; }
  1890. X   else if (f < 0.0) { sign = -1; f = -f; }
  1891. X   else sign = 1;
  1892. X   log_value = log(f);
  1893. X}
  1894. X
  1895. XLogAndSign DiagonalMatrix::LogDeterminant() const
  1896. X{
  1897. X   REPORT
  1898. X   int i = nrows; LogAndSign sum; Real* s = store;
  1899. X   while (i--) sum *= *s++;
  1900. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1901. X}
  1902. X
  1903. XLogAndSign LowerTriangularMatrix::LogDeterminant() const
  1904. X{
  1905. X   REPORT
  1906. X   int i = nrows; LogAndSign sum; Real* s = store; int j = 2;
  1907. X   while (i--) { sum *= *s; s += j++; }
  1908. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1909. X}
  1910. X
  1911. XLogAndSign UpperTriangularMatrix::LogDeterminant() const
  1912. X{
  1913. X   REPORT
  1914. X   int i = nrows; LogAndSign sum; Real* s = store;
  1915. X   while (i) { sum *= *s; s += i--; }
  1916. X   ((GeneralMatrix&)*this).tDelete(); return sum;
  1917. X}
  1918. X
  1919. XLogAndSign BaseMatrix::LogDeterminant() const
  1920. X{
  1921. X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  1922. X   LogAndSign sum = gm->LogDeterminant(); return sum;
  1923. X}
  1924. X
  1925. XLogAndSign GeneralMatrix::LogDeterminant() const
  1926. X{
  1927. X   REPORT
  1928. X   Tracer tr("Determinant");
  1929. X   if (nrows != ncols) Throw(NotSquareException(*this));
  1930. X   CroutMatrix C(*this); return C.LogDeterminant();
  1931. X}
  1932. X
  1933. XLogAndSign CroutMatrix::LogDeterminant() const
  1934. X{
  1935. X   REPORT
  1936. X   if (sing) return 0.0;
  1937. X   int i = nrows; int dd = i+1; LogAndSign sum; Real* s = store;
  1938. X   while (i--) { sum *= *s; s += dd; }
  1939. X   if (!d) sum.ChangeSign(); return sum;
  1940. X
  1941. X}
  1942. X
  1943. XLinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
  1944. X: gm( ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver() )
  1945. X{
  1946. X   if (gm==&bm) { REPORT  gm = gm->Image(); } 
  1947. X   // want a copy if  *gm is actually bm
  1948. X   else { REPORT  gm->Protect(); }
  1949. X}
  1950. X
  1951. END_OF_FILE
  1952. if test 8288 -ne `wc -c <'newmat8.cxx'`; then
  1953.     echo shar: \"'newmat8.cxx'\" unpacked with wrong size!
  1954. fi
  1955. # end of 'newmat8.cxx'
  1956. fi
  1957. if test -f 'newmat9.cxx' -a "${1}" != "-c" ; then 
  1958.   echo shar: Will not clobber existing file \"'newmat9.cxx'\"
  1959. else
  1960. echo shar: Extracting \"'newmat9.cxx'\" \(954 characters\)
  1961. sed "s/^X//" >'newmat9.cxx' <<'END_OF_FILE'
  1962. X//$$ newmat9.cxx         Input and output
  1963. X
  1964. X// Copyright (C) 1991,2,3: R B Davies
  1965. X
  1966. X
  1967. X#define WANT_STREAM
  1968. X
  1969. X#include "include.h"
  1970. X
  1971. X#include "newmat.h"
  1972. X#include "newmatrc.h"
  1973. X#include "newmatio.h"
  1974. X
  1975. X//#define REPORT { static ExeCounter ExeCount(__LINE__,9); ExeCount++; }
  1976. X
  1977. X#define REPORT {}
  1978. X
  1979. Xostream& operator<<(ostream& s, const BaseMatrix& X)
  1980. X{
  1981. X   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); operator<<(s, *gm);
  1982. X   gm->tDelete(); return s;
  1983. X}
  1984. X
  1985. X
  1986. Xostream& operator<<(ostream& s, const GeneralMatrix& X)
  1987. X{
  1988. X   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
  1989. X   int w = s.width();  int nr = X.Nrows();  long f = s.flags();
  1990. X   s.setf(ios::fixed, ios::floatfield);
  1991. X   for (int i=1; i<=nr; i++)
  1992. X   {
  1993. X      int skip = mr.skip;  int storage = mr.storage;
  1994. X      Real* store = mr.store+skip;  skip *= w+1;
  1995. X      while (skip--) s << " ";
  1996. X      while (storage--) s << setw(w) << *store++ << " ";
  1997. X      mr.Next();  s << "\n";
  1998. X   }
  1999. X   s << flush;  s.flags(f); return s;
  2000. X}
  2001. X
  2002. END_OF_FILE
  2003. if test 954 -ne `wc -c <'newmat9.cxx'`; then
  2004.     echo shar: \"'newmat9.cxx'\" unpacked with wrong size!
  2005. fi
  2006. # end of 'newmat9.cxx'
  2007. fi
  2008. if test -f 'newmatrm.cxx' -a "${1}" != "-c" ; then 
  2009.   echo shar: Will not clobber existing file \"'newmatrm.cxx'\"
  2010. else
  2011. echo shar: Extracting \"'newmatrm.cxx'\" \(3388 characters\)
  2012. sed "s/^X//" >'newmatrm.cxx' <<'END_OF_FILE'
  2013. X//$$newmatrm.cxx                         rectangular matrix operations
  2014. X
  2015. X// Copyright (C) 1991,2,3: R B Davies
  2016. X
  2017. X#include "include.h"
  2018. X
  2019. X#include "newmat.h"
  2020. X#include "newmatrm.h"
  2021. X
  2022. X
  2023. X// operations on rectangular matrices
  2024. X
  2025. X
  2026. Xvoid RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
  2027. X{
  2028. X   RectMatrixRowCol::Reset
  2029. X      ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
  2030. X}
  2031. X
  2032. Xvoid RectMatrixRow::Reset (const Matrix& M, int row)
  2033. X{
  2034. X   RectMatrixRowCol::Reset( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() );
  2035. X}
  2036. X
  2037. Xvoid RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
  2038. X{
  2039. X   RectMatrixRowCol::Reset
  2040. X      ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
  2041. X}
  2042. X
  2043. Xvoid RectMatrixCol::Reset (const Matrix& M, int col)
  2044. X   { RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 ); }
  2045. X
  2046. X
  2047. XReal RectMatrixRowCol::SumSquare() const
  2048. X{
  2049. X   long_Real sum = 0.0; int i = n; Real* s = store; int d = spacing;
  2050. X   while (i--) { sum += (long_Real)*s * *s; s += d; }
  2051. X   return (Real)sum;
  2052. X}
  2053. X
  2054. XReal RectMatrixRowCol::operator*(const RectMatrixRowCol& rmrc) const
  2055. X{
  2056. X   long_Real sum = 0.0; int i = n;
  2057. X   Real* s = store; int d = spacing;
  2058. X   Real* s1 = rmrc.store; int d1 = rmrc.spacing;
  2059. X   if (i!=rmrc.n)
  2060. X   {
  2061. X      Tracer tr("newmatrm");
  2062. X      Throw(InternalException("Dimensions differ in *"));
  2063. X   }
  2064. X   while (i--) { sum += (long_Real)*s * *s1; s += d; s1 += d1; }
  2065. X   return (Real)sum;
  2066. X}
  2067. X
  2068. Xvoid RectMatrixRowCol::AddScaled(const RectMatrixRowCol& rmrc, Real r)
  2069. X{
  2070. X   int i = n; Real* s = store; int d = spacing;
  2071. X   Real* s1 = rmrc.store; int d1 = rmrc.spacing;
  2072. X   if (i!=rmrc.n)
  2073. X   {
  2074. X      Tracer tr("newmatrm");
  2075. X      Throw(InternalException("Dimensions differ in AddScaled"));
  2076. X   }
  2077. X   while (i--) { *s += *s1 * r; s += d; s1 += d1; }
  2078. X}
  2079. X
  2080. Xvoid RectMatrixRowCol::Divide(const RectMatrixRowCol& rmrc, Real r)
  2081. X{
  2082. X   int i = n; Real* s = store; int d = spacing;
  2083. X   Real* s1 = rmrc.store; int d1 = rmrc.spacing;
  2084. X   if (i!=rmrc.n)
  2085. X   {
  2086. X      Tracer tr("newmatrm");
  2087. X      Throw(InternalException("Dimensions differ in Divide"));
  2088. X   }
  2089. X   while (i--) { *s = *s1 / r; s += d; s1 += d1; }
  2090. X}
  2091. X
  2092. Xvoid RectMatrixRowCol::Divide(Real r)
  2093. X{
  2094. X   int i = n; Real* s = store; int d = spacing;
  2095. X   while (i--) { *s /= r; s += d; }
  2096. X}
  2097. X
  2098. Xvoid RectMatrixRowCol::Negate()
  2099. X{
  2100. X   int i = n; Real* s = store; int d = spacing;
  2101. X   while (i--) { *s = - *s; s += d; }
  2102. X}
  2103. X
  2104. Xvoid RectMatrixRowCol::Zero()
  2105. X{
  2106. X   int i = n; Real* s = store; int d = spacing;
  2107. X   while (i--) { *s = 0.0; s += d; }
  2108. X}
  2109. X
  2110. Xvoid ComplexScale(RectMatrixCol& U, RectMatrixCol& V, Real x, Real y)
  2111. X{
  2112. X   int n = U.n;
  2113. X   if (n != V.n)
  2114. X   {
  2115. X      Tracer tr("newmatrm");
  2116. X      Throw(InternalException("Dimensions differ in ComplexScale"));
  2117. X   }
  2118. X   Real* u = U.store; Real* v = V.store; 
  2119. X   int su = U.spacing; int sv = V.spacing;
  2120. X   while (n--)
  2121. X   {
  2122. X      Real z = *u * x - *v * y;  *v =  *u * y + *v * x;  *u = z;
  2123. X      u += su;  v += sv;
  2124. X   }
  2125. X}
  2126. X
  2127. Xvoid Rotate(RectMatrixCol& U, RectMatrixCol& V, Real tau, Real s)
  2128. X{
  2129. X   //  (U, V) = (U, V) * (c, s)  where  tau = s/(1+c), c^2 + s^2 = 1
  2130. X   int n = U.n;
  2131. X   if (n != V.n)
  2132. X   {
  2133. X      Tracer tr("newmatrm");
  2134. X      Throw(InternalException("Dimensions differ in Rotate"));
  2135. X   }
  2136. X   Real* u = U.store; Real* v = V.store; 
  2137. X   int su = U.spacing; int sv = V.spacing;
  2138. X   while (n--)
  2139. X   {
  2140. X      Real zu = *u; Real zv = *v;
  2141. X      *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
  2142. X      u += su;  v += sv;
  2143. X   }
  2144. X}
  2145. X
  2146. X
  2147. END_OF_FILE
  2148. if test 3388 -ne `wc -c <'newmatrm.cxx'`; then
  2149.     echo shar: \"'newmatrm.cxx'\" unpacked with wrong size!
  2150. fi
  2151. # end of 'newmatrm.cxx'
  2152. fi
  2153. echo shar: End of archive 5 \(of 8\).
  2154. cp /dev/null ark5isdone
  2155. MISSING=""
  2156. for I in 1 2 3 4 5 6 7 8 ; do
  2157.     if test ! -f ark${I}isdone ; then
  2158.     MISSING="${MISSING} ${I}"
  2159.     fi
  2160. done
  2161. if test "${MISSING}" = "" ; then
  2162.     echo You have unpacked all 8 archives.
  2163.     rm -f ark[1-9]isdone
  2164. else
  2165.     echo You still need to unpack the following archives:
  2166.     echo "        " ${MISSING}
  2167. fi
  2168. ##  End of shell archive.
  2169. exit 0
  2170.  
  2171. exit 0 # Just in case...
  2172.