home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1994 March / Source_Code_CD-ROM_Walnut_Creek_March_1994.iso / compsrcs / misc / volume34 / newmat07 / part04 < prev    next >
Encoding:
Text File  |  1993-01-10  |  58.6 KB  |  2,083 lines
  1. Newsgroups: comp.sources.misc
  2. From: robertd@kauri.vuw.ac.nz (Robert Davies)
  3. Subject: v34i110:  newmat07 - A matrix package in C++, Part04/08
  4. Message-ID: <1993Jan11.153136.2384@sparky.imd.sterling.com>
  5. X-Md4-Signature: 3ad844ba263547ee914107dc1613b003
  6. Date: Mon, 11 Jan 1993 15:31:36 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 110
  11. Archive-name: newmat07/part04
  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 4 (of 8)."
  22. # Contents:  newmat1.cxx newmat2.cxx newmat3.cxx newmat4.cxx
  23. #   newmatex.cxx
  24. # Wrapped by robert@kea on Sun Jan 10 23:57:52 1993
  25. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  26. if test -f 'newmat1.cxx' -a "${1}" != "-c" ; then 
  27.   echo shar: Will not clobber existing file \"'newmat1.cxx'\"
  28. else
  29. echo shar: Extracting \"'newmat1.cxx'\" \(3276 characters\)
  30. sed "s/^X//" >'newmat1.cxx' <<'END_OF_FILE'
  31. X//$$ newmat1.cxx   Matrix type functions
  32. X
  33. X// Copyright (C) 1991,2,3: R B Davies
  34. X
  35. X
  36. X#include "include.h"
  37. X
  38. X#include "newmat.h"
  39. X
  40. X
  41. X
  42. X/************************* MatrixType functions *****************************/
  43. X
  44. X
  45. XMatrixType MatrixType::operator*(const MatrixType& mt) const
  46. X{
  47. X   int a = attribute & mt.attribute;
  48. X   if ((a & (Upper+Lower)) == (Upper+Lower)) return Dg;
  49. X   else return a & ~Symmetric;
  50. X}
  51. X
  52. XMatrixType MatrixType::t() const
  53. X{
  54. X   int a = attribute & ~(Upper+Lower);
  55. X   if (attribute & Upper) a |= Lower;
  56. X   if (attribute & Lower) a |= Upper;
  57. X   return a;
  58. X}
  59. X
  60. XMatrixType MatrixType::MultRHS() const
  61. X{
  62. X   if ((attribute & (Upper+Lower)) == (Upper+Lower)) return Dg;
  63. X   else return attribute & ~Symmetric;
  64. X}
  65. X
  66. X
  67. X
  68. XMatrixType::operator char*() const
  69. X{
  70. X// make a string with the name of matrix with the given attributes
  71. X   switch (attribute)
  72. X   {
  73. X   case Valid:                              return "Rect ";
  74. X   case Valid+Symmetric:                    return "Sym  ";
  75. X   case Valid+Band:                         return "Band ";
  76. X   case Valid+Symmetric+Band:               return "SmBnd";
  77. X   case Valid+Upper:                        return "UT   ";
  78. X   case Valid+Upper+Lower:
  79. X   case Valid+Band+Upper+Lower:
  80. X   case Valid+Symmetric+Upper:
  81. X   case Valid+Symmetric+Band+Upper:
  82. X   case Valid+Symmetric+Lower:
  83. X   case Valid+Symmetric+Band+Lower:
  84. X   case Valid+Symmetric+Upper+Lower:
  85. X   case Valid+Symmetric+Band+Upper+Lower:   return "Diag ";
  86. X   case Valid+Band+Upper:                   return "UpBnd";
  87. X   case Valid+Lower:                        return "LT   ";
  88. X   case Valid+Band+Lower:                   return "LwBnd";
  89. X   default:
  90. X      if (!(attribute & Valid))             return "UnSp ";
  91. X      if (attribute & LUDeco)
  92. X         return (attribute & Band) ?     "BndLU" : "Crout";
  93. X                                            return "?????";
  94. X   }
  95. X}
  96. X
  97. X
  98. XGeneralMatrix* MatrixType::New(int nr, int nc, BaseMatrix* bm) const
  99. X{
  100. X// make a new matrix with the given attributes
  101. X
  102. X   Tracer tr("New"); GeneralMatrix* gm;
  103. X   switch (attribute)
  104. X   {
  105. X   case Valid:
  106. X      if (nc==1) { gm = new ColumnVector(nr); break; }
  107. X      if (nr==1) { gm = new RowVector(nc); break; }
  108. X      gm = new Matrix(nr, nc); break;
  109. X
  110. X   case Valid+Symmetric:
  111. X      gm = new SymmetricMatrix(nr); break;
  112. X
  113. X   case Valid+Band:
  114. X      {
  115. X         MatrixBandWidth bw = bm->BandWidth();
  116. X         gm = new BandMatrix(nr,bw.lower,bw.upper); break;
  117. X      }
  118. X
  119. X   case Valid+Symmetric+Band:
  120. X      gm = new SymmetricBandMatrix(nr,bm->BandWidth().lower); break;
  121. X
  122. X   case Valid+Upper:
  123. X      gm = new UpperTriangularMatrix(nr); break;
  124. X
  125. X   case Valid+Upper+Lower:
  126. X   case Valid+Band+Upper+Lower:
  127. X   case Valid+Symmetric+Upper:
  128. X   case Valid+Symmetric+Band+Upper:
  129. X   case Valid+Symmetric+Lower:
  130. X   case Valid+Symmetric+Band+Lower:
  131. X   case Valid+Symmetric+Upper+Lower:
  132. X   case Valid+Symmetric+Band+Upper+Lower:
  133. X      gm = new DiagonalMatrix(nr); break;
  134. X
  135. X   case Valid+Band+Upper:
  136. X      gm = new UpperBandMatrix(nr,bm->BandWidth().upper); break;
  137. X
  138. X   case Valid+Lower:
  139. X      gm = new LowerTriangularMatrix(nr); break;
  140. X
  141. X   case Valid+Band+Lower:
  142. X      gm = new LowerBandMatrix(nr,bm->BandWidth().lower); break;
  143. X
  144. X   default:
  145. X      Throw(ProgramException("Invalid matrix type"));
  146. X   }
  147. X   
  148. X   MatrixErrorNoSpace(gm); gm->Protect(); return gm;
  149. X}
  150. X
  151. END_OF_FILE
  152. if test 3276 -ne `wc -c <'newmat1.cxx'`; then
  153.     echo shar: \"'newmat1.cxx'\" unpacked with wrong size!
  154. fi
  155. # end of 'newmat1.cxx'
  156. fi
  157. if test -f 'newmat2.cxx' -a "${1}" != "-c" ; then 
  158.   echo shar: Will not clobber existing file \"'newmat2.cxx'\"
  159. else
  160. echo shar: Extracting \"'newmat2.cxx'\" \(10973 characters\)
  161. sed "s/^X//" >'newmat2.cxx' <<'END_OF_FILE'
  162. X//$$ newmat2.cxx      Matrix row and column operations
  163. X
  164. X// Copyright (C) 1991,2,3: R B Davies
  165. X
  166. X#define WANT_MATH
  167. X
  168. X#include "include.h"
  169. X
  170. X#include "newmat.h"
  171. X#include "newmatrc.h"
  172. X
  173. X//#define REPORT { static ExeCounter ExeCount(__LINE__,2); ExeCount++; }
  174. X
  175. X#define REPORT {}
  176. X
  177. X//#define MONITOR(what,storage,store) \
  178. X//   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  179. X
  180. X#define MONITOR(what,store,storage) {}
  181. X
  182. X/************************** Matrix Row/Col functions ************************/
  183. X
  184. Xvoid MatrixRowCol::Add(const MatrixRowCol& mrc)
  185. X{
  186. X   REPORT
  187. X   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  188. X   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  189. X   if (l<=0) return;
  190. X   Real* elx=store+f; Real* el=mrc.store+f;
  191. X   while (l--) *elx++ += *el++;
  192. X}
  193. X
  194. Xvoid MatrixRowCol::AddScaled(const MatrixRowCol& mrc, Real x)
  195. X{
  196. X   REPORT
  197. X   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  198. X   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  199. X   if (l<=0) return;
  200. X   Real* elx=store+f; Real* el=mrc.store+f;
  201. X   while (l--) *elx++ += *el++ * x;
  202. X}
  203. X
  204. Xvoid MatrixRowCol::Sub(const MatrixRowCol& mrc)
  205. X{
  206. X   REPORT
  207. X   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  208. X   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  209. X   if (l<=0) return;
  210. X   Real* elx=store+f; Real* el=mrc.store+f;
  211. X   while (l--) *elx++ -= *el++;
  212. X}
  213. X
  214. Xvoid MatrixRowCol::Inject(const MatrixRowCol& mrc)
  215. X// copy stored elements only
  216. X{
  217. X   REPORT
  218. X   int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  219. X   if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  220. X   if (l<=0) return;
  221. X   Real* elx = store+f; Real* ely = mrc.store+f;
  222. X   while (l--) *elx++ = *ely++;
  223. X}
  224. X
  225. XReal DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  226. X{
  227. X   REPORT                                         // not accessed
  228. X   int f = mrc1.skip; int f2 = mrc2.skip;
  229. X   int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
  230. X   if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
  231. X   if (l<=0) return 0.0;
  232. X   
  233. X   Real* el1=mrc1.store+f; Real* el2=mrc2.store+f;
  234. X   Real sum = 0.0;
  235. X   while (l--) sum += *el1++ * *el2++;
  236. X   return sum;
  237. X}
  238. X
  239. Xvoid MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  240. X{
  241. X   int f = skip; int l = skip + storage;
  242. X   int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  243. X   if (f1<f) f1=f; if (l1>l) l1=l;
  244. X   int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  245. X   if (f2<f) f2=f; if (l2>l) l2=l;
  246. X   Real* el = store + f;
  247. X   Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
  248. X   if (f1<f2)
  249. X   {
  250. X      register int i = f1-f; while (i--) *el++ = 0.0;
  251. X      if (l1<=f2)                              // disjoint
  252. X      {
  253. X     REPORT                                // not accessed
  254. X         i = l1-f1;     while (i--) *el++ = *el1++;
  255. X         i = f2-l1;     while (i--) *el++ = 0.0;
  256. X         i = l2-f2;     while (i--) *el++ = *el2++;
  257. X         i = l-l2;      while (i--) *el++ = 0.0;
  258. X      }
  259. X      else
  260. X      {
  261. X         i = f2-f1;    while (i--) *el++ = *el1++;
  262. X         if (l1<=l2)
  263. X         {
  264. X            REPORT
  265. X            i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
  266. X            i = l2-l1; while (i--) *el++ = *el2++;
  267. X            i = l-l2;  while (i--) *el++ = 0.0;
  268. X         }
  269. X         else
  270. X         {
  271. X            REPORT
  272. X            i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
  273. X            i = l1-l2; while (i--) *el++ = *el1++;
  274. X            i = l-l1;  while (i--) *el++ = 0.0;
  275. X         }
  276. X      }
  277. X   }
  278. X   else
  279. X   {
  280. X      register int i = f2-f; while (i--) *el++ = 0.0;
  281. X      if (l2<=f1)                              // disjoint
  282. X      {
  283. X     REPORT                                // not accessed
  284. X         i = l2-f2;     while (i--) *el++ = *el2++;
  285. X         i = f1-l2;     while (i--) *el++ = 0.0;
  286. X     i = l1-f1;     while (i--) *el++ = *el1++;
  287. X     i = l-l1;      while (i--) *el++ = 0.0;
  288. X      }
  289. X      else
  290. X      {
  291. X         i = f1-f2;    while (i--) *el++ = *el2++;
  292. X         if (l2<=l1)
  293. X         {
  294. X        REPORT
  295. X            i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
  296. X            i = l1-l2; while (i--) *el++ = *el1++;
  297. X            i = l-l1;  while (i--) *el++ = 0.0;
  298. X         }
  299. X         else
  300. X         {
  301. X        REPORT
  302. X            i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
  303. X            i = l2-l1; while (i--) *el++ = *el2++;
  304. X            i = l-l2;  while (i--) *el++ = 0.0;
  305. X         }
  306. X      }
  307. X   }
  308. X}
  309. X
  310. Xvoid MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  311. X{
  312. X   int f = skip; int l = skip + storage;
  313. X   int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  314. X   if (f1<f) f1=f; if (l1>l) l1=l;
  315. X   int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  316. X   if (f2<f) f2=f; if (l2>l) l2=l;
  317. X   Real* el = store + f;
  318. X   Real* el1 = mrc1.store+f1; Real* el2 = mrc2.store+f2;
  319. X   if (f1<f2)
  320. X   {
  321. X      register int i = f1-f; while (i--) *el++ = 0.0;
  322. X      if (l1<=f2)                              // disjoint
  323. X      {
  324. X     REPORT                                // not accessed
  325. X         i = l1-f1;     while (i--) *el++ = *el1++;
  326. X         i = f2-l1;     while (i--) *el++ = 0.0;
  327. X         i = l2-f2;     while (i--) *el++ = - *el2++;
  328. X         i = l-l2;      while (i--) *el++ = 0.0;
  329. X      }
  330. X      else
  331. X      {
  332. X         i = f2-f1;    while (i--) *el++ = *el1++;
  333. X         if (l1<=l2)
  334. X         {
  335. X        REPORT
  336. X            i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
  337. X            i = l2-l1; while (i--) *el++ = - *el2++;
  338. X            i = l-l2;  while (i--) *el++ = 0.0;
  339. X         }
  340. X         else
  341. X         {
  342. X        REPORT
  343. X            i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
  344. X            i = l1-l2; while (i--) *el++ = *el1++;
  345. X            i = l-l1;  while (i--) *el++ = 0.0;
  346. X         }
  347. X      }
  348. X   }
  349. X   else
  350. X   {
  351. X      register int i = f2-f; while (i--) *el++ = 0.0;
  352. X      if (l2<=f1)                              // disjoint
  353. X      {
  354. X     REPORT                                // not accessed
  355. X         i = l2-f2;     while (i--) *el++ = - *el2++;
  356. X         i = f1-l2;     while (i--) *el++ = 0.0;
  357. X         i = l1-f1;     while (i--) *el++ = *el1++;
  358. X         i = l-l1;      while (i--) *el++ = 0.0;
  359. X      }
  360. X      else
  361. X      {
  362. X         i = f1-f2;    while (i--) *el++ = - *el2++;
  363. X         if (l2<=l1)
  364. X         {
  365. X        REPORT
  366. X            i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
  367. X            i = l1-l2; while (i--) *el++ = *el1++;
  368. X            i = l-l1;  while (i--) *el++ = 0.0;
  369. X         }
  370. X         else
  371. X         {
  372. X            REPORT
  373. X            i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
  374. X            i = l2-l1; while (i--) *el++ = - *el2++;
  375. X            i = l-l2;  while (i--) *el++ = 0.0;
  376. X         }
  377. X      }
  378. X   }
  379. X}
  380. X
  381. X
  382. Xvoid MatrixRowCol::Add(const MatrixRowCol& mrc1, Real x)
  383. X{
  384. X   REPORT
  385. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  386. X   if (f < skip) { f = skip; if (l < f) l = f; }
  387. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  388. X
  389. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  390. X
  391. X   int l1 = f-skip;  while (l1--) *elx++ = x;
  392. X       l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
  393. X       lx -= l;      while (lx--) *elx++ = x;
  394. X}
  395. X
  396. Xvoid MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
  397. X{
  398. X   REPORT
  399. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  400. X   if (f < skip) { f = skip; if (l < f) l = f; }
  401. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  402. X
  403. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  404. X
  405. X   int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
  406. X       l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
  407. X       lx -= l;      while (lx--) { *elx = - *elx; elx++; }
  408. X}
  409. X
  410. Xvoid MatrixRowCol::Copy(const MatrixRowCol& mrc1)
  411. X{
  412. X   REPORT
  413. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  414. X   if (f < skip) { f = skip; if (l < f) l = f; }
  415. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  416. X
  417. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  418. X
  419. X   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  420. X       l1 = l-f;     while (l1--) *elx++ = *ely++;
  421. X       lx -= l;      while (lx--) *elx++ = 0.0;
  422. X}
  423. X
  424. Xvoid MatrixRowCol::CopyCheck(const MatrixRowCol& mrc1)
  425. X// Throw an exception this would lead to a loss of data
  426. X{
  427. X   REPORT
  428. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  429. X   if (f < skip || l > lx) Throw(ProgramException("Illegal Conversion"));
  430. X
  431. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  432. X
  433. X   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  434. X       l1 = l-f;     while (l1--) *elx++ = *ely++;
  435. X       lx -= l;      while (lx--) *elx++ = 0.0;
  436. X}
  437. X
  438. Xvoid MatrixRowCol::Negate(const MatrixRowCol& mrc1)
  439. X{
  440. X   REPORT
  441. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  442. X   if (f < skip) { f = skip; if (l < f) l = f; }
  443. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  444. X
  445. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  446. X
  447. X   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  448. X       l1 = l-f;     while (l1--) *elx++ = - *ely++;
  449. X       lx -= l;      while (lx--) *elx++ = 0.0;
  450. X}
  451. X
  452. Xvoid MatrixRowCol::Multiply(const MatrixRowCol& mrc1, Real s)
  453. X{
  454. X   REPORT
  455. X   int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  456. X   if (f < skip) { f = skip; if (l < f) l = f; }
  457. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  458. X
  459. X   Real* elx = store+skip; Real* ely = mrc1.store+f;
  460. X
  461. X   int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  462. X       l1 = l-f;     while (l1--) *elx++ = *ely++ * s;
  463. X       lx -= l;      while (lx--) *elx++ = 0.0;
  464. X}
  465. X
  466. Xvoid DiagonalMatrix::Solver(MatrixRowCol& mrc, const MatrixRowCol& mrc1)
  467. X{
  468. X   REPORT
  469. X   int f = mrc1.skip; int f0 = mrc.skip;
  470. X   int l = f + mrc1.storage; int lx = f0 + mrc.storage;
  471. X   if (f < f0) { f = f0; if (l < f) l = f; }
  472. X   if (l > lx) { l = lx; if (f > lx) f = lx; }
  473. X
  474. X   Real* elx = mrc.store+f0; Real* eld = store+f;
  475. X
  476. X   int l1 = f-f0;    while (l1--) *elx++ = 0.0;
  477. X       l1 = l-f;     while (l1--) *elx++ /= *eld++;
  478. X       lx -= l;      while (lx--) *elx++ = 0.0;
  479. X   // Solver makes sure input and output point to same memory
  480. X}
  481. X
  482. Xvoid MatrixRowCol::Copy(const Real*& r)
  483. X{
  484. X   REPORT
  485. X   Real* elx = store+skip; const Real* ely = r+skip; r += length;
  486. X   int l = storage; while (l--) *elx++ = *ely++;
  487. X}
  488. X
  489. Xvoid MatrixRowCol::Copy(Real r)
  490. X{
  491. X   REPORT
  492. X   Real* elx = store+skip; int l = storage; while (l--) *elx++ = r;
  493. X}
  494. X
  495. XReal MatrixRowCol::SumAbsoluteValue()
  496. X{
  497. X   REPORT
  498. X   Real sum = 0.0; Real* elx = store+skip; int l = storage;
  499. X   while (l--) sum += fabs(*elx++);
  500. X   return sum;
  501. X}
  502. X
  503. Xvoid MatrixRowCol::SubRowCol(MatrixRowCol& mrc, int skip1, int l1) const
  504. X{
  505. X   mrc.length = l1;  mrc.store = store + skip1;  int d = skip - skip1;
  506. X   mrc.skip = (d < 0) ? 0 : d;  d = skip + storage - skip1;
  507. X   d = ((l1 < d) ? l1 : d) - mrc.skip;  mrc.storage = (d < 0) ? 0 : d;
  508. X   mrc.cw = 0;
  509. X}
  510. X
  511. XMatrixRowCol::~MatrixRowCol()
  512. X{
  513. X   int t1 = +(cw*IsACopy); int t2 = !(cw*StoreHere);
  514. X   if (t1 && t2)
  515. X   {
  516. X      Real* f = store+skip;
  517. X      MONITOR_REAL_DELETE("Free    (RowCol)",storage,f) 
  518. X#ifdef Version21
  519. X      delete [] f;
  520. X#else
  521. X      delete [storage] f;
  522. X#endif
  523. X   }
  524. X}
  525. X
  526. XMatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }
  527. X
  528. XMatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }
  529. X
  530. END_OF_FILE
  531. if test 10973 -ne `wc -c <'newmat2.cxx'`; then
  532.     echo shar: \"'newmat2.cxx'\" unpacked with wrong size!
  533. fi
  534. # end of 'newmat2.cxx'
  535. fi
  536. if test -f 'newmat3.cxx' -a "${1}" != "-c" ; then 
  537.   echo shar: Will not clobber existing file \"'newmat3.cxx'\"
  538. else
  539. echo shar: Extracting \"'newmat3.cxx'\" \(14262 characters\)
  540. sed "s/^X//" >'newmat3.cxx' <<'END_OF_FILE'
  541. X//$$ newmat3.cxx        Matrix get and restore rows and columns
  542. X
  543. X// Copyright (C) 1991,2,3: R B Davies
  544. X
  545. X
  546. X#include "include.h"
  547. X
  548. X#include "newmat.h"
  549. X#include "newmatrc.h"
  550. X
  551. X//#define REPORT { static ExeCounter ExeCount(__LINE__,3); ExeCount++; }
  552. X
  553. X#define REPORT {}
  554. X
  555. X//#define MONITOR(what,storage,store) \
  556. X//   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  557. X
  558. X#define MONITOR(what,store,storage) {}
  559. X
  560. X
  561. X// Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
  562. X//
  563. X// LoadOnEntry:
  564. X//    Load data into MatrixRow or Col dummy array under GetRow or GetCol
  565. X// StoreOnExit:
  566. X//    Restore data to original matrix under RestoreRow or RestoreCol
  567. X// IsACopy:
  568. X//    Set by GetRow/Col: MatrixRow or Col array is a copy
  569. X// DirectPart:
  570. X//    Load or restore only part directly stored; must be set with StoreOnExit
  571. X//    Still have decide  how to handle this with symmetric
  572. X// StoreHere:
  573. X//    used in columns only - store data at supplied storage address, adjusted
  574. X//    for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
  575. X//    zeros.
  576. X
  577. X
  578. X// These will work as a default
  579. X// but need to override NextRow for efficiency
  580. X
  581. X// Assume pointer arithmetic works for pointers out of range - not strict C++.
  582. X
  583. X
  584. Xvoid GeneralMatrix::NextRow(MatrixRowCol& mrc)
  585. X{
  586. X   REPORT
  587. X   if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreRow(mrc); }
  588. X   if (+(mrc.cw*IsACopy))
  589. X   {
  590. X      REPORT
  591. X      Real* s = mrc.store + mrc.skip;
  592. X      MONITOR_REAL_DELETE("Free   (NextRow)",mrc.storage,s)
  593. X#ifdef Version21
  594. X      delete [] s;
  595. X#else
  596. X      delete [mrc.storage] s;
  597. X#endif
  598. X   }
  599. X   mrc.rowcol++;
  600. X   if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
  601. X   else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  602. X}
  603. X
  604. Xvoid GeneralMatrix::NextCol(MatrixRowCol& mrc)
  605. X{
  606. X   REPORT                                                // 423
  607. X   if (+(mrc.cw*StoreOnExit)) { REPORT this->RestoreCol(mrc); }
  608. X   int t1 = +(mrc.cw*IsACopy); int t2 = !(mrc.cw*StoreHere);
  609. X   if ( t1 && t2 )
  610. X   {
  611. X      REPORT                                             // not accessed
  612. X      Real* s = mrc.store + mrc.skip;
  613. X      MONITOR_REAL_DELETE("Free   (NextCol)",mrc.storage,s) 
  614. X#ifdef Version21
  615. X      delete [] s;
  616. X#else
  617. X      delete [mrc.storage] s;
  618. X#endif
  619. X   }
  620. X   mrc.rowcol++;
  621. X   if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  622. X   else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  623. X}
  624. X
  625. X
  626. X// routines for matrix
  627. X
  628. Xvoid Matrix::GetRow(MatrixRowCol& mrc)
  629. X{
  630. X   REPORT
  631. X   mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
  632. X   mrc.store=store+mrc.rowcol*ncols;
  633. X}
  634. X
  635. X
  636. Xvoid Matrix::GetCol(MatrixRowCol& mrc)
  637. X{
  638. X   REPORT
  639. X   mrc.skip=0; mrc.storage=nrows; int t1 = !(mrc.cw*StoreHere);
  640. X   if ( ncols==1 && t1 )
  641. X      { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
  642. X   else
  643. X   {
  644. X      mrc.cw+=IsACopy; Real* ColCopy;
  645. X      if ( !(mrc.cw*StoreHere) )
  646. X      {
  647. X         REPORT
  648. X         ColCopy = new Real [nrows]; MatrixErrorNoSpace(ColCopy);
  649. X         MONITOR_REAL_NEW("Make (MatGetCol)",nrows,ColCopy)
  650. X         mrc.store = ColCopy;
  651. X      }
  652. X      else { REPORT ColCopy = mrc.store; }
  653. X      if (+(mrc.cw*LoadOnEntry))
  654. X      {
  655. X         REPORT
  656. X         Real* Mstore = store+mrc.rowcol; int i=nrows;
  657. X         while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  658. X      }
  659. X   }
  660. X}
  661. X
  662. Xvoid Matrix::RestoreCol(MatrixRowCol& mrc)
  663. X{
  664. X//  if (mrc.cw*StoreOnExit)
  665. X   REPORT                                   // 429
  666. X   if (+(mrc.cw*IsACopy))
  667. X   {
  668. X      REPORT                                // 426
  669. X      Real* Mstore = store+mrc.rowcol; int i=nrows; Real* Cstore = mrc.store;
  670. X      while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  671. X   }
  672. X}
  673. X
  674. Xvoid Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  675. X
  676. Xvoid Matrix::NextCol(MatrixRowCol& mrc)
  677. X{
  678. X   REPORT                                        // 632
  679. X   if (+(mrc.cw*StoreOnExit)) { REPORT RestoreCol(mrc); }
  680. X   mrc.rowcol++;
  681. X   if (mrc.rowcol<ncols)
  682. X   {
  683. X      if (+(mrc.cw*LoadOnEntry))
  684. X      {
  685. X     REPORT
  686. X         Real* ColCopy = mrc.store;
  687. X         Real* Mstore = store+mrc.rowcol; int i=nrows;
  688. X         while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  689. X      }
  690. X   }
  691. X   else { REPORT mrc.cw -= StoreOnExit; }
  692. X}
  693. X
  694. X// routines for diagonal matrix
  695. X
  696. Xvoid DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  697. X{
  698. X   REPORT
  699. X   mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  700. X}
  701. X
  702. Xvoid DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  703. X{
  704. X   REPORT 
  705. X   mrc.skip=mrc.rowcol; mrc.storage=1;
  706. X   if (+(mrc.cw*StoreHere))
  707. X      { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  708. X   else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  709. X}
  710. X
  711. Xvoid DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  712. X                              // 800
  713. Xvoid DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  714. X{
  715. X   REPORT
  716. X   if (+(mrc.cw*StoreHere))
  717. X   {
  718. X      if (+(mrc.cw*StoreOnExit))
  719. X         { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  720. X      mrc.IncrDiag();
  721. X      int t1 = +(mrc.cw*LoadOnEntry);
  722. X      if (t1 && mrc.rowcol < ncols)
  723. X         { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  724. X   }
  725. X   else { REPORT mrc.IncrDiag(); }                     // not accessed
  726. X}
  727. X
  728. X// routines for upper triangular matrix
  729. X
  730. Xvoid UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  731. X{
  732. X   REPORT
  733. X   int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
  734. X   mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  735. X}
  736. X
  737. X
  738. Xvoid UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  739. X{
  740. X   REPORT
  741. X   mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  742. X   Real* ColCopy;
  743. X   if ( !(mrc.cw*StoreHere) )
  744. X   {
  745. X      REPORT                                              // not accessed
  746. X      ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  747. X      MONITOR_REAL_NEW("Make (UT GetCol)",i,ColCopy) 
  748. X      mrc.store = ColCopy;
  749. X   }
  750. X   else { REPORT ColCopy = mrc.store; }
  751. X   if (+(mrc.cw*LoadOnEntry))
  752. X   {
  753. X      REPORT
  754. X      Real* Mstore = store+mrc.rowcol; int j = ncols;
  755. X      while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  756. X   }
  757. X}
  758. X
  759. Xvoid UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  760. X{
  761. X//  if (mrc.cw*StoreOnExit)
  762. X  {
  763. X     REPORT
  764. X     Real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  765. X     Real* Cstore = mrc.store;
  766. X     while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  767. X  }
  768. X}
  769. X
  770. Xvoid UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  771. X                              // 722
  772. X
  773. X// routines for lower triangular matrix
  774. X
  775. Xvoid LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  776. X{
  777. X   REPORT
  778. X   int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  779. X   mrc.store=store+(row*(row+1))/2;
  780. X}
  781. X
  782. Xvoid LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  783. X{
  784. X   REPORT
  785. X   int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
  786. X   int i=nrows-col; mrc.storage=i; Real* ColCopy;
  787. X   if ( !(mrc.cw*StoreHere) )
  788. X   {
  789. X      REPORT                                            // not accessed
  790. X      ColCopy = new Real [i]; MatrixErrorNoSpace(ColCopy);
  791. X      MONITOR_REAL_NEW("Make (LT GetCol)",i,ColCopy) 
  792. X      mrc.store = ColCopy-col;
  793. X   }
  794. X   else { REPORT ColCopy = mrc.store+col; }
  795. X   if (+(mrc.cw*LoadOnEntry))
  796. X   {
  797. X      REPORT
  798. X      Real* Mstore = store+(col*(col+3))/2;
  799. X      while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  800. X   }
  801. X}
  802. X
  803. Xvoid LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  804. X{
  805. X//  if (mrc.cw*StoreOnExit)
  806. X   {
  807. X      REPORT
  808. X      int col=mrc.rowcol; Real* Cstore = mrc.store+col;
  809. X      Real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  810. X      while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  811. X   }
  812. X}
  813. X
  814. Xvoid LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  815. X                                     //712
  816. X// routines for symmetric matrix
  817. X
  818. Xvoid SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  819. X{
  820. X   REPORT                                                //571
  821. X   mrc.skip=0; int row=mrc.rowcol;
  822. X   if (+(mrc.cw*DirectPart))
  823. X   {
  824. X      REPORT
  825. X      mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  826. X   }
  827. X   else
  828. X   {
  829. X      mrc.cw+=IsACopy; mrc.storage=ncols;
  830. X      Real* RowCopy = new Real [ncols]; MatrixErrorNoSpace(RowCopy);
  831. X      MONITOR_REAL_NEW("Make (SymGetRow)",ncols,RowCopy) 
  832. X      mrc.store = RowCopy;
  833. X      if (+(mrc.cw*LoadOnEntry))
  834. X      {
  835. X     REPORT                                         // 544
  836. X         Real* Mstore = store+(row*(row+1))/2; int i = row;
  837. X         while (i--) *RowCopy++ = *Mstore++;
  838. X         i = ncols-row;
  839. X     while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
  840. X      }
  841. X   }
  842. X}
  843. X
  844. X// need to check this out under StoreHere
  845. X
  846. Xvoid SymmetricMatrix::GetCol(MatrixRowCol& mrc)
  847. X{
  848. X   REPORT
  849. X   mrc.skip=0; int col=mrc.rowcol;
  850. X   if (+(mrc.cw*DirectPart))
  851. X   {
  852. X      REPORT                                         // not accessed
  853. X      mrc.cw-=IsACopy; mrc.storage=col+1; mrc.store=store+(col*(col+1))/2;
  854. X   }
  855. X   else
  856. X   {
  857. X      mrc.cw+=IsACopy; mrc.storage=ncols; Real* ColCopy;
  858. X      if ( !(mrc.cw*StoreHere) )
  859. X      {
  860. X         REPORT                                      // not accessed
  861. X         ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
  862. X         MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy) 
  863. X         mrc.store = ColCopy;
  864. X      }
  865. X      else { REPORT ColCopy = mrc.store; }
  866. X      if (+(mrc.cw*LoadOnEntry))
  867. X      {
  868. X         REPORT
  869. X         Real* Mstore = store+(col*(col+1))/2; int i = col;
  870. X         while (i--) *ColCopy++ = *Mstore++;
  871. X         i = ncols-col;
  872. X     while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  873. X      }
  874. X   }
  875. X}
  876. X
  877. X//void SymmetricMatrix::RestoreRow(int row, Real* Rstore)
  878. X//{
  879. X////   if (cw*IsACopy && cw*StoreOnExit)
  880. X//   {
  881. X//      Real* Mstore = store+(row*(row+1))/2; int i = row+1;
  882. X//      while (i--) *Mstore++ = *Rstore++;
  883. X//   }
  884. X//}
  885. X
  886. X//void SymmetricMatrix::RestoreCol(int col, Real* Cstore)
  887. X//{
  888. X////   if (cw*IsACopy && cw*StoreOnExit)
  889. X//   {
  890. X//      Real* Mstore = store+(col*(col+3))/2;
  891. X//      int i = nrows-col; int j = col;
  892. X//      while (i--) { *Mstore = *Cstore++; Mstore+= ++j; }
  893. X//   }
  894. X//}
  895. X
  896. X// routines for row vector
  897. X
  898. Xvoid RowVector::GetCol(MatrixRowCol& mrc)
  899. X{
  900. X   REPORT 
  901. X   mrc.skip=0; mrc.storage=1;
  902. X   if (mrc.cw >= StoreHere)
  903. X   {
  904. X      if (mrc.cw >= LoadOnEntry) { REPORT *(mrc.store) = *(store+mrc.rowcol); }
  905. X      mrc.cw+=IsACopy;
  906. X   }
  907. X   else  { REPORT mrc.store = store+mrc.rowcol; mrc.cw-=IsACopy; }
  908. X                                                         // not accessed
  909. X}
  910. X
  911. Xvoid RowVector::NextCol(MatrixRowCol& mrc) 
  912. X{
  913. X   REPORT
  914. X   if (+(mrc.cw*StoreHere))
  915. X   {
  916. X      if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  917. X                             // not accessed
  918. X      mrc.rowcol++;
  919. X      if (mrc.rowcol < ncols)
  920. X      {
  921. X     if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.store)=*(store+mrc.rowcol); }
  922. X      }
  923. X      else { REPORT mrc.cw -= StoreOnExit; }
  924. X   }
  925. X   else  { REPORT mrc.rowcol++; mrc.store++; }             // not accessed
  926. X}
  927. X
  928. Xvoid RowVector::RestoreCol(MatrixRowCol& mrc)
  929. X{
  930. X   REPORT                                            // not accessed
  931. X   if (mrc.cw>=IsACopy)  { REPORT *(store+mrc.rowcol)=*(mrc.store); }
  932. X}
  933. X
  934. X
  935. X// routines for band matrices
  936. X
  937. Xvoid BandMatrix::GetRow(MatrixRowCol& mrc)
  938. X{
  939. X   REPORT
  940. X   mrc.cw -= IsACopy; int r = mrc.rowcol; int w = lower+1+upper;
  941. X   int s = r-lower; mrc.store = store+(r*w-s); if (s<0) { w += s; s = 0; }
  942. X   mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
  943. X}
  944. X
  945. X// make special versions of this for upper and lower band matrices
  946. X
  947. Xvoid BandMatrix::NextRow(MatrixRowCol& mrc)
  948. X{
  949. X   REPORT
  950. X   int r = ++mrc.rowcol; mrc.store += lower+upper;
  951. X   if (r<=lower) mrc.storage++; else mrc.skip++;
  952. X   if (r>=ncols-upper) mrc.storage--;
  953. X}
  954. X
  955. Xvoid BandMatrix::GetCol(MatrixRowCol& mrc)
  956. X{
  957. X   REPORT
  958. X   mrc.cw += IsACopy; int c = mrc.rowcol; int n = lower+upper; int w = n+1;
  959. X   int b; int s = c-upper; Real* ColCopy;
  960. X   if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
  961. X   mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
  962. X   if ( !(mrc.cw*StoreHere) )
  963. X   {
  964. X      REPORT
  965. X      ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  966. X      MONITOR_REAL_NEW("Make (BMGetCol)",w,ColCopy)
  967. X      mrc.store = ColCopy-mrc.skip;
  968. X   }
  969. X   else { REPORT ColCopy = mrc.store+mrc.skip; }
  970. X   if (+(mrc.cw*LoadOnEntry))
  971. X   {
  972. X      REPORT
  973. X      Real* Mstore = store+b;
  974. X      while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
  975. X   }
  976. X}
  977. X
  978. Xvoid BandMatrix::RestoreCol(MatrixRowCol& mrc)
  979. X{
  980. X//  if (mrc.cw*StoreOnExit)
  981. X   REPORT
  982. X   int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
  983. X   Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
  984. X   Real* Cstore = mrc.store+mrc.skip; int w = mrc.storage;
  985. X   while (w--) { *Mstore = *Cstore++; Mstore += n; }
  986. X}
  987. X
  988. X// routines for symmetric band matrix
  989. X
  990. Xvoid SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
  991. X{
  992. X   REPORT
  993. X   int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
  994. X   if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  995. X
  996. X   if (+(mrc.cw*DirectPart))
  997. X      { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  998. X   else
  999. X   {
  1000. X      mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  1001. X      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
  1002. X      Real* RowCopy = new Real [w]; MatrixErrorNoSpace(RowCopy);
  1003. X      MONITOR_REAL_NEW("Make (SmBGetRow)",w,RowCopy) 
  1004. X      mrc.store = RowCopy-mrc.skip;
  1005. X      if (+(mrc.cw*LoadOnEntry))
  1006. X      {
  1007. X     REPORT
  1008. X         Real* Mstore = store+o;
  1009. X         while (w1--) *RowCopy++ = *Mstore++;   Mstore--;
  1010. X         while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
  1011. X      }
  1012. X   }
  1013. X}
  1014. X
  1015. X// need to check this out under StoreHere
  1016. X
  1017. Xvoid SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
  1018. X{
  1019. X   REPORT
  1020. X   int c=mrc.rowcol; int s = c-lower; int w1 = lower+1; int o = c*w1;
  1021. X   if (s<0) { w1 += s; o -= s; s = 0; }  mrc.skip = s;
  1022. X
  1023. X   if (+(mrc.cw*DirectPart))
  1024. X      { REPORT  mrc.cw -= IsACopy; mrc.store = store+o-s; mrc.storage = w1; }
  1025. X   else
  1026. X   {
  1027. X      mrc.cw += IsACopy; int w = w1+lower; s += w-ncols;
  1028. X      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1; Real* ColCopy;
  1029. X      if ( !(mrc.cw*StoreHere) )
  1030. X      {
  1031. X         ColCopy = new Real [w]; MatrixErrorNoSpace(ColCopy);
  1032. X         MONITOR_REAL_NEW("Make (SmBGetCol)",w,ColCopy) 
  1033. X         mrc.store = ColCopy-mrc.skip;
  1034. X      }
  1035. X      else { REPORT ColCopy = mrc.store+mrc.skip; }
  1036. X      if (+(mrc.cw*LoadOnEntry))
  1037. X      {
  1038. X     REPORT
  1039. X         Real* Mstore = store+o;
  1040. X         while (w1--) *ColCopy++ = *Mstore++;   Mstore--;
  1041. X         while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
  1042. X      }
  1043. X   }
  1044. X}
  1045. X
  1046. END_OF_FILE
  1047. if test 14262 -ne `wc -c <'newmat3.cxx'`; then
  1048.     echo shar: \"'newmat3.cxx'\" unpacked with wrong size!
  1049. fi
  1050. # end of 'newmat3.cxx'
  1051. fi
  1052. if test -f 'newmat4.cxx' -a "${1}" != "-c" ; then 
  1053.   echo shar: Will not clobber existing file \"'newmat4.cxx'\"
  1054. else
  1055. echo shar: Extracting \"'newmat4.cxx'\" \(16990 characters\)
  1056. sed "s/^X//" >'newmat4.cxx' <<'END_OF_FILE'
  1057. X//$$ newmat4.cxx       Constructors, ReDimension, basic utilities
  1058. X
  1059. X// Copyright (C) 1991,2,3: R B Davies
  1060. X
  1061. X#include "include.h"
  1062. X
  1063. X#include "newmat.h"
  1064. X#include "newmatrc.h"
  1065. X
  1066. X//#define REPORT { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  1067. X
  1068. X#define REPORT {}
  1069. X
  1070. X//#define REPORT1 { static ExeCounter ExeCount(__LINE__,4); ExeCount++; }
  1071. X
  1072. X// REPORT1 constructors only - doesn't work in turbo and Borland C++
  1073. X
  1074. X#define REPORT1 {}
  1075. X
  1076. X
  1077. X/*************************** general utilities *************************/
  1078. X
  1079. Xstatic int tristore(int n)                      // els in triangular matrix
  1080. X{ return (n*(n+1))/2; }
  1081. X
  1082. X
  1083. X/****************************** constructors ***************************/
  1084. X
  1085. XGeneralMatrix::GeneralMatrix()
  1086. X{ store=0; storage=0; nrows=0; ncols=0; tag=-1; }
  1087. X
  1088. XGeneralMatrix::GeneralMatrix(ArrayLengthSpecifier s)
  1089. X{
  1090. X   REPORT1
  1091. X   storage=s.Value(); tag=-1;
  1092. X   if (storage)
  1093. X   {
  1094. X      store = new Real [storage]; MatrixErrorNoSpace(store);
  1095. X      MONITOR_REAL_NEW("Make (GenMatrix)",storage,store)
  1096. X   }
  1097. X   else store = 0;
  1098. X}
  1099. X
  1100. XMatrix::Matrix(int m, int n) : GeneralMatrix(m*n)
  1101. X{ REPORT1 nrows=m; ncols=n; }
  1102. X
  1103. XSymmetricMatrix::SymmetricMatrix(ArrayLengthSpecifier n)
  1104. X   : GeneralMatrix(tristore(n.Value()))
  1105. X{ REPORT1 nrows=n.Value(); ncols=n.Value(); }
  1106. X
  1107. XUpperTriangularMatrix::UpperTriangularMatrix(ArrayLengthSpecifier n)
  1108. X   : GeneralMatrix(tristore(n.Value()))
  1109. X{ REPORT1 nrows=n.Value(); ncols=n.Value(); }
  1110. X
  1111. XLowerTriangularMatrix::LowerTriangularMatrix(ArrayLengthSpecifier n)
  1112. X   : GeneralMatrix(tristore(n.Value()))
  1113. X{ REPORT1 nrows=n.Value(); ncols=n.Value(); }
  1114. X
  1115. XDiagonalMatrix::DiagonalMatrix(ArrayLengthSpecifier m) : GeneralMatrix(m)
  1116. X{ REPORT1 nrows=m.Value(); ncols=m.Value(); }
  1117. X
  1118. XMatrix::Matrix(const BaseMatrix& M)
  1119. X{
  1120. X   REPORT1 // CheckConversion(M);
  1121. X   MatrixConversionCheck mcc;
  1122. X   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Rt);
  1123. X   GetMatrix(gmx);
  1124. X}
  1125. X
  1126. XRowVector::RowVector(const BaseMatrix& M) : Matrix(M)
  1127. X{
  1128. X   if (nrows!=1)
  1129. X   {
  1130. X      Tracer tr("RowVector");
  1131. X      Throw(VectorException(*this));
  1132. X   }
  1133. X}
  1134. X
  1135. XColumnVector::ColumnVector(const BaseMatrix& M) : Matrix(M)
  1136. X{
  1137. X   if (ncols!=1)
  1138. X   {
  1139. X      Tracer tr("ColumnVector");
  1140. X      Throw(VectorException(*this));
  1141. X   }
  1142. X}
  1143. X
  1144. XSymmetricMatrix::SymmetricMatrix(const BaseMatrix& M)
  1145. X{
  1146. X   REPORT1  // CheckConversion(M);
  1147. X   MatrixConversionCheck mcc;
  1148. X   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Sm);
  1149. X   GetMatrix(gmx);
  1150. X}
  1151. X
  1152. XUpperTriangularMatrix::UpperTriangularMatrix(const BaseMatrix& M)
  1153. X{
  1154. X   REPORT1 // CheckConversion(M);
  1155. X   MatrixConversionCheck mcc;
  1156. X   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::UT);
  1157. X   GetMatrix(gmx);
  1158. X}
  1159. X
  1160. XLowerTriangularMatrix::LowerTriangularMatrix(const BaseMatrix& M)
  1161. X{
  1162. X   REPORT1 // CheckConversion(M);
  1163. X   MatrixConversionCheck mcc;
  1164. X   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::LT);
  1165. X   GetMatrix(gmx);
  1166. X}
  1167. X
  1168. XDiagonalMatrix::DiagonalMatrix(const BaseMatrix& M)
  1169. X{
  1170. X   REPORT1 //CheckConversion(M);
  1171. X   MatrixConversionCheck mcc;
  1172. X   GeneralMatrix* gmx=((BaseMatrix&)M).Evaluate(MatrixType::Dg);
  1173. X   GetMatrix(gmx);
  1174. X}
  1175. X
  1176. XGeneralMatrix::~GeneralMatrix()
  1177. X{
  1178. X   if (store)
  1179. X   {
  1180. X      MONITOR_REAL_DELETE("Free (GenMatrix)",storage,store)
  1181. X#ifdef Version21
  1182. X      delete [] store;
  1183. X#else
  1184. X      delete [storage] store;
  1185. X#endif
  1186. X   }
  1187. X}
  1188. X
  1189. XCroutMatrix::CroutMatrix(const BaseMatrix& m)
  1190. X{
  1191. X   REPORT1
  1192. X   Tracer tr("CroutMatrix");
  1193. X   GeneralMatrix* gm = ((BaseMatrix&)m).Evaluate(MatrixType::Rt);
  1194. X   GetMatrix(gm);
  1195. X   if (nrows!=ncols) Throw(NotSquareException(*this));
  1196. X   d=TRUE; sing=FALSE;
  1197. X   indx=new int [nrows]; MatrixErrorNoSpace(indx);
  1198. X   MONITOR_INT_NEW("Index (CroutMat)",nrows,indx)
  1199. X   ludcmp();
  1200. X}
  1201. X
  1202. XCroutMatrix::~CroutMatrix()
  1203. X{
  1204. X   MONITOR_INT_DELETE("Index (CroutMat)",nrows,indx)
  1205. X#ifdef Version21
  1206. X   delete [] indx;
  1207. X#else
  1208. X   delete [nrows] indx;
  1209. X#endif
  1210. X}
  1211. X
  1212. X//ReturnMatrixX::ReturnMatrixX(GeneralMatrix& gmx)
  1213. X//{
  1214. X//   REPORT1
  1215. X//   gm = gmx.Image(); gm->ReleaseAndDelete();
  1216. X//}
  1217. X
  1218. X#ifndef TEMPS_DESTROYED_QUICKLY
  1219. X
  1220. XGeneralMatrix::operator ReturnMatrixX() const
  1221. X{
  1222. X   REPORT
  1223. X   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete(); 
  1224. X   return ReturnMatrixX(gm);
  1225. X}
  1226. X
  1227. X#else
  1228. X
  1229. XGeneralMatrix::operator ReturnMatrixX&() const
  1230. X{
  1231. X   REPORT
  1232. X   GeneralMatrix* gm = Image(); gm->ReleaseAndDelete();
  1233. X   ReturnMatrixX* x = new ReturnMatrixX(gm);
  1234. X   MatrixErrorNoSpace(x); return *x;
  1235. X}
  1236. X
  1237. X#endif
  1238. X
  1239. X/**************************** ReDimension matrices ***************************/
  1240. X
  1241. Xvoid GeneralMatrix::ReDimension(int nr, int nc, int s)
  1242. X{
  1243. X   REPORT 
  1244. X   if (store)
  1245. X   {
  1246. X      MONITOR_REAL_DELETE("Free (ReDimensi)",storage,store)
  1247. X#ifdef Version21
  1248. X      delete [] store;
  1249. X#else
  1250. X      delete [storage] store;
  1251. X#endif
  1252. X   }
  1253. X   storage=s; nrows=nr; ncols=nc; tag=-1;
  1254. X   if (s)
  1255. X   {
  1256. X      store = new Real [storage]; MatrixErrorNoSpace(store);
  1257. X      MONITOR_REAL_NEW("Make (ReDimensi)",storage,store)
  1258. X   }
  1259. X   else store = 0;
  1260. X}
  1261. X
  1262. Xvoid Matrix::ReDimension(int nr, int nc)
  1263. X{ REPORT GeneralMatrix::ReDimension(nr,nc,nr*nc); }
  1264. X
  1265. Xvoid SymmetricMatrix::ReDimension(int nr)
  1266. X{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  1267. X
  1268. Xvoid UpperTriangularMatrix::ReDimension(int nr)
  1269. X{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  1270. X
  1271. Xvoid LowerTriangularMatrix::ReDimension(int nr)
  1272. X{ REPORT GeneralMatrix::ReDimension(nr,nr,tristore(nr)); }
  1273. X
  1274. Xvoid DiagonalMatrix::ReDimension(int nr)
  1275. X{ REPORT GeneralMatrix::ReDimension(nr,nr,nr); }
  1276. X
  1277. Xvoid RowVector::ReDimension(int nc)
  1278. X{ REPORT GeneralMatrix::ReDimension(1,nc,nc); }
  1279. X
  1280. Xvoid ColumnVector::ReDimension(int nr)
  1281. X{ REPORT GeneralMatrix::ReDimension(nr,1,nr); }
  1282. X
  1283. X
  1284. X/********************* manipulate types, storage **************************/
  1285. X
  1286. Xint GeneralMatrix::search(const BaseMatrix* s) const
  1287. X{ REPORT return (s==this) ? 1 : 0; }
  1288. X
  1289. Xint MultipliedMatrix::search(const BaseMatrix* s) const
  1290. X{ REPORT return bm1->search(s) + bm2->search(s); }
  1291. X
  1292. Xint ShiftedMatrix::search(const BaseMatrix* s) const
  1293. X{ REPORT return bm->search(s); }
  1294. X
  1295. Xint NegatedMatrix::search(const BaseMatrix* s) const
  1296. X{ REPORT return bm->search(s); }
  1297. X
  1298. Xint ConstMatrix::search(const BaseMatrix* s) const
  1299. X{ REPORT return (s==cgm) ? 1 : 0; }
  1300. X
  1301. Xint ReturnMatrixX::search(const BaseMatrix* s) const
  1302. X{ REPORT return (s==gm) ? 1 : 0; }
  1303. X
  1304. XMatrixType Matrix::Type() const { return MatrixType::Rt; }
  1305. XMatrixType SymmetricMatrix::Type() const { return MatrixType::Sm; }
  1306. XMatrixType UpperTriangularMatrix::Type() const { return MatrixType::UT; }
  1307. XMatrixType LowerTriangularMatrix::Type() const { return MatrixType::LT; }
  1308. XMatrixType DiagonalMatrix::Type() const { return MatrixType::Dg; }
  1309. XMatrixType RowVector::Type() const { return MatrixType::RV; }
  1310. XMatrixType ColumnVector::Type() const { return MatrixType::CV; }
  1311. XMatrixType CroutMatrix::Type() const { return MatrixType::Ct; }
  1312. XMatrixType BandMatrix::Type() const { return MatrixType::BM; }
  1313. XMatrixType UpperBandMatrix::Type() const { return MatrixType::UB; }
  1314. XMatrixType LowerBandMatrix::Type() const { return MatrixType::LB; }
  1315. XMatrixType SymmetricBandMatrix::Type() const { return MatrixType::SB; }
  1316. X
  1317. XMatrixBandWidth BaseMatrix::BandWidth() const { return -1; }
  1318. XMatrixBandWidth DiagonalMatrix::BandWidth() const { return 0; }
  1319. X
  1320. XMatrixBandWidth BandMatrix::BandWidth() const
  1321. X   { return MatrixBandWidth(lower,upper); }
  1322. X
  1323. XMatrixBandWidth AddedMatrix::BandWidth() const
  1324. X{ return gm1->BandWidth() + gm2->BandWidth(); }
  1325. X
  1326. XMatrixBandWidth MultipliedMatrix::BandWidth() const
  1327. X{ return gm1->BandWidth() * gm2->BandWidth(); }
  1328. X
  1329. XMatrixBandWidth SolvedMatrix::BandWidth() const { return -1; }
  1330. XMatrixBandWidth ScaledMatrix::BandWidth() const { return gm->BandWidth(); }
  1331. XMatrixBandWidth NegatedMatrix::BandWidth() const { return gm->BandWidth(); }
  1332. X
  1333. XMatrixBandWidth TransposedMatrix::BandWidth() const
  1334. X{ return gm->BandWidth().t(); }
  1335. X
  1336. XMatrixBandWidth InvertedMatrix::BandWidth() const { return -1; }
  1337. XMatrixBandWidth RowedMatrix::BandWidth() const { return -1; }
  1338. XMatrixBandWidth ColedMatrix::BandWidth() const { return -1; }
  1339. XMatrixBandWidth DiagedMatrix::BandWidth() const { return 0; }
  1340. XMatrixBandWidth MatedMatrix::BandWidth() const { return -1; }
  1341. XMatrixBandWidth ConstMatrix::BandWidth() const { return cgm->BandWidth(); }
  1342. XMatrixBandWidth ReturnMatrixX::BandWidth() const { return gm->BandWidth(); }
  1343. X
  1344. XMatrixBandWidth GetSubMatrix::BandWidth() const
  1345. X{
  1346. X
  1347. X   if (row_skip==col_skip && row_number==col_number) return gm->BandWidth();
  1348. X   else return MatrixBandWidth(-1);
  1349. X}
  1350. X
  1351. X/************************ the memory managment tools **********************/
  1352. X
  1353. X//  Rules regarding tDelete, reuse, GetStore
  1354. X//    All matrices processed during expression evaluation must be subject
  1355. X//    to exactly one of reuse(), tDelete(), GetStore() or BorrowStore().
  1356. X//    If reuse returns TRUE the matrix must be reused.
  1357. X//    GetMatrix(gm) always calls gm->GetStore()
  1358. X//    gm->Evaluate obeys rules
  1359. X//    bm->Evaluate obeys rules for matrices in bm structure
  1360. X
  1361. Xvoid GeneralMatrix::tDelete()
  1362. X{
  1363. X   if (tag<0)
  1364. X   {
  1365. X      if (tag<-1) { REPORT store=0; delete this; return; }  // borrowed
  1366. X      else { REPORT return; }
  1367. X   }
  1368. X   if (tag==1)
  1369. X   {
  1370. X      REPORT  MONITOR_REAL_DELETE("Free   (tDelete)",storage,store)
  1371. X#ifdef Version21
  1372. X      if (store) delete [] store;
  1373. X#else
  1374. X      if (store) delete [storage] store;
  1375. X#endif
  1376. X      store=0; tag=-1; return;
  1377. X   }
  1378. X   if (tag==0) { REPORT delete this; return; }
  1379. X   REPORT tag--; return;
  1380. X}
  1381. X
  1382. Xstatic void BlockCopy(int n, Real* from, Real* to)
  1383. X{
  1384. X   REPORT
  1385. X   int i = (n >> 3);
  1386. X   while (i--)
  1387. X   {
  1388. X      *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
  1389. X      *to++ = *from++; *to++ = *from++; *to++ = *from++; *to++ = *from++;
  1390. X   }
  1391. X   i = n & 7; while (i--) *to++ = *from++;
  1392. X}
  1393. X
  1394. XBoolean GeneralMatrix::reuse()
  1395. X{
  1396. X   if (tag<-1)
  1397. X   {
  1398. X      REPORT
  1399. X      Real* s = new Real [storage]; MatrixErrorNoSpace(s);
  1400. X      MONITOR_REAL_NEW("Make     (reuse)",storage,s)
  1401. X      BlockCopy(storage, store, s); store=s; tag=0; return TRUE;
  1402. X   }
  1403. X   if (tag<0) { REPORT return FALSE; }
  1404. X   if (tag<=1)  { REPORT return TRUE; }
  1405. X   REPORT tag--; return FALSE;
  1406. X}
  1407. X
  1408. XReal* GeneralMatrix::GetStore()
  1409. X{
  1410. X   if (tag<0 || tag>1)
  1411. X   {
  1412. X      Real* s;
  1413. X      if (storage)
  1414. X      {
  1415. X         s = new Real [storage]; MatrixErrorNoSpace(s);
  1416. X         MONITOR_REAL_NEW("Make  (GetStore)",storage,s)
  1417. X         BlockCopy(storage, store, s);
  1418. X      }
  1419. X      else s = 0;
  1420. X      if (tag>1) { REPORT tag--; }
  1421. X      else if (tag < -1) { REPORT store=0; delete this; } // borrowed store
  1422. X      else { REPORT }
  1423. X      return s;
  1424. X   }
  1425. X   Real* s=store; store=0;
  1426. X   if (tag==0) { REPORT delete this; }
  1427. X   else { REPORT tag=-1; }
  1428. X   return s;
  1429. X}
  1430. X
  1431. X/*
  1432. X#ifndef __ZTC__
  1433. Xvoid GeneralMatrix::GetMatrixC(const GeneralMatrix* gmx)
  1434. X{
  1435. X   REPORT tag=-1;
  1436. X   nrows=gmx->nrows; ncols=gmx->ncols; storage=gmx->storage;
  1437. X   SetParameters(gmx); 
  1438. X   store = new Real [storage]; MatrixErrorNoSpace(store);
  1439. X   MONITOR_REAL_NEW("Make (GetMatrix)",storage,store)
  1440. X   BlockCopy(storage, gmx->store, store);
  1441. X}
  1442. X#endif
  1443. X*/
  1444. X
  1445. Xvoid GeneralMatrix::GetMatrix(const GeneralMatrix* gmx)
  1446. X{
  1447. X   REPORT  tag=-1; nrows=gmx->Nrows(); ncols=gmx->Ncols();
  1448. X   storage=gmx->storage; SetParameters(gmx);
  1449. X   store=((GeneralMatrix*)gmx)->GetStore();
  1450. X}
  1451. X
  1452. XGeneralMatrix* GeneralMatrix::BorrowStore(GeneralMatrix* gmx, MatrixType mt)
  1453. X// Copy storage of *this to storage of *gmx. Then convert to type mt.
  1454. X// If mt == 0 just let *gmx point to storage of *this if tag==-1.
  1455. X{
  1456. X   if (!mt)
  1457. X   {
  1458. X      if (tag == -1) { REPORT gmx->tag = -2; gmx->store = store; }
  1459. X      else { REPORT gmx->tag = 0; gmx->store = GetStore(); }
  1460. X   }
  1461. X   else if (Compare(gmx->Type(),mt))
  1462. X   { REPORT  gmx->tag = 0; gmx->store = GetStore(); }
  1463. X   else
  1464. X   {
  1465. X      REPORT gmx->tag = -2; gmx->store = store;
  1466. X      gmx = gmx->Evaluate(mt); gmx->tag = 0; tDelete();
  1467. X   }
  1468. X
  1469. X   return gmx;
  1470. X}
  1471. X
  1472. Xvoid GeneralMatrix::Eq(const BaseMatrix& X, MatrixType mt)
  1473. X// Count number of references to this in X.
  1474. X// If zero delete storage in X;
  1475. X// otherwise tag X to show when storage can be deleted
  1476. X// evaluate X and copy to current object
  1477. X{
  1478. X   int counter=X.search(this);
  1479. X   if (counter==0)
  1480. X   {
  1481. X      REPORT
  1482. X      if (store)
  1483. X      {
  1484. X         MONITOR_REAL_DELETE("Free (operator=)",storage,store)
  1485. X#ifdef Version21
  1486. X         REPORT delete [] store; storage=0;
  1487. X#else
  1488. X         REPORT delete [storage] store; storage=0;
  1489. X#endif
  1490. X      }
  1491. X   }
  1492. X   else { REPORT Release(counter); }
  1493. X   GeneralMatrix* gmx = ((BaseMatrix&)X).Evaluate(mt);
  1494. X   if (gmx!=this) { REPORT GetMatrix(gmx); }
  1495. X   else { REPORT }
  1496. X   Protect();
  1497. X}
  1498. X
  1499. Xvoid GeneralMatrix::Inject(const GeneralMatrix& X)
  1500. X// copy stored values of X; otherwise leave els of *this unchanged
  1501. X{
  1502. X   REPORT
  1503. X   Tracer tr("Inject");
  1504. X   if (nrows != X.nrows || ncols != X.ncols)
  1505. X      Throw(IncompatibleDimensionsException());
  1506. X   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
  1507. X   MatrixRow mrx(this, LoadOnEntry+StoreOnExit+DirectPart);
  1508. X   int i=nrows;
  1509. X   while (i--) { mrx.Inject(mr); mrx.Next(); mr.Next(); }
  1510. X}  
  1511. X
  1512. X/*************** checking for data loss during conversion *******************/
  1513. X
  1514. X//void GeneralMatrix::CheckConversion(const BaseMatrix& M)
  1515. X//{
  1516. X//   if (!(this->Type() >= M.Type()))
  1517. X//      Throw(ProgramException("Illegal Conversion"));
  1518. X//}
  1519. X
  1520. XBoolean MatrixConversionCheck::DoCheck = FALSE;
  1521. X
  1522. Xvoid MatrixConversionCheck::DataLoss()
  1523. X   { if (DoCheck) Throw(ProgramException("Illegal Conversion")); }
  1524. X
  1525. XBoolean Compare(const MatrixType& source, MatrixType& destination)
  1526. X{
  1527. X   if (!destination) { destination=source; return TRUE; }
  1528. X   if (destination==source) return TRUE;
  1529. X   if (MatrixConversionCheck::IsOn() && !(destination>=source))
  1530. X      Throw(ProgramException("Illegal Conversion"));
  1531. X   return FALSE;
  1532. X}
  1533. X
  1534. X/*************** Make a copy of a matrix on the heap *********************/
  1535. X
  1536. XGeneralMatrix* Matrix::Image() const
  1537. X{
  1538. X   REPORT
  1539. X   GeneralMatrix* gm = new Matrix(*this); MatrixErrorNoSpace(gm);
  1540. X   return gm;
  1541. X}
  1542. X
  1543. XGeneralMatrix* SymmetricMatrix::Image() const
  1544. X{
  1545. X   REPORT
  1546. X   GeneralMatrix* gm = new SymmetricMatrix(*this); MatrixErrorNoSpace(gm);
  1547. X   return gm;
  1548. X}
  1549. X
  1550. XGeneralMatrix* UpperTriangularMatrix::Image() const
  1551. X{
  1552. X   REPORT
  1553. X   GeneralMatrix* gm = new UpperTriangularMatrix(*this);
  1554. X   MatrixErrorNoSpace(gm); return gm;
  1555. X}
  1556. X
  1557. XGeneralMatrix* LowerTriangularMatrix::Image() const
  1558. X{
  1559. X   REPORT
  1560. X   GeneralMatrix* gm = new LowerTriangularMatrix(*this);
  1561. X   MatrixErrorNoSpace(gm); return gm;
  1562. X}
  1563. X
  1564. XGeneralMatrix* DiagonalMatrix::Image() const
  1565. X{
  1566. X   REPORT
  1567. X   GeneralMatrix* gm = new DiagonalMatrix(*this); MatrixErrorNoSpace(gm);
  1568. X   return gm;
  1569. X}
  1570. X
  1571. XGeneralMatrix* RowVector::Image() const
  1572. X{
  1573. X   REPORT
  1574. X   GeneralMatrix* gm = new RowVector(*this); MatrixErrorNoSpace(gm);
  1575. X   return gm;
  1576. X}
  1577. X
  1578. XGeneralMatrix* ColumnVector::Image() const
  1579. X{
  1580. X   REPORT
  1581. X   GeneralMatrix* gm = new ColumnVector(*this); MatrixErrorNoSpace(gm);
  1582. X   return gm;
  1583. X}
  1584. X
  1585. XGeneralMatrix* BandMatrix::Image() const
  1586. X{
  1587. X   REPORT
  1588. X   GeneralMatrix* gm = new BandMatrix(*this); MatrixErrorNoSpace(gm);
  1589. X   return gm;
  1590. X}
  1591. X
  1592. XGeneralMatrix* UpperBandMatrix::Image() const
  1593. X{
  1594. X   REPORT
  1595. X   GeneralMatrix* gm = new UpperBandMatrix(*this); MatrixErrorNoSpace(gm);
  1596. X   return gm;
  1597. X}
  1598. X
  1599. XGeneralMatrix* LowerBandMatrix::Image() const
  1600. X{
  1601. X   REPORT
  1602. X   GeneralMatrix* gm = new LowerBandMatrix(*this); MatrixErrorNoSpace(gm);
  1603. X   return gm;
  1604. X}
  1605. X
  1606. XGeneralMatrix* SymmetricBandMatrix::Image() const
  1607. X{
  1608. X   REPORT
  1609. X   GeneralMatrix* gm = new SymmetricBandMatrix(*this); MatrixErrorNoSpace(gm);
  1610. X   return gm;
  1611. X}
  1612. X
  1613. XGeneralMatrix* nricMatrix::Image() const
  1614. X{
  1615. X   REPORT
  1616. X   GeneralMatrix* gm = new nricMatrix(*this); MatrixErrorNoSpace(gm);
  1617. X   return gm;
  1618. X}
  1619. X
  1620. XGeneralMatrix* GeneralMatrix::Image() const
  1621. X{
  1622. X   REPORT
  1623. X   Throw(InternalException("Cannot apply Image to this matrix type"));
  1624. X   return 0;
  1625. X}
  1626. X
  1627. X
  1628. X/************************* nricMatrix routines *****************************/
  1629. X
  1630. Xvoid nricMatrix::MakeRowPointer()
  1631. X{
  1632. X   row_pointer = new Real* [nrows]; MatrixErrorNoSpace(row_pointer);
  1633. X   Real* s = Store() - 1; int i = nrows; Real** rp = row_pointer;
  1634. X   while (i--) { *rp++ = s; s+=ncols; }
  1635. X}
  1636. X
  1637. Xvoid nricMatrix::DeleteRowPointer()
  1638. X#ifdef Version21
  1639. X{ if (nrows) delete [] row_pointer; }
  1640. X#else
  1641. X{ if (nrows) delete [nrows] row_pointer; }
  1642. X#endif
  1643. X
  1644. Xvoid GeneralMatrix::CheckStore() const
  1645. X{
  1646. X   if (!store) 
  1647. X      Throw(ProgramException("NRIC accessing matrix with unset dimensions"));
  1648. X}
  1649. X
  1650. X
  1651. X/***************************** CleanUp routines *****************************/
  1652. X
  1653. Xvoid GeneralMatrix::CleanUp()
  1654. X{
  1655. X   // set matrix dimensions to zero, delete storage
  1656. X   REPORT
  1657. X   if (store && storage)
  1658. X   {
  1659. X      MONITOR_REAL_DELETE("Free (CleanUp)    ",storage,store)
  1660. X#ifdef Version21
  1661. X         REPORT delete [] store;
  1662. X#else
  1663. X         REPORT delete [storage] store;
  1664. X#endif
  1665. X   }
  1666. X   store=0; storage=0; nrows=0; ncols=0;
  1667. X}
  1668. X
  1669. Xvoid nricMatrix::CleanUp()
  1670. X{ DeleteRowPointer(); GeneralMatrix::CleanUp(); }
  1671. X
  1672. Xvoid RowVector::CleanUp()
  1673. X{ GeneralMatrix::CleanUp(); nrows=1; }
  1674. X
  1675. Xvoid ColumnVector::CleanUp()
  1676. X{ GeneralMatrix::CleanUp(); ncols=1; }
  1677. X
  1678. Xvoid CroutMatrix::CleanUp()
  1679. X{ 
  1680. X#ifdef Version21
  1681. X   if (nrows) delete [] indx;
  1682. X#else
  1683. X   if (nrows) delete [nrows] indx;
  1684. X#endif
  1685. X   GeneralMatrix::CleanUp();
  1686. X}
  1687. X
  1688. Xvoid BandLUMatrix::CleanUp()
  1689. X{ 
  1690. X#ifdef Version21
  1691. X   if (nrows) delete [] indx;
  1692. X   if (storage2) delete [] store2;
  1693. X#else
  1694. X   if (nrows) delete [nrows] indx;
  1695. X   if (storage2) delete [storage2] store2;
  1696. X#endif
  1697. X   GeneralMatrix::CleanUp();
  1698. X}
  1699. X
  1700. X
  1701. END_OF_FILE
  1702. if test 16990 -ne `wc -c <'newmat4.cxx'`; then
  1703.     echo shar: \"'newmat4.cxx'\" unpacked with wrong size!
  1704. fi
  1705. # end of 'newmat4.cxx'
  1706. fi
  1707. if test -f 'newmatex.cxx' -a "${1}" != "-c" ; then 
  1708.   echo shar: Will not clobber existing file \"'newmatex.cxx'\"
  1709. else
  1710. echo shar: Extracting \"'newmatex.cxx'\" \(9162 characters\)
  1711. sed "s/^X//" >'newmatex.cxx' <<'END_OF_FILE'
  1712. X//$$ newmatex.cxx                    Exception handler
  1713. X
  1714. X// Copyright (C) 1992: R B Davies
  1715. X
  1716. X#define WANT_STREAM                  // include.h will get stream fns
  1717. X
  1718. X#include "include.h"                 // include standard files
  1719. X#include "newmat.h"
  1720. X
  1721. X
  1722. X// action = -1    print message and exit(1)
  1723. X//           0    no message if handler available
  1724. X//           1    print message and use handler
  1725. X
  1726. X
  1727. X
  1728. Xint SpaceException::action = 1;
  1729. Xint DataException::action = 1;
  1730. Xint ConvergenceException::action = 1;
  1731. Xint ProgramException::action = 1;
  1732. Xint InternalException::action = 1;
  1733. X
  1734. X
  1735. Xstatic inline iabs(int i) { return i >= 0 ? i : -i; }
  1736. X
  1737. X
  1738. XMatrixDetails::MatrixDetails(const GeneralMatrix& A)
  1739. X   : type(A.Type()), nrows(A.Nrows()), ncols(A.Ncols())
  1740. X{ MatrixBandWidth bw = A.BandWidth(); ubw = bw.upper; lbw = bw.lower; }
  1741. X
  1742. Xvoid MatrixDetails::PrintOut()
  1743. X{
  1744. X   cout << "MatrixType = " << (char*)type;
  1745. X   cout << "  # Rows = " << nrows;
  1746. X   cout << "; # Cols = " << ncols;
  1747. X   if (lbw >=0) cout << "; lower BW = " << lbw;
  1748. X   if (ubw >=0) cout << "; upper BW = " << ubw;
  1749. X   cout << "\n";
  1750. X}
  1751. X
  1752. X
  1753. X
  1754. XSpaceException::SpaceException() : Exception(iabs(action))
  1755. X{
  1756. X   if (action) cout << "Out of space on heap\n";
  1757. X   if (action < 0) exit(1);
  1758. X}
  1759. X
  1760. XMatrixException::MatrixException(int action) : Exception(iabs(action))
  1761. X{ if (action) cout << "The exception is from newmat.\n"; }
  1762. X
  1763. XMatrixException::MatrixException(int action, const GeneralMatrix& A)
  1764. X   : Exception(iabs(action))
  1765. X{
  1766. X   if (action)
  1767. X   {
  1768. X      cout << "The exception is from newmat: details of matrix follow:\n";
  1769. X      MatrixDetails(A).PrintOut();
  1770. X   }
  1771. X}
  1772. X
  1773. XMatrixException::MatrixException(int action, const GeneralMatrix& A,
  1774. X   const GeneralMatrix& B) : Exception(iabs(action))
  1775. X{
  1776. X   if (action)
  1777. X   {
  1778. X      cout << "The exception is from newmat: details of matrices follow:\n";
  1779. X      MatrixDetails(A).PrintOut();
  1780. X      MatrixDetails(B).PrintOut();
  1781. X   }
  1782. X}
  1783. X
  1784. XDataException::DataException(const GeneralMatrix& A)
  1785. X   : MatrixException(action, A) {}
  1786. X
  1787. XNPDException::NPDException(const GeneralMatrix& A)
  1788. X   : DataException(A)
  1789. X{
  1790. X   if (action) cout << "The matrix is not positive definite\n\n";
  1791. X   if (action < 0) exit(1);
  1792. X}
  1793. X
  1794. XSingularException::SingularException(const GeneralMatrix& A)
  1795. X   : DataException(A)
  1796. X{
  1797. X   if (action) cout << "The matrix is singular\n\n";
  1798. X   if (action < 0) exit(1);
  1799. X}
  1800. X
  1801. XConvergenceException::ConvergenceException(const GeneralMatrix& A)
  1802. X   : MatrixException(action,A)
  1803. X{
  1804. X   if (action) cout << "Process fails to converge\n\n";
  1805. X   if (action < 0) exit(1);
  1806. X}
  1807. X
  1808. XProgramException::ProgramException(char* c) : MatrixException(action)
  1809. X{
  1810. X   if (action) cout << c << "\n\n";
  1811. X   if (action < 0) exit(1);
  1812. X}
  1813. X
  1814. XProgramException::ProgramException(char* c, const GeneralMatrix& A)
  1815. X   : MatrixException(action,A)
  1816. X{
  1817. X   if (action) cout << c << "\n\n";
  1818. X   if (action < 0) exit(1);
  1819. X}
  1820. X
  1821. XProgramException::ProgramException(char* c, const GeneralMatrix& A,
  1822. X   const GeneralMatrix& B) : MatrixException(action,A,B)
  1823. X{
  1824. X   if (action) cout << c << "\n\n";
  1825. X   if (action < 0) exit(1);
  1826. X}
  1827. X
  1828. XProgramException::ProgramException(const GeneralMatrix& A)
  1829. X   : MatrixException(action, A) {}
  1830. X
  1831. XProgramException::ProgramException() : MatrixException(action) {}
  1832. X
  1833. XVectorException::VectorException() : ProgramException()
  1834. X{
  1835. X   if (action) cout << "Cannot convert matrix to vector\n\n";
  1836. X   if (action < 0) exit(1);
  1837. X}
  1838. X
  1839. XVectorException::VectorException(const GeneralMatrix& A)
  1840. X   : ProgramException(A)
  1841. X{
  1842. X   if (action) cout << "Cannot convert matrix to vector\n\n";
  1843. X   if (action < 0) exit(1);
  1844. X}
  1845. X
  1846. XNotSquareException::NotSquareException(const GeneralMatrix& A)
  1847. X   : ProgramException(A)
  1848. X{
  1849. X   if (action) cout << "Matrix is not square\n\n";
  1850. X   if (action < 0) exit(1);
  1851. X}
  1852. X
  1853. XSubMatrixDimensionException::SubMatrixDimensionException()
  1854. X   : ProgramException()
  1855. X{
  1856. X   if (action) cout << "Incompatible submatrix dimension\n\n";
  1857. X   if (action < 0) exit(1);
  1858. X}
  1859. X
  1860. XIncompatibleDimensionsException::IncompatibleDimensionsException()
  1861. X   : ProgramException()
  1862. X{
  1863. X   if (action) cout << "Incompatible dimensions\n\n";
  1864. X   if (action < 0) exit(1);
  1865. X}
  1866. X
  1867. XNotDefinedException::NotDefinedException(char* op, char* matrix)
  1868. X   : ProgramException()
  1869. X{
  1870. X   if (action)
  1871. X      cout << "Operation " << op << " not defined for " << matrix << "\n\n";
  1872. X   if (action < 0) exit(1);
  1873. X}
  1874. X
  1875. XCannotBuildException::CannotBuildException(char* matrix)
  1876. X   : ProgramException()
  1877. X{
  1878. X   if (action)
  1879. X      cout << "Cannot build matrix type " << matrix << "\n\n";
  1880. X   if (action < 0) exit(1);
  1881. X}
  1882. X
  1883. XIndexException::IndexException(int i, const GeneralMatrix& A)
  1884. X   : ProgramException(A)
  1885. X{
  1886. X   if (action)
  1887. X      { cout << "Index error: requested index = " << i << "\n\n"; }
  1888. X   if (action < 0) exit(1);
  1889. X}
  1890. X
  1891. XIndexException::IndexException(int i, int j, const GeneralMatrix& A)
  1892. X   : ProgramException(A)
  1893. X{
  1894. X   if (action)
  1895. X   {
  1896. X      cout << "Index error: requested indices = " << i << ", " << j << "\n\n";
  1897. X   }
  1898. X   if (action < 0) exit(1);
  1899. X}
  1900. X
  1901. X
  1902. XIndexException::IndexException(int i, const GeneralMatrix& A, Boolean)
  1903. X   : ProgramException(A)
  1904. X{
  1905. X   if (action)
  1906. X      { cout << "Element error: requested index (wrt 0) = " << i << "\n\n"; }
  1907. X   if (action < 0) exit(1);
  1908. X}
  1909. X
  1910. XIndexException::IndexException(int i, int j, const GeneralMatrix& A, Boolean)
  1911. X   : ProgramException(A)
  1912. X{
  1913. X   if (action)
  1914. X   {
  1915. X      cout << "Element error: requested indices (wrt 0) = " 
  1916. X         << i << ", " << j << "\n\n";
  1917. X   }
  1918. X   if (action < 0) exit(1);
  1919. X}
  1920. X
  1921. XInternalException::InternalException(char* c) : MatrixException(action)
  1922. X{
  1923. X   if (action) cout << c << "\n\n";
  1924. X   if (action < 0) exit(1);
  1925. X}
  1926. X
  1927. X
  1928. X
  1929. X
  1930. X/************************* ExeCounter functions *****************************/
  1931. X
  1932. X
  1933. X
  1934. Xint ExeCounter::nreports = 0;
  1935. X
  1936. XExeCounter::ExeCounter(int xl, int xf) : line(xl), fileid(xf), nexe(0) {}
  1937. X
  1938. XExeCounter::~ExeCounter()
  1939. X{
  1940. X   nreports++;
  1941. X   cout << nreports << "  " << fileid << "  " << line << "  " << nexe << "\n";
  1942. X}
  1943. X
  1944. X
  1945. X 
  1946. X/**************************** error handler *******************************/
  1947. X
  1948. Xvoid MatrixErrorNoSpace(void* v) { if (!v) Throw(SpaceException()); }
  1949. X// throw exception if v is null
  1950. X
  1951. X
  1952. X
  1953. X/************************* test type manipulation **************************/
  1954. X
  1955. X
  1956. X
  1957. X
  1958. X// These functions may cause problems for Glockenspiel 2.0c; they are used
  1959. X// only for testing so you can delete them
  1960. X
  1961. X
  1962. Xvoid TestTypeAdd()
  1963. X{
  1964. X   MatrixType list[] = { MatrixType::UT,
  1965. X                         MatrixType::LT,
  1966. X                         MatrixType::Rt,
  1967. X                         MatrixType::Sm,
  1968. X             MatrixType::Dg,
  1969. X                         MatrixType::BM,
  1970. X                         MatrixType::UB,
  1971. X             MatrixType::LB,
  1972. X             MatrixType::SB };
  1973. X
  1974. X   cout << "+     ";
  1975. X   for (int i=0; i<MatrixType::nTypes(); i++) cout << (char*)list[i] << " ";
  1976. X   cout << "\n";
  1977. X   for (i=0; i<MatrixType::nTypes(); i++)
  1978. X   {
  1979. X      cout << (char*)(list[i]) << " ";
  1980. X      for (int j=0; j<MatrixType::nTypes(); j++)
  1981. X     cout << (char*)(list[j]+list[i]) << " ";
  1982. X      cout << "\n";
  1983. X   }
  1984. X   cout << "\n";
  1985. X}
  1986. X
  1987. Xvoid TestTypeMult()
  1988. X{
  1989. X   MatrixType list[] = { MatrixType::UT,
  1990. X                         MatrixType::LT,
  1991. X                         MatrixType::Rt,
  1992. X                         MatrixType::Sm,
  1993. X                         MatrixType::Dg,
  1994. X                         MatrixType::BM,
  1995. X                         MatrixType::UB,
  1996. X                         MatrixType::LB,
  1997. X                         MatrixType::SB };
  1998. X   cout << "*     ";
  1999. X   for (int i=0; i<MatrixType::nTypes(); i++)
  2000. X      cout << (char*)list[i] << " ";
  2001. X   cout << "\n";
  2002. X   for (i=0; i<MatrixType::nTypes(); i++)
  2003. X   {
  2004. X      cout << (char*)list[i] << " ";
  2005. X      for (int j=0; j<MatrixType::nTypes(); j++)
  2006. X     cout << (char*)(list[j]*list[i]) << " ";
  2007. X      cout << "\n";
  2008. X   }
  2009. X   cout << "\n";
  2010. X}
  2011. X
  2012. Xvoid TestTypeOrder()
  2013. X{
  2014. X   MatrixType list[] = { MatrixType::UT,
  2015. X                         MatrixType::LT,
  2016. X                         MatrixType::Rt,
  2017. X                         MatrixType::Sm,
  2018. X                         MatrixType::Dg,
  2019. X                         MatrixType::BM,
  2020. X                         MatrixType::UB,
  2021. X                         MatrixType::LB,
  2022. X                         MatrixType::SB };
  2023. X   cout << ">=    ";
  2024. X   for (int i = 0; i<MatrixType::nTypes(); i++)
  2025. X      cout << (char*)list[i] << " ";
  2026. X   cout << "\n";
  2027. X   for (i=0; i<MatrixType::nTypes(); i++)
  2028. X   {
  2029. X      cout << (char*)list[i] << " ";
  2030. X      for (int j=0; j<MatrixType::nTypes(); j++)
  2031. X     cout << ((list[j]>=list[i]) ? "Yes   " : "No    ");
  2032. X      cout << "\n";
  2033. X   }
  2034. X   cout << "\n";
  2035. X}
  2036. X
  2037. X
  2038. X/************************* miscellanous errors ***************************/
  2039. X
  2040. X
  2041. Xvoid CroutMatrix::GetRow(MatrixRowCol&)
  2042. X   { Throw(NotDefinedException("GetRow","Crout")); }
  2043. Xvoid CroutMatrix::GetCol(MatrixRowCol&)
  2044. X   { Throw(NotDefinedException("GetCol","Crout")); }
  2045. Xvoid CroutMatrix::operator=(const BaseMatrix&)
  2046. X   { Throw(NotDefinedException("=","Crout")); }
  2047. Xvoid BandLUMatrix::GetRow(MatrixRowCol&)
  2048. X   { Throw(NotDefinedException("GetRow","BandLUMatrix")); }
  2049. Xvoid BandLUMatrix::GetCol(MatrixRowCol&)
  2050. X   { Throw(NotDefinedException("GetCol","BandLUMatrix")); }
  2051. Xvoid BandLUMatrix::operator=(const BaseMatrix&)
  2052. X   { Throw(NotDefinedException("=","BandLUMatrix")); }
  2053. X#ifdef TEMPS_DESTROYED_QUICKLY
  2054. X   ReturnMatrixX::ReturnMatrixX(const ReturnMatrixX& tm)
  2055. X     : gm(tm.gm) { Throw(ProgramException("ReturnMatrixX error")); }
  2056. X#endif
  2057. X
  2058. END_OF_FILE
  2059. if test 9162 -ne `wc -c <'newmatex.cxx'`; then
  2060.     echo shar: \"'newmatex.cxx'\" unpacked with wrong size!
  2061. fi
  2062. # end of 'newmatex.cxx'
  2063. fi
  2064. echo shar: End of archive 4 \(of 8\).
  2065. cp /dev/null ark4isdone
  2066. MISSING=""
  2067. for I in 1 2 3 4 5 6 7 8 ; do
  2068.     if test ! -f ark${I}isdone ; then
  2069.     MISSING="${MISSING} ${I}"
  2070.     fi
  2071. done
  2072. if test "${MISSING}" = "" ; then
  2073.     echo You have unpacked all 8 archives.
  2074.     rm -f ark[1-9]isdone
  2075. else
  2076.     echo You still need to unpack the following archives:
  2077.     echo "        " ${MISSING}
  2078. fi
  2079. ##  End of shell archive.
  2080. exit 0
  2081.  
  2082. exit 0 # Just in case...
  2083.