home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1994 March / Source_Code_CD-ROM_Walnut_Creek_March_1994.iso / compsrcs / misc / volume34 / newmat07 / part03 < prev    next >
Encoding:
Text File  |  1993-01-10  |  61.4 KB  |  1,885 lines
  1. Newsgroups: comp.sources.misc
  2. From: robertd@kauri.vuw.ac.nz (Robert Davies)
  3. Subject: v34i109:  newmat07 - A matrix package in C++, Part03/08
  4. Message-ID: <1993Jan11.153101.2305@sparky.imd.sterling.com>
  5. X-Md4-Signature: 808bc6c5a49855ab8bd0e2b9325ffdb1
  6. Date: Mon, 11 Jan 1993 15:31:01 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: robertd@kauri.vuw.ac.nz (Robert Davies)
  10. Posting-number: Volume 34, Issue 109
  11. Archive-name: newmat07/part03
  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 3 (of 8)."
  22. # Contents:  ex_ms.mak ex_z.mak newmat.h newmatc.txt
  23. # Wrapped by robert@kea on Sun Jan 10 23:57:39 1993
  24. PATH=/bin:/usr/bin:/usr/ucb ; export PATH
  25. if test -f 'ex_ms.mak' -a "${1}" != "-c" ; then 
  26.   echo shar: Will not clobber existing file \"'ex_ms.mak'\"
  27. else
  28. echo shar: Extracting \"'ex_ms.mak'\" \(6300 characters\)
  29. sed "s/^X//" >'ex_ms.mak' <<'END_OF_FILE'
  30. XORIGIN = PWB
  31. XORIGIN_VER = 2.0
  32. XPROJ = EX_MS
  33. XPROJFILE = EX_MS.MAK
  34. XDEBUG = 0
  35. X
  36. XCC  = cl
  37. XCFLAGS_G  = /W2 /BATCH
  38. XCFLAGS_D  = /f /Zi /Od
  39. XCFLAGS_R  = /f- /Ot /Oi /Ol /Oe /Og /Gs
  40. XCXX  = cl
  41. XCXXFLAGS_G  = /AL /W2 /BATCH
  42. XCXXFLAGS_D  = /f /Od /Zi
  43. XCXXFLAGS_R  = /f- /Ot /Gs
  44. XMAPFILE_D  = NUL
  45. XMAPFILE_R  = NUL
  46. XLFLAGS_G  = /NOI /BATCH /ONERROR:NOEXE
  47. XLFLAGS_D  = /CO /FAR /PACKC
  48. XLFLAGS_R  = /EXE /FAR /PACKC
  49. XLINKER    = link
  50. XILINK  = ilink
  51. XLRF  = echo > NUL
  52. XILFLAGS  = /a /e
  53. X
  54. XFILES  = BANDMAT.CXX CHOLESKY.CXX EXAMPLE.CXX EXCEPT.CXX HHOLDER.CXX\
  55. X    NEWMAT1.CXX NEWMAT2.CXX NEWMAT3.CXX NEWMAT4.CXX NEWMAT5.CXX\
  56. X    NEWMAT6.CXX NEWMAT7.CXX NEWMAT8.CXX NEWMATEX.CXX SUBMAT.CXX SVD.CXX\
  57. X    NEWMATRM.CXX
  58. XOBJS  = BANDMAT.obj CHOLESKY.obj EXAMPLE.obj EXCEPT.obj HHOLDER.obj\
  59. X    NEWMAT1.obj NEWMAT2.obj NEWMAT3.obj NEWMAT4.obj NEWMAT5.obj\
  60. X    NEWMAT6.obj NEWMAT7.obj NEWMAT8.obj NEWMATEX.obj SUBMAT.obj SVD.obj\
  61. X    NEWMATRM.obj
  62. X
  63. Xall: $(PROJ).exe
  64. X
  65. X.SUFFIXES:
  66. X.SUFFIXES:
  67. X.SUFFIXES: .obj .cxx
  68. X
  69. XBANDMAT.obj : BANDMAT.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  70. X    controlw.h
  71. X!IF $(DEBUG)
  72. X    @$(CXX) @<<$(PROJ).rsp
  73. X/c $(CXXFLAGS_G)
  74. X$(CXXFLAGS_D) /FoBANDMAT.obj BANDMAT.CXX
  75. X<<
  76. X!ELSE
  77. X    @$(CXX) @<<$(PROJ).rsp
  78. X/c $(CXXFLAGS_G)
  79. X$(CXXFLAGS_R) /FoBANDMAT.obj BANDMAT.CXX
  80. X<<
  81. X!ENDIF
  82. X
  83. XCHOLESKY.obj : CHOLESKY.CXX include.h newmat.h boolean.h except.h
  84. X!IF $(DEBUG)
  85. X    @$(CXX) @<<$(PROJ).rsp
  86. X/c $(CXXFLAGS_G)
  87. X$(CXXFLAGS_D) /FoCHOLESKY.obj CHOLESKY.CXX
  88. X<<
  89. X!ELSE
  90. X    @$(CXX) @<<$(PROJ).rsp
  91. X/c $(CXXFLAGS_G)
  92. X$(CXXFLAGS_R) /FoCHOLESKY.obj CHOLESKY.CXX
  93. X<<
  94. X!ENDIF
  95. X
  96. XEXAMPLE.obj : EXAMPLE.CXX include.h newmatap.h newmat.h boolean.h except.h
  97. X!IF $(DEBUG)
  98. X    @$(CXX) @<<$(PROJ).rsp
  99. X/c $(CXXFLAGS_G)
  100. X$(CXXFLAGS_D) /FoEXAMPLE.obj EXAMPLE.CXX
  101. X<<
  102. X!ELSE
  103. X    @$(CXX) @<<$(PROJ).rsp
  104. X/c $(CXXFLAGS_G)
  105. X$(CXXFLAGS_R) /FoEXAMPLE.obj EXAMPLE.CXX
  106. X<<
  107. X!ENDIF
  108. X
  109. XEXCEPT.obj : EXCEPT.CXX include.h boolean.h except.h
  110. X!IF $(DEBUG)
  111. X    @$(CXX) @<<$(PROJ).rsp
  112. X/c $(CXXFLAGS_G)
  113. X$(CXXFLAGS_D) /FoEXCEPT.obj EXCEPT.CXX
  114. X<<
  115. X!ELSE
  116. X    @$(CXX) @<<$(PROJ).rsp
  117. X/c $(CXXFLAGS_G)
  118. X$(CXXFLAGS_R) /FoEXCEPT.obj EXCEPT.CXX
  119. X<<
  120. X!ENDIF
  121. X
  122. XHHOLDER.obj : HHOLDER.CXX include.h newmatap.h newmat.h boolean.h except.h
  123. X!IF $(DEBUG)
  124. X    @$(CXX) @<<$(PROJ).rsp
  125. X/c $(CXXFLAGS_G)
  126. X$(CXXFLAGS_D) /FoHHOLDER.obj HHOLDER.CXX
  127. X<<
  128. X!ELSE
  129. X    @$(CXX) @<<$(PROJ).rsp
  130. X/c $(CXXFLAGS_G)
  131. X$(CXXFLAGS_R) /FoHHOLDER.obj HHOLDER.CXX
  132. X<<
  133. X!ENDIF
  134. X
  135. XNEWMAT1.obj : NEWMAT1.CXX include.h newmat.h boolean.h except.h
  136. X!IF $(DEBUG)
  137. X    @$(CXX) @<<$(PROJ).rsp
  138. X/c $(CXXFLAGS_G)
  139. X$(CXXFLAGS_D) /FoNEWMAT1.obj NEWMAT1.CXX
  140. X<<
  141. X!ELSE
  142. X    @$(CXX) @<<$(PROJ).rsp
  143. X/c $(CXXFLAGS_G)
  144. X$(CXXFLAGS_R) /FoNEWMAT1.obj NEWMAT1.CXX
  145. X<<
  146. X!ENDIF
  147. X
  148. XNEWMAT2.obj : NEWMAT2.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  149. X    controlw.h
  150. X!IF $(DEBUG)
  151. X    @$(CXX) @<<$(PROJ).rsp
  152. X/c $(CXXFLAGS_G)
  153. X$(CXXFLAGS_D) /FoNEWMAT2.obj NEWMAT2.CXX
  154. X<<
  155. X!ELSE
  156. X    @$(CXX) @<<$(PROJ).rsp
  157. X/c $(CXXFLAGS_G)
  158. X$(CXXFLAGS_R) /FoNEWMAT2.obj NEWMAT2.CXX
  159. X<<
  160. X!ENDIF
  161. X
  162. XNEWMAT3.obj : NEWMAT3.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  163. X    controlw.h
  164. X!IF $(DEBUG)
  165. X    @$(CXX) @<<$(PROJ).rsp
  166. X/c $(CXXFLAGS_G)
  167. X$(CXXFLAGS_D) /FoNEWMAT3.obj NEWMAT3.CXX
  168. X<<
  169. X!ELSE
  170. X    @$(CXX) @<<$(PROJ).rsp
  171. X/c $(CXXFLAGS_G)
  172. X$(CXXFLAGS_R) /FoNEWMAT3.obj NEWMAT3.CXX
  173. X<<
  174. X!ENDIF
  175. X
  176. XNEWMAT4.obj : NEWMAT4.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  177. X    controlw.h
  178. X!IF $(DEBUG)
  179. X    @$(CXX) @<<$(PROJ).rsp
  180. X/c $(CXXFLAGS_G)
  181. X$(CXXFLAGS_D) /FoNEWMAT4.obj NEWMAT4.CXX
  182. X<<
  183. X!ELSE
  184. X    @$(CXX) @<<$(PROJ).rsp
  185. X/c $(CXXFLAGS_G)
  186. X$(CXXFLAGS_R) /FoNEWMAT4.obj NEWMAT4.CXX
  187. X<<
  188. X!ENDIF
  189. X
  190. XNEWMAT5.obj : NEWMAT5.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  191. X    controlw.h
  192. X!IF $(DEBUG)
  193. X    @$(CXX) @<<$(PROJ).rsp
  194. X/c $(CXXFLAGS_G)
  195. X$(CXXFLAGS_D) /FoNEWMAT5.obj NEWMAT5.CXX
  196. X<<
  197. X!ELSE
  198. X    @$(CXX) @<<$(PROJ).rsp
  199. X/c $(CXXFLAGS_G)
  200. X$(CXXFLAGS_R) /FoNEWMAT5.obj NEWMAT5.CXX
  201. X<<
  202. X!ENDIF
  203. X
  204. XNEWMAT6.obj : NEWMAT6.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  205. X    controlw.h
  206. X!IF $(DEBUG)
  207. X    @$(CXX) @<<$(PROJ).rsp
  208. X/c $(CXXFLAGS_G)
  209. X$(CXXFLAGS_D) /FoNEWMAT6.obj NEWMAT6.CXX
  210. X<<
  211. X!ELSE
  212. X    @$(CXX) @<<$(PROJ).rsp
  213. X/c $(CXXFLAGS_G)
  214. X$(CXXFLAGS_R) /FoNEWMAT6.obj NEWMAT6.CXX
  215. X<<
  216. X!ENDIF
  217. X
  218. XNEWMAT7.obj : NEWMAT7.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  219. X    controlw.h
  220. X!IF $(DEBUG)
  221. X    @$(CXX) @<<$(PROJ).rsp
  222. X/c $(CXXFLAGS_G)
  223. X$(CXXFLAGS_D) /FoNEWMAT7.obj NEWMAT7.CXX
  224. X<<
  225. X!ELSE
  226. X    @$(CXX) @<<$(PROJ).rsp
  227. X/c $(CXXFLAGS_G)
  228. X$(CXXFLAGS_R) /FoNEWMAT7.obj NEWMAT7.CXX
  229. X<<
  230. X!ENDIF
  231. X
  232. XNEWMAT8.obj : NEWMAT8.CXX include.h newmatap.h newmat.h boolean.h except.h
  233. X!IF $(DEBUG)
  234. X    @$(CXX) @<<$(PROJ).rsp
  235. X/c $(CXXFLAGS_G)
  236. X$(CXXFLAGS_D) /FoNEWMAT8.obj NEWMAT8.CXX
  237. X<<
  238. X!ELSE
  239. X    @$(CXX) @<<$(PROJ).rsp
  240. X/c $(CXXFLAGS_G)
  241. X$(CXXFLAGS_R) /FoNEWMAT8.obj NEWMAT8.CXX
  242. X<<
  243. X!ENDIF
  244. X
  245. XNEWMATEX.obj : NEWMATEX.CXX include.h newmat.h boolean.h except.h
  246. X!IF $(DEBUG)
  247. X    @$(CXX) @<<$(PROJ).rsp
  248. X/c $(CXXFLAGS_G)
  249. X$(CXXFLAGS_D) /FoNEWMATEX.obj NEWMATEX.CXX
  250. X<<
  251. X!ELSE
  252. X    @$(CXX) @<<$(PROJ).rsp
  253. X/c $(CXXFLAGS_G)
  254. X$(CXXFLAGS_R) /FoNEWMATEX.obj NEWMATEX.CXX
  255. X<<
  256. X!ENDIF
  257. X
  258. XSUBMAT.obj : SUBMAT.CXX include.h newmat.h newmatrc.h boolean.h except.h\
  259. X    controlw.h
  260. X!IF $(DEBUG)
  261. X    @$(CXX) @<<$(PROJ).rsp
  262. X/c $(CXXFLAGS_G)
  263. X$(CXXFLAGS_D) /FoSUBMAT.obj SUBMAT.CXX
  264. X<<
  265. X!ELSE
  266. X    @$(CXX) @<<$(PROJ).rsp
  267. X/c $(CXXFLAGS_G)
  268. X$(CXXFLAGS_R) /FoSUBMAT.obj SUBMAT.CXX
  269. X<<
  270. X!ENDIF
  271. X
  272. XSVD.obj : SVD.CXX include.h newmat.h newmatrm.h precisio.h boolean.h except.h
  273. X!IF $(DEBUG)
  274. X    @$(CXX) @<<$(PROJ).rsp
  275. X/c $(CXXFLAGS_G)
  276. X$(CXXFLAGS_D) /FoSVD.obj SVD.CXX
  277. X<<
  278. X!ELSE
  279. X    @$(CXX) @<<$(PROJ).rsp
  280. X/c $(CXXFLAGS_G)
  281. X$(CXXFLAGS_R) /FoSVD.obj SVD.CXX
  282. X<<
  283. X!ENDIF
  284. X
  285. XNEWMATRM.obj : NEWMATRM.CXX include.h newmat.h newmatrm.h boolean.h except.h
  286. X!IF $(DEBUG)
  287. X    @$(CXX) @<<$(PROJ).rsp
  288. X/c $(CXXFLAGS_G)
  289. X$(CXXFLAGS_D) /FoNEWMATRM.obj NEWMATRM.CXX
  290. X<<
  291. X!ELSE
  292. X    @$(CXX) @<<$(PROJ).rsp
  293. X/c $(CXXFLAGS_G)
  294. X$(CXXFLAGS_R) /FoNEWMATRM.obj NEWMATRM.CXX
  295. X<<
  296. X!ENDIF
  297. X
  298. X
  299. X$(PROJ).exe : $(OBJS)
  300. X!IF $(DEBUG)
  301. X    $(LRF) @<<$(PROJ).lrf
  302. X$(RT_OBJS: = +^
  303. X) $(OBJS: = +^
  304. X)
  305. X$@
  306. X$(MAPFILE_D)
  307. X$(LIBS: = +^
  308. X) +
  309. X$(LLIBS_G: = +^
  310. X) +
  311. X$(LLIBS_D: = +^
  312. X)
  313. X$(DEF_FILE) $(LFLAGS_G) $(LFLAGS_D);
  314. X<<
  315. X!ELSE
  316. X    $(LRF) @<<$(PROJ).lrf
  317. X$(RT_OBJS: = +^
  318. X) $(OBJS: = +^
  319. X)
  320. X$@
  321. X$(MAPFILE_R)
  322. X$(LIBS: = +^
  323. X) +
  324. X$(LLIBS_G: = +^
  325. X) +
  326. X$(LLIBS_R: = +^
  327. X)
  328. X$(DEF_FILE) $(LFLAGS_G) $(LFLAGS_R);
  329. X<<
  330. X!ENDIF
  331. X    $(LINKER) @$(PROJ).lrf
  332. X
  333. X
  334. X.cxx.obj :
  335. X!IF $(DEBUG)
  336. X    @$(CXX) @<<$(PROJ).rsp
  337. X/c $(CXXFLAGS_G)
  338. X$(CXXFLAGS_D) /Fo$@ $<
  339. X<<
  340. X!ELSE
  341. X    @$(CXX) @<<$(PROJ).rsp
  342. X/c $(CXXFLAGS_G)
  343. X$(CXXFLAGS_R) /Fo$@ $<
  344. X<<
  345. X!ENDIF
  346. X
  347. X
  348. Xrun: $(PROJ).exe
  349. X    $(PROJ).exe $(RUNFLAGS)
  350. X
  351. Xdebug: $(PROJ).exe
  352. X    CV $(CVFLAGS) $(PROJ).exe $(RUNFLAGS)
  353. END_OF_FILE
  354. if test 6300 -ne `wc -c <'ex_ms.mak'`; then
  355.     echo shar: \"'ex_ms.mak'\" unpacked with wrong size!
  356. fi
  357. # end of 'ex_ms.mak'
  358. fi
  359. if test -f 'ex_z.mak' -a "${1}" != "-c" ; then 
  360.   echo shar: Will not clobber existing file \"'ex_z.mak'\"
  361. else
  362. echo shar: Extracting \"'ex_z.mak'\" \(2558 characters\)
  363. sed "s/^X//" >'ex_z.mak' <<'END_OF_FILE'
  364. XC = ztc -c -ml $*.cxx
  365. X
  366. XOBJ = fft.obj evalue.obj submat.obj cholesky.obj hholder.obj          \
  367. X  sort.obj newmatrm.obj jacobi.obj svd.obj example.obj                \
  368. X  newmat8.obj newmat7.obj newmat6.obj                                 \
  369. X  newmat5.obj newmat3.obj newmat4.obj newmat2.obj newmat1.obj         \
  370. X  bandmat.obj except.obj newmatex.obj
  371. X
  372. Xex_z.exe:       $(OBJ) ex_z.lnk
  373. X                blink @ex_z.lnk
  374. X
  375. Xex_z.lnk:       ex_z.mak
  376. X            echo newmat1.obj+newmat2.obj+newmat3.obj+    > ex_z.lnk
  377. X            echo newmat4.obj+svd.obj+newmat5.obj+       >> ex_z.lnk
  378. X            echo newmat6.obj+newmat7.obj+newmat8.obj+   >> ex_z.lnk
  379. X            echo cholesky.obj+hholder.obj+sort.obj+     >> ex_z.lnk
  380. X            echo submat.obj+jacobi.obj+newmatrm.obj+    >> ex_z.lnk
  381. X            echo fft.obj+evalue.obj+bandmat.obj+        >> ex_z.lnk
  382. X            echo newmatex.obj+except.obj+example.obj    >> ex_z.lnk
  383. X            echo ex_z.exe                               >> ex_z.lnk
  384. X
  385. Xnewmatxx:       include.h newmat.h boolean.h except.h
  386. X            echo "main .h files uptodate?" > newmatxx
  387. X
  388. Xexcept.obj:     except.h except.cxx
  389. X            $C
  390. X
  391. Xnewmatex.obj:   newmatxx newmatex.cxx
  392. X            $C
  393. X
  394. Xexample.obj:    newmatxx newmatap.h example.cxx
  395. X            $C
  396. X
  397. Xcholesky.obj:   newmatxx cholesky.cxx
  398. X            $C
  399. X
  400. Xevalue.obj:     newmatxx newmatrm.h precisio.h evalue.cxx
  401. X            $C
  402. X
  403. Xfft.obj:        newmatxx newmatap.h fft.cxx
  404. X            $C
  405. X
  406. Xhholder.obj:    newmatxx newmatap.h hholder.cxx
  407. X            $C
  408. X
  409. Xjacobi.obj:     newmatxx precisio.h newmatrm.h jacobi.cxx
  410. X            $C
  411. X
  412. Xbandmat.obj:    newmatxx newmatrc.h controlw.h bandmat.cxx
  413. X            $C
  414. X
  415. Xnewmat1.obj:    newmatxx newmat1.cxx
  416. X            $C
  417. X
  418. Xnewmat2.obj:    newmatxx newmatrc.h controlw.h newmat2.cxx
  419. X            $C
  420. X
  421. Xnewmat3.obj:    newmatxx newmatrc.h controlw.h newmat3.cxx
  422. X            $C
  423. X
  424. Xnewmat4.obj:    newmatxx newmatrc.h controlw.h newmat4.cxx
  425. X            $C
  426. X
  427. Xnewmat5.obj:    newmatxx newmatrc.h controlw.h newmat5.cxx
  428. X            $C
  429. X
  430. Xnewmat6.obj:    newmatxx newmatrc.h controlw.h newmat6.cxx
  431. X            $C
  432. X
  433. Xnewmat7.obj:    newmatxx newmatrc.h controlw.h newmat7.cxx
  434. X            $C
  435. X
  436. Xnewmat8.obj:    newmatxx newmatap.h newmat8.cxx
  437. X            $C
  438. X
  439. Xnewmat9.obj:    newmatxx newmatrc.h controlw.h newmatio.h newmat9.cxx
  440. X            $C
  441. X
  442. Xnewmatrm.obj:   newmatxx newmatrm.h newmatrm.cxx
  443. X            $C
  444. X
  445. Xsort.obj:       newmatxx newmatap.h sort.cxx
  446. X            $C
  447. X
  448. Xsubmat.obj:     newmatxx newmatrc.h controlw.h submat.cxx
  449. X            $C
  450. X
  451. Xsvd.obj:        newmatxx newmatrm.h precisio.h svd.cxx
  452. X            $C
  453. X
  454. Xexample.obj:    newmatxx newmatap.h example.cxx 
  455. X            $C
  456. X
  457. END_OF_FILE
  458. if test 2558 -ne `wc -c <'ex_z.mak'`; then
  459.     echo shar: \"'ex_z.mak'\" unpacked with wrong size!
  460. fi
  461. # end of 'ex_z.mak'
  462. fi
  463. if test -f 'newmat.h' -a "${1}" != "-c" ; then 
  464.   echo shar: Will not clobber existing file \"'newmat.h'\"
  465. else
  466. echo shar: Extracting \"'newmat.h'\" \(47360 characters\)
  467. sed "s/^X//" >'newmat.h' <<'END_OF_FILE'
  468. X//$$ newmat.h           definition file for new version of matrix package
  469. X
  470. X// Copyright (C) 1991,2,3: R B Davies
  471. X
  472. X#ifndef MATRIX_LIB
  473. X#define MATRIX_LIB 0
  474. X
  475. X#ifdef NO_LONG_NAMES
  476. X#define UpperTriangularMatrix UTMatrix
  477. X#define LowerTriangularMatrix LTMatrix
  478. X#define SymmetricMatrix SMatrix
  479. X#define DiagonalMatrix DMatrix
  480. X#endif
  481. X
  482. X#ifndef TEMPS_DESTROYED_QUICKLY
  483. X#define ReturnMatrix ReturnMatrixX
  484. X#else
  485. X#define ReturnMatrix ReturnMatrixX&
  486. X#endif
  487. X
  488. X#include "boolean.h"
  489. X#include "except.h"
  490. X
  491. X/**************************** general utilities ****************************/
  492. X
  493. Xclass GeneralMatrix;
  494. Xvoid MatrixErrorNoSpace(void*);                 // no space handler
  495. X
  496. Xclass LogAndSign
  497. X// Return from LogDeterminant function
  498. X//    - value of the log plus the sign (+, - or 0)
  499. X{
  500. X   Real log_value;
  501. X   int sign;
  502. Xpublic:
  503. X   LogAndSign() { log_value=0.0; sign=1; }
  504. X   LogAndSign(Real);
  505. X   void operator*=(Real);
  506. X   void ChangeSign() { sign = -sign; }
  507. X   Real LogValue() const { return log_value; }
  508. X   int Sign() const { return sign; }
  509. X   Real Value() const;
  510. X   FREE_CHECK(LogAndSign)
  511. X};
  512. X
  513. X// the following class is for counting the number of times a piece of code
  514. X// is executed. It is used for locating any code not executed by test
  515. X// routines. Use turbo GREP locate all places this code is called and
  516. X// check which ones are not accessed.
  517. X// Somewhat implementation dependent as it relies on "cout" still being
  518. X// present when ExeCounter objects are destructed.
  519. X
  520. Xclass ExeCounter
  521. X{
  522. X   int line;                                    // code line number
  523. X   int fileid;                                  // file identifier
  524. X   long nexe;                                   // number of executions
  525. X   static int nreports;                         // number of reports
  526. Xpublic:
  527. X   ExeCounter(int,int);
  528. X   void operator++() { nexe++; }
  529. X   ~ExeCounter();                               // prints out reports
  530. X};
  531. X
  532. X
  533. X/************** Class to show whether to check for loss of data ************/
  534. X
  535. Xclass MatrixConversionCheck : public Janitor
  536. X{
  537. X   static Boolean DoCheck;
  538. Xpublic:
  539. X   MatrixConversionCheck() { DoCheck=TRUE; }          // turn check on
  540. X   ~MatrixConversionCheck() { DoCheck=FALSE; }        // turn check off
  541. X   void CleanUp() { DoCheck=FALSE; }
  542. X   static Boolean IsOn() { return DoCheck; }
  543. X   static void DataLoss();
  544. X};
  545. X
  546. X/**************************** class MatrixType *****************************/
  547. X
  548. X// Is used for finding the type of a matrix resulting from the binary operations
  549. X// +, -, * and identifying what conversions are permissible.
  550. X// This class must be updated when new matrix types are added.
  551. X
  552. Xclass GeneralMatrix;                            // defined later
  553. Xclass BaseMatrix;                               // defined later
  554. Xclass MatrixInput;                              // defined later
  555. X
  556. Xclass MatrixType
  557. X{
  558. Xpublic:
  559. X   enum Attribute {  Valid     = 1,
  560. X                     Symmetric = 2,
  561. X                     Band      = 4,
  562. X                     Upper     = 8,
  563. X                     Lower     = 16,
  564. X                     LUDeco    = 32 };
  565. X
  566. X   enum            { US = 0,
  567. X                     UT = Valid + Upper,
  568. X                     LT = Valid + Lower,
  569. X                     Rt = Valid,
  570. X                     Sm = Valid + Symmetric,
  571. X                     Dg = Valid + Band + Lower + Upper + Symmetric,
  572. X             RV = Valid,     //   don't separate out
  573. X             CV = Valid,     //   vectors
  574. X             BM = Valid + Band,
  575. X             UB = Valid + Band + Upper,
  576. X             LB = Valid + Band + Lower,
  577. X             SB = Valid + Band + Symmetric,
  578. X             Ct = Valid + LUDeco,
  579. X             BC = Valid + Band + LUDeco,
  580. X                   };
  581. X
  582. X
  583. X   static nTypes() { return 9; }               // number of different types
  584. X                           // exclude Ct, US, BC
  585. Xpublic:
  586. X   int attribute;
  587. Xpublic:
  588. X   MatrixType () {}
  589. X   MatrixType (int i) : attribute(i) {}
  590. X   int operator+() const { return attribute; }
  591. X   MatrixType operator+(MatrixType mt) const
  592. X      { return MatrixType(attribute & mt.attribute); }
  593. X   MatrixType operator*(const MatrixType&) const;
  594. X   Boolean operator>=(MatrixType mt) const
  595. X      { return ( attribute & mt.attribute ) == attribute; }
  596. X   Boolean operator==(MatrixType t) const
  597. X      { return (attribute == t.attribute); }
  598. X   Boolean operator!=(MatrixType t) const
  599. X      { return (attribute != t.attribute); }
  600. X   Boolean operator!() const { return (attribute & Valid) == 0; }
  601. X   MatrixType i() const                        // type of inverse
  602. X      { return MatrixType(attribute & ~(Band+LUDeco)); }
  603. X   MatrixType t() const;                       // type of transpose
  604. X   MatrixType AddEqualEl() const               // Add constant to matrix
  605. X      { return MatrixType(attribute & (Valid + Symmetric)); }
  606. X   MatrixType MultRHS() const;                 // type for rhs of multiply
  607. X   MatrixType sub() const                      // type of submatrix
  608. X      { return MatrixType(attribute & Valid); }
  609. X   MatrixType ssub() const                     // type of sym submatrix
  610. X      { return MatrixType(attribute); }        // not for selection matrix
  611. X   GeneralMatrix* New() const;                 // new matrix of given type
  612. X   GeneralMatrix* New(int,int,BaseMatrix*) const;
  613. X                                               // new matrix of given type
  614. X   operator char*() const;                     // to print type
  615. X   friend Boolean Rectangular(MatrixType a, MatrixType b, MatrixType c)
  616. X      { return ((a.attribute | b.attribute | c.attribute)
  617. X          & ~MatrixType::Valid) == 0; }
  618. X   friend Boolean Compare(const MatrixType&, MatrixType&);
  619. X                                               // compare and check conv.
  620. X   FREE_CHECK(MatrixType)
  621. X};
  622. X
  623. Xvoid TestTypeAdd();                            // test +
  624. Xvoid TestTypeMult();                           // test *
  625. Xvoid TestTypeOrder();                          // test >=
  626. X
  627. X
  628. X/************************* class MatrixBandWidth ***********************/
  629. X
  630. Xclass MatrixBandWidth
  631. X{
  632. Xpublic:
  633. X   int lower;
  634. X   int upper;
  635. X   MatrixBandWidth(const int l, const int u) : lower(l), upper (u) {}
  636. X   MatrixBandWidth(const int i) : lower(i), upper(i) {}
  637. X   MatrixBandWidth operator+(const MatrixBandWidth&) const;
  638. X   MatrixBandWidth operator*(const MatrixBandWidth&) const;
  639. X   MatrixBandWidth t() const { return MatrixBandWidth(upper,lower); }
  640. X   Boolean operator==(const MatrixBandWidth& bw) const
  641. X      { return (lower == bw.lower) && (upper == bw.upper); }
  642. X   FREE_CHECK(MatrixBandWidth)
  643. X};
  644. X
  645. X
  646. X/*********************** Array length specifier ************************/
  647. X
  648. X// This class is introduced to avoid constructors such as
  649. X//   ColumnVector(int)
  650. X// being used for conversions
  651. X
  652. Xclass ArrayLengthSpecifier
  653. X{
  654. X   int value;
  655. Xpublic:
  656. X   int Value() const { return value; }
  657. X   ArrayLengthSpecifier(int l) : value(l) {}
  658. X   FREE_CHECK(ArrayLengthSpecifier)
  659. X};
  660. X
  661. X/*************************** Matrix routines ***************************/
  662. X
  663. X
  664. Xclass MatrixRowCol;                             // defined later
  665. Xclass MatrixRow;
  666. Xclass MatrixCol;
  667. X
  668. Xclass GeneralMatrix;                            // defined later
  669. Xclass AddedMatrix;
  670. Xclass MultipliedMatrix;
  671. Xclass SubtractedMatrix;
  672. Xclass SolvedMatrix;
  673. Xclass ShiftedMatrix;
  674. Xclass ScaledMatrix;
  675. Xclass TransposedMatrix;
  676. Xclass NegatedMatrix;
  677. Xclass InvertedMatrix;
  678. Xclass RowedMatrix;
  679. Xclass ColedMatrix;
  680. Xclass DiagedMatrix;
  681. Xclass MatedMatrix;
  682. Xclass GetSubMatrix;
  683. Xclass ConstMatrix;
  684. Xclass ReturnMatrixX;
  685. Xclass Matrix;
  686. Xclass nricMatrix;
  687. Xclass RowVector;
  688. Xclass ColumnVector;
  689. Xclass SymmetricMatrix;
  690. Xclass UpperTriangularMatrix;
  691. Xclass LowerTriangularMatrix;
  692. Xclass DiagonalMatrix;
  693. Xclass CroutMatrix;
  694. Xclass BandMatrix;
  695. Xclass LowerBandMatrix;
  696. Xclass UpperBandMatrix;
  697. Xclass SymmetricBandMatrix;
  698. Xclass LinearEquationSolver;
  699. X
  700. X
  701. X
  702. Xstatic MatrixType MatrixTypeUnSp(MatrixType::US);
  703. X                        // AT&T needs this
  704. X
  705. Xclass BaseMatrix : public Janitor               // base of all matrix classes
  706. X{
  707. Xprotected:
  708. X   virtual int search(const BaseMatrix*) const = 0;
  709. X                        // count number of times matrix
  710. X                        // is referred to
  711. Xpublic:
  712. X#ifndef __GNUG__
  713. X   virtual GeneralMatrix* Evaluate(MatrixType mt=MatrixTypeUnSp) = 0;
  714. X                        // evaluate temporary
  715. X#else
  716. X   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
  717. X   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
  718. X#endif
  719. X#ifndef TEMPS_DESTROYED_QUICKLY
  720. X   AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
  721. X   MultipliedMatrix operator*(const BaseMatrix&) const;
  722. X   SubtractedMatrix operator-(const BaseMatrix&) const;
  723. X   ShiftedMatrix operator+(Real) const;
  724. X   ScaledMatrix operator*(Real) const;
  725. X   ScaledMatrix operator/(Real) const;
  726. X   ShiftedMatrix operator-(Real) const;
  727. X   TransposedMatrix t() const;
  728. X//   TransposedMatrix t;
  729. X   NegatedMatrix operator-() const;                   // change sign of elements
  730. X   InvertedMatrix i() const;
  731. X//   InvertedMatrix i;
  732. X   RowedMatrix AsRow() const;
  733. X   ColedMatrix AsColumn() const;
  734. X   DiagedMatrix AsDiagonal() const;
  735. X   MatedMatrix AsMatrix(int,int) const;
  736. X   GetSubMatrix SubMatrix(int,int,int,int) const;
  737. X   GetSubMatrix SymSubMatrix(int,int) const;
  738. X   GetSubMatrix Row(int) const;
  739. X   GetSubMatrix Rows(int,int) const;
  740. X   GetSubMatrix Column(int) const;
  741. X   GetSubMatrix Columns(int,int) const;
  742. X#else
  743. X   AddedMatrix& operator+(const BaseMatrix&) const;    // results of operations
  744. X   MultipliedMatrix& operator*(const BaseMatrix&) const;
  745. X   SubtractedMatrix& operator-(const BaseMatrix&) const;
  746. X   ShiftedMatrix& operator+(Real) const;
  747. X   ScaledMatrix& operator*(Real) const;
  748. X   ScaledMatrix& operator/(Real) const;
  749. X   ShiftedMatrix& operator-(Real) const;
  750. X   TransposedMatrix& t() const;
  751. X//   TransposedMatrix& t;
  752. X   NegatedMatrix& operator-() const;                   // change sign of elements
  753. X   InvertedMatrix& i() const;
  754. X//   InvertedMatrix& i;
  755. X   RowedMatrix& AsRow() const;
  756. X   ColedMatrix& AsColumn() const;
  757. X   DiagedMatrix& AsDiagonal() const;
  758. X   MatedMatrix& AsMatrix(int,int) const;
  759. X   GetSubMatrix& SubMatrix(int,int,int,int) const;
  760. X   GetSubMatrix& SymSubMatrix(int,int) const;
  761. X   GetSubMatrix& Row(int) const;
  762. X   GetSubMatrix& Rows(int,int) const;
  763. X   GetSubMatrix& Column(int) const;
  764. X   GetSubMatrix& Columns(int,int) const;
  765. X#endif
  766. X   Real AsScalar() const;                      // conversion of 1 x 1 matrix
  767. X   virtual LogAndSign LogDeterminant() const;
  768. X   virtual Real SumSquare() const;
  769. X   virtual Real SumAbsoluteValue() const;
  770. X   virtual Real MaximumAbsoluteValue() const;
  771. X   virtual Real Trace() const;
  772. X   Real Norm1() const;
  773. X   Real NormInfinity() const;
  774. X   virtual MatrixBandWidth BandWidth() const;  // bandwidths of band matrix
  775. X   virtual void CleanUp() {}                     // to clear store
  776. X//protected:
  777. X//   BaseMatrix() : t(this), i(this) {}
  778. X
  779. X   friend class GeneralMatrix;
  780. X   friend class Matrix;
  781. X   friend class nricMatrix;
  782. X   friend class RowVector;
  783. X   friend class ColumnVector;
  784. X   friend class SymmetricMatrix;
  785. X   friend class UpperTriangularMatrix;
  786. X   friend class LowerTriangularMatrix;
  787. X   friend class DiagonalMatrix;
  788. X   friend class CroutMatrix;
  789. X   friend class BandMatrix;
  790. X   friend class LowerBandMatrix;
  791. X   friend class UpperBandMatrix;
  792. X   friend class SymmetricBandMatrix;
  793. X   friend class AddedMatrix;
  794. X   friend class MultipliedMatrix;
  795. X   friend class SubtractedMatrix;
  796. X   friend class SolvedMatrix;
  797. X   friend class ShiftedMatrix;
  798. X   friend class ScaledMatrix;
  799. X   friend class TransposedMatrix;
  800. X   friend class NegatedMatrix;
  801. X   friend class InvertedMatrix;
  802. X   friend class RowedMatrix;
  803. X   friend class ColedMatrix;
  804. X   friend class DiagedMatrix;
  805. X   friend class MatedMatrix;
  806. X   friend class GetSubMatrix;
  807. X   friend class ConstMatrix;
  808. X   friend class ReturnMatrixX;
  809. X   friend class LinearEquationSolver;
  810. X   NEW_DELETE(BaseMatrix)
  811. X};
  812. X
  813. X
  814. X/******************************* working classes **************************/
  815. X
  816. Xclass GeneralMatrix : public BaseMatrix         // declarable matrix types
  817. X{
  818. X   virtual GeneralMatrix* Image() const;        // copy of matrix
  819. Xprotected:
  820. X   int tag;                                     // shows whether can reuse
  821. X   int nrows, ncols;                            // dimensions
  822. X   int storage;                                 // total store required
  823. X   Real* store;                                 // point to store (0=not set)
  824. X   GeneralMatrix();                             // initialise with no store
  825. X   GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
  826. X   void Add(GeneralMatrix*, Real);              // sum of GM and Real
  827. X   void Add(Real);                              // add Real to this
  828. X   void Multiply(GeneralMatrix*, Real);         // product of GM and Real
  829. X   void Multiply(Real);                         // multiply this by Real
  830. X   void Negate(GeneralMatrix*);                 // change sign
  831. X   void Negate();                               // change sign
  832. X   void operator=(Real);                        // set matrix to constant
  833. X   Real* GetStore();                            // get store or copy
  834. X   GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
  835. X                                                // temporarily access store
  836. X   void GetMatrix(const GeneralMatrix*);        // used by = and initialise
  837. X   void Eq(const BaseMatrix&, MatrixType);      // used by =
  838. X   int search(const BaseMatrix*) const;
  839. X   virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  840. X   void CheckConversion(const BaseMatrix&);     // check conversion OK
  841. X   void ReDimension(int, int, int);             // change dimensions
  842. Xpublic:
  843. X   GeneralMatrix* Evaluate(MatrixType);
  844. X   virtual MatrixType Type() const = 0;       // type of a matrix
  845. X   int Nrows() const { return nrows; }          // get dimensions
  846. X   int Ncols() const { return ncols; }
  847. X   int Storage() const { return storage; }
  848. X   Real* Store() const { return store; }
  849. X   virtual ~GeneralMatrix();                    // delete store if set
  850. X   void tDelete();                              // delete if tag permits
  851. X   Boolean reuse();                                // TRUE if tag allows reuse
  852. X   void Protect() { tag=-1; }                   // can't delete or reuse
  853. X   int Tag() const { return tag; }
  854. X   Boolean IsZero() const;                         // test matrix has all zeros
  855. X   void Release() { tag=1; }                    // del store after next use
  856. X   void Release(int t) { tag=t; }               // del store after t accesses
  857. X   void ReleaseAndDelete() { tag=0; }           // delete matrix after use
  858. X   void operator<<(const Real*);                // assignment from an array
  859. X   void operator<<(const BaseMatrix& X) { Eq(X,this->Type()); }
  860. X                                                // = without checking type
  861. X   void Inject(const GeneralMatrix&);           // copy stored els only
  862. X   virtual GeneralMatrix* MakeSolver();         // for solving
  863. X   virtual void Solver(MatrixRowCol&, const MatrixRowCol&) {}
  864. X   virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
  865. X   virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
  866. X   virtual void NextRow(MatrixRowCol&);         // Go to next row
  867. X   virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
  868. X   virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
  869. X   virtual void NextCol(MatrixRowCol&);         // Go to next col
  870. X   Real SumSquare() const;
  871. X   Real SumAbsoluteValue() const;
  872. X   Real MaximumAbsoluteValue() const;
  873. X   LogAndSign LogDeterminant() const;
  874. X#ifndef TEMPS_DESTROYED_QUICKLY
  875. X   ConstMatrix c() const;                       // to access constant matrices
  876. X#else
  877. X   ConstMatrix& c() const;                      // to access constant matrices
  878. X#endif
  879. X   void CheckStore() const;                     // check store is non-zero
  880. X   virtual void SetParameters(const GeneralMatrix*) {}
  881. X                                                // set parameters in GetMatrix
  882. X#ifndef TEMPS_DESTROYED_QUICKLY
  883. X   operator ReturnMatrixX() const;              // for building a ReturnMatrix
  884. X#else
  885. X   operator ReturnMatrixX&() const;             // for building a ReturnMatrix
  886. X#endif
  887. X   MatrixInput operator<<(Real);                // for loading a list
  888. X   void CleanUp();                              // to clear store
  889. X
  890. X   friend class Matrix;
  891. X   friend class nricMatrix;
  892. X   friend class SymmetricMatrix;
  893. X   friend class UpperTriangularMatrix;
  894. X   friend class LowerTriangularMatrix;
  895. X   friend class DiagonalMatrix;
  896. X   friend class CroutMatrix;
  897. X   friend class RowVector;
  898. X   friend class ColumnVector;
  899. X   friend class BandMatrix;
  900. X   friend class LowerBandMatrix;
  901. X   friend class UpperBandMatrix;
  902. X   friend class SymmetricBandMatrix;
  903. X   friend class BaseMatrix;
  904. X   friend class AddedMatrix;
  905. X   friend class MultipliedMatrix;
  906. X   friend class SubtractedMatrix;
  907. X   friend class SolvedMatrix;
  908. X   friend class ShiftedMatrix;
  909. X   friend class ScaledMatrix;
  910. X   friend class TransposedMatrix;
  911. X   friend class NegatedMatrix;
  912. X   friend class InvertedMatrix;
  913. X   friend class RowedMatrix;
  914. X   friend class ColedMatrix;
  915. X   friend class DiagedMatrix;
  916. X   friend class MatedMatrix;
  917. X   friend class GetSubMatrix;
  918. X   friend class ConstMatrix;
  919. X   friend class ReturnMatrixX;
  920. X   friend class LinearEquationSolver;
  921. X   NEW_DELETE(GeneralMatrix)
  922. X};
  923. X
  924. Xclass Matrix : public GeneralMatrix             // usual rectangular matrix
  925. X{
  926. X   GeneralMatrix* Image() const;                // copy of matrix
  927. Xpublic:
  928. X   Matrix() {}
  929. X   Matrix(int, int);                            // standard declaration
  930. X   Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
  931. X   void operator=(const BaseMatrix&);
  932. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  933. X   void operator=(const Matrix& m) { operator=((const BaseMatrix&)m); }
  934. X   MatrixType Type() const;
  935. X   Real& operator()(int, int);                  // access element
  936. X   Real& element(int, int);                     // access element
  937. X#ifdef SETUP_C_SUBSCRIPTS
  938. X   Real* operator[](int m) { return store+m*ncols; }
  939. X   const Real* operator[](int m) const { return store+m*ncols; }
  940. X#endif
  941. X   Matrix(const Matrix& gm) { GetMatrix(&gm); }
  942. X#ifndef __ZTC__
  943. X   Real operator()(int, int) const;             // access element
  944. X#endif
  945. X   GeneralMatrix* MakeSolver();
  946. X   Real Trace() const;
  947. X   void GetRow(MatrixRowCol&);
  948. X   void GetCol(MatrixRowCol&);
  949. X   void RestoreCol(MatrixRowCol&);
  950. X   void NextRow(MatrixRowCol&);
  951. X   void NextCol(MatrixRowCol&);
  952. X   void ReDimension(int,int);                   // change dimensions
  953. X   NEW_DELETE(Matrix)
  954. X};
  955. X
  956. Xclass nricMatrix : public Matrix                // for use with Numerical
  957. X                                                // Recipes in C
  958. X{
  959. X   GeneralMatrix* Image() const;                // copy of matrix
  960. X   Real** row_pointer;                          // points to rows
  961. X   void MakeRowPointer();                       // build rowpointer
  962. X   void DeleteRowPointer();
  963. Xpublic:
  964. X   nricMatrix() {}
  965. X   nricMatrix(int m, int n)                     // standard declaration
  966. X      :  Matrix(m,n) { MakeRowPointer(); }
  967. X   nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
  968. X      :  Matrix(bm) { MakeRowPointer(); }
  969. X   void operator=(const BaseMatrix& bm)
  970. X      { DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer(); }
  971. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  972. X   void operator=(const nricMatrix& m) { operator=((const BaseMatrix&)m); }
  973. X   void operator<<(const BaseMatrix& X)
  974. X      { DeleteRowPointer(); Eq(X,this->Type()); MakeRowPointer(); }
  975. X   nricMatrix(const nricMatrix& gm) { GetMatrix(&gm); MakeRowPointer(); }
  976. X   void ReDimension(int m, int n)               // change dimensions
  977. X      { DeleteRowPointer(); Matrix::ReDimension(m,n); MakeRowPointer(); }
  978. X   ~nricMatrix() { DeleteRowPointer(); }
  979. X#ifndef __ZTC__
  980. X   Real** nric() const { CheckStore(); return row_pointer-1; }
  981. X#endif
  982. X   void CleanUp();                                // to clear store
  983. X   NEW_DELETE(nricMatrix)
  984. X};
  985. X
  986. Xclass SymmetricMatrix : public GeneralMatrix
  987. X{
  988. X   GeneralMatrix* Image() const;                // copy of matrix
  989. Xpublic:
  990. X   SymmetricMatrix() {}
  991. X   SymmetricMatrix(ArrayLengthSpecifier);
  992. X   SymmetricMatrix(const BaseMatrix&);
  993. X   void operator=(const BaseMatrix&);
  994. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  995. X   void operator=(const SymmetricMatrix& m) { operator=((const BaseMatrix&)m); }
  996. X   Real& operator()(int, int);                  // access element
  997. X   Real& element(int, int);                     // access element
  998. X#ifdef SETUP_C_SUBSCRIPTS
  999. X   Real* operator[](int m) { return store+(m*(m+1))/2; }
  1000. X   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  1001. X#endif
  1002. X   MatrixType Type() const;
  1003. X   SymmetricMatrix(const SymmetricMatrix& gm) { GetMatrix(&gm); }
  1004. X#ifndef __ZTC__
  1005. X   Real operator()(int, int) const;             // access element
  1006. X#endif
  1007. X   Real SumSquare() const;
  1008. X   Real SumAbsoluteValue() const;
  1009. X   Real Trace() const;
  1010. X   void GetRow(MatrixRowCol&);
  1011. X   void GetCol(MatrixRowCol&);
  1012. X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  1013. X   void ReDimension(int);                       // change dimensions
  1014. X   NEW_DELETE(SymmetricMatrix)
  1015. X};
  1016. X
  1017. Xclass UpperTriangularMatrix : public GeneralMatrix
  1018. X{
  1019. X   GeneralMatrix* Image() const;                // copy of matrix
  1020. Xpublic:
  1021. X   UpperTriangularMatrix() {}
  1022. X   UpperTriangularMatrix(ArrayLengthSpecifier);
  1023. X   void operator=(const BaseMatrix&);
  1024. X   void operator=(const UpperTriangularMatrix& m)
  1025. X      { operator=((const BaseMatrix&)m); }
  1026. X   UpperTriangularMatrix(const BaseMatrix&);
  1027. X   UpperTriangularMatrix(const UpperTriangularMatrix& gm) { GetMatrix(&gm); }
  1028. X#ifndef __ZTC__
  1029. X   Real operator()(int, int) const;             // access element
  1030. X#endif
  1031. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1032. X   Real& operator()(int, int);                  // access element
  1033. X   Real& element(int, int);                     // access element
  1034. X#ifdef SETUP_C_SUBSCRIPTS
  1035. X   Real* operator[](int m) { return store+m*ncols-(m*(m+1))/2; }
  1036. X   const Real* operator[](int m) const { return store+m*ncols-(m*(m+1))/2; }
  1037. X#endif
  1038. X   MatrixType Type() const;
  1039. X   GeneralMatrix* MakeSolver() { return this; } // for solving
  1040. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1041. X   LogAndSign LogDeterminant() const;
  1042. X   Real Trace() const;
  1043. X   void GetRow(MatrixRowCol&);
  1044. X   void GetCol(MatrixRowCol&);
  1045. X   void RestoreCol(MatrixRowCol&);
  1046. X   void NextRow(MatrixRowCol&);
  1047. X   void ReDimension(int);                       // change dimensions
  1048. X   NEW_DELETE(UpperTriangularMatrix)
  1049. X};
  1050. X
  1051. Xclass LowerTriangularMatrix : public GeneralMatrix
  1052. X{
  1053. X   GeneralMatrix* Image() const;                // copy of matrix
  1054. Xpublic:
  1055. X   LowerTriangularMatrix() {}
  1056. X   LowerTriangularMatrix(ArrayLengthSpecifier);
  1057. X   LowerTriangularMatrix(const LowerTriangularMatrix& gm) { GetMatrix(&gm); }
  1058. X#ifndef __ZTC__
  1059. X   Real operator()(int, int) const;             // access element
  1060. X#endif
  1061. X   LowerTriangularMatrix(const BaseMatrix& M);
  1062. X   void operator=(const BaseMatrix&);
  1063. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1064. X   void operator=(const LowerTriangularMatrix& m)
  1065. X      { operator=((const BaseMatrix&)m); }
  1066. X   Real& operator()(int, int);                  // access element
  1067. X   Real& element(int, int);                     // access element
  1068. X#ifdef SETUP_C_SUBSCRIPTS
  1069. X   Real* operator[](int m) { return store+(m*(m+1))/2; }
  1070. X   const Real* operator[](int m) const { return store+(m*(m+1))/2; }
  1071. X#endif
  1072. X   MatrixType Type() const;
  1073. X   GeneralMatrix* MakeSolver() { return this; } // for solving
  1074. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1075. X   LogAndSign LogDeterminant() const;
  1076. X   Real Trace() const;
  1077. X   void GetRow(MatrixRowCol&);
  1078. X   void GetCol(MatrixRowCol&);
  1079. X   void RestoreCol(MatrixRowCol&);
  1080. X   void NextRow(MatrixRowCol&);
  1081. X   void ReDimension(int);                       // change dimensions
  1082. X   NEW_DELETE(LowerTriangularMatrix)
  1083. X};
  1084. X
  1085. Xclass DiagonalMatrix : public GeneralMatrix
  1086. X{
  1087. X   GeneralMatrix* Image() const;                // copy of matrix
  1088. Xpublic:
  1089. X   DiagonalMatrix() {}
  1090. X   DiagonalMatrix(ArrayLengthSpecifier);
  1091. X   DiagonalMatrix(const BaseMatrix&);
  1092. X   DiagonalMatrix(const DiagonalMatrix& gm) { GetMatrix(&gm); }
  1093. X#ifndef __ZTC__
  1094. X   Real operator()(int, int) const;             // access element
  1095. X   Real operator()(int) const;
  1096. X#endif
  1097. X   void operator=(const BaseMatrix&);
  1098. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1099. X   void operator=(const DiagonalMatrix& m) { operator=((const BaseMatrix&)m); }
  1100. X   Real& operator()(int, int);                  // access element
  1101. X   Real& operator()(int);                       // access element
  1102. X   Real& element(int, int);                     // access element
  1103. X   Real& element(int);                          // access element
  1104. X#ifdef SETUP_C_SUBSCRIPTS
  1105. X   Real& operator[](int m) { return store[m]; }
  1106. X   const Real& operator[](int m) const { return store[m]; }
  1107. X#endif
  1108. X   MatrixType Type() const;
  1109. X
  1110. X   LogAndSign LogDeterminant() const;
  1111. X   Real Trace() const;
  1112. X   void GetRow(MatrixRowCol&);
  1113. X   void GetCol(MatrixRowCol&);
  1114. X   void NextRow(MatrixRowCol&);
  1115. X   void NextCol(MatrixRowCol&);
  1116. X   GeneralMatrix* MakeSolver() { return this; } // for solving
  1117. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1118. X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  1119. X   void ReDimension(int);                       // change dimensions
  1120. X#ifndef __ZTC__
  1121. X   Real* nric() const
  1122. X      { CheckStore(); return store-1; }         // for use by NRIC
  1123. X#endif
  1124. X   MatrixBandWidth BandWidth() const;
  1125. X   NEW_DELETE(DiagonalMatrix)
  1126. X};
  1127. X
  1128. Xclass RowVector : public Matrix
  1129. X{
  1130. X   GeneralMatrix* Image() const;                // copy of matrix
  1131. Xpublic:
  1132. X   RowVector() {}
  1133. X   RowVector(ArrayLengthSpecifier n) : Matrix(1,n.Value()) {}
  1134. X   RowVector(const BaseMatrix&);
  1135. X   RowVector(const RowVector& gm) { GetMatrix(&gm); }
  1136. X#ifndef __ZTC__
  1137. X   Real operator()(int) const;                  // access element
  1138. X#endif
  1139. X   void operator=(const BaseMatrix&);
  1140. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1141. X   void operator=(const RowVector& m) { operator=((const BaseMatrix&)m); }
  1142. X   Real& operator()(int);                       // access element
  1143. X   Real& element(int);                          // access element
  1144. X#ifdef SETUP_C_SUBSCRIPTS
  1145. X   Real& operator[](int m) { return store[m]; }
  1146. X   const Real& operator[](int m) const { return store[m]; }
  1147. X#endif
  1148. X   MatrixType Type() const;
  1149. X   void GetCol(MatrixRowCol&);
  1150. X   void NextCol(MatrixRowCol&);
  1151. X   void RestoreCol(MatrixRowCol&);
  1152. X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  1153. X   void ReDimension(int);                       // change dimensions
  1154. X#ifndef __ZTC__
  1155. X   Real* nric() const
  1156. X      { CheckStore(); return store-1; }         // for use by NRIC
  1157. X#endif
  1158. X   void CleanUp();                                // to clear store
  1159. X   NEW_DELETE(RowVector)
  1160. X};
  1161. X
  1162. Xclass ColumnVector : public Matrix
  1163. X{
  1164. X   GeneralMatrix* Image() const;                // copy of matrix
  1165. Xpublic:
  1166. X   ColumnVector() {}
  1167. X   ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(),1) {}
  1168. X   ColumnVector(const BaseMatrix&);
  1169. X   ColumnVector(const ColumnVector& gm) { GetMatrix(&gm); }
  1170. X#ifndef __ZTC__
  1171. X   Real operator()(int) const;                  // access element
  1172. X#endif
  1173. X   void operator=(const BaseMatrix&);
  1174. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1175. X   void operator=(const ColumnVector& m) { operator=((const BaseMatrix&)m); }
  1176. X   Real& operator()(int);                       // access element
  1177. X   Real& element(int);                          // access element
  1178. X#ifdef SETUP_C_SUBSCRIPTS
  1179. X   Real& operator[](int m) { return store[m]; }
  1180. X   const Real& operator[](int m) const { return store[m]; }
  1181. X#endif
  1182. X   MatrixType Type() const;
  1183. X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  1184. X   void ReDimension(int);                       // change dimensions
  1185. X#ifndef __ZTC__
  1186. X   Real* nric() const
  1187. X      { CheckStore(); return store-1; }         // for use by NRIC
  1188. X#endif
  1189. X   void CleanUp();                                // to clear store
  1190. X   NEW_DELETE(ColumnVector)
  1191. X};
  1192. X
  1193. Xclass CroutMatrix : public GeneralMatrix        // for LU decomposition
  1194. X{
  1195. X   int* indx;
  1196. X   Boolean d;
  1197. X   Boolean sing;
  1198. X   void ludcmp();
  1199. Xpublic:
  1200. X   CroutMatrix(const BaseMatrix&);
  1201. X   MatrixType Type() const;
  1202. X   void lubksb(Real*, int=0);
  1203. X   ~CroutMatrix();
  1204. X   GeneralMatrix* MakeSolver() { return this; } // for solving
  1205. X   LogAndSign LogDeterminant() const;
  1206. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1207. X   void GetRow(MatrixRowCol&);
  1208. X   void GetCol(MatrixRowCol&);
  1209. X   void operator=(const BaseMatrix&);
  1210. X   void operator=(const CroutMatrix& m) { operator=((const BaseMatrix&)m); }
  1211. X   void CleanUp();                                // to clear store
  1212. X   NEW_DELETE(CroutMatrix)
  1213. X};
  1214. X
  1215. X/******************************* band matrices ***************************/
  1216. X
  1217. Xclass BandMatrix : public GeneralMatrix         // band matrix
  1218. X{
  1219. X   GeneralMatrix* Image() const;                // copy of matrix
  1220. Xprotected:
  1221. X   void CornerClear() const;                    // set unused elements to zero
  1222. Xpublic:
  1223. X   int lower, upper;                            // band widths
  1224. X   BandMatrix() { lower=0; upper=0; CornerClear(); }
  1225. X   BandMatrix(int n,int lb,int ub) { ReDimension(n,lb,ub); CornerClear(); }
  1226. X                                                // standard declaration
  1227. X   BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
  1228. X   void operator=(const BaseMatrix&);
  1229. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1230. X   void operator=(const BandMatrix& m) { operator=((const BaseMatrix&)m); }
  1231. X   MatrixType Type() const;
  1232. X   Real& operator()(int, int);                  // access element
  1233. X   Real& element(int, int);                     // access element
  1234. X#ifdef SETUP_C_SUBSCRIPTS
  1235. X   Real* operator[](int m) { return store+(upper+lower)*m+lower; }
  1236. X   const Real* operator[](int m) const { return store+(upper+lower)*m+lower; }
  1237. X#endif
  1238. X   BandMatrix(const BandMatrix& gm) { GetMatrix(&gm); }
  1239. X#ifndef __ZTC__
  1240. X   Real operator()(int, int) const;             // access element
  1241. X#endif
  1242. X   LogAndSign LogDeterminant() const;
  1243. X   GeneralMatrix* MakeSolver();
  1244. X   Real Trace() const;
  1245. X   Real SumSquare() const { CornerClear(); return GeneralMatrix::SumSquare(); }
  1246. X   Real SumAbsoluteValue() const
  1247. X      { CornerClear(); return GeneralMatrix::SumAbsoluteValue(); }
  1248. X   Real MaximumAbsoluteValue() const
  1249. X      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  1250. X   void GetRow(MatrixRowCol&);
  1251. X   void GetCol(MatrixRowCol&);
  1252. X   void RestoreCol(MatrixRowCol&);
  1253. X   void NextRow(MatrixRowCol&);
  1254. X   void ReDimension(int, int, int);             // change dimensions
  1255. X   MatrixBandWidth BandWidth() const;
  1256. X   void SetParameters(const GeneralMatrix*);
  1257. X   MatrixInput operator<<(Real);                // will give error
  1258. X   void operator<<(const Real* r);              // will give error
  1259. X      // the next is included because Zortech and Borland
  1260. X      // can't find the copy in GeneralMatrix
  1261. X   void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
  1262. X   NEW_DELETE(BandMatrix)
  1263. X};
  1264. X
  1265. Xclass UpperBandMatrix : public BandMatrix       // upper band matrix
  1266. X{
  1267. X   GeneralMatrix* Image() const;                // copy of matrix
  1268. Xpublic:
  1269. X   UpperBandMatrix() {}
  1270. X   UpperBandMatrix(int n, int ubw)              // standard declaration
  1271. X      : BandMatrix(n, 0, ubw) {}
  1272. X   UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  1273. X   void operator=(const BaseMatrix&);
  1274. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1275. X   void operator=(const UpperBandMatrix& m)
  1276. X      { operator=((const BaseMatrix&)m); }
  1277. X   MatrixType Type() const;
  1278. X   UpperBandMatrix(const UpperBandMatrix& gm) { GetMatrix(&gm); }
  1279. X   GeneralMatrix* MakeSolver() { return this; }
  1280. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1281. X   LogAndSign LogDeterminant() const;
  1282. X   void ReDimension(int n,int ubw)              // change dimensions
  1283. X      { BandMatrix::ReDimension(n,0,ubw); }
  1284. X   Real& operator()(int, int);
  1285. X   Real& element(int, int);
  1286. X#ifdef SETUP_C_SUBSCRIPTS
  1287. X   Real* operator[](int m) { return store+upper*m; }
  1288. X   const Real* operator[](int m) const { return store+upper*m; }
  1289. X#endif
  1290. X   NEW_DELETE(UpperBandMatrix)
  1291. X};
  1292. X
  1293. Xclass LowerBandMatrix : public BandMatrix       // upper band matrix
  1294. X{
  1295. X   GeneralMatrix* Image() const;                // copy of matrix
  1296. Xpublic:
  1297. X   LowerBandMatrix() {}
  1298. X   LowerBandMatrix(int n, int lbw)              // standard declaration
  1299. X      : BandMatrix(n, lbw, 0) {}
  1300. X   LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
  1301. X   void operator=(const BaseMatrix&);
  1302. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1303. X   void operator=(const LowerBandMatrix& m)
  1304. X      { operator=((const BaseMatrix&)m); }
  1305. X   MatrixType Type() const;
  1306. X   LowerBandMatrix(const LowerBandMatrix& gm) { GetMatrix(&gm); }
  1307. X   GeneralMatrix* MakeSolver() { return this; }
  1308. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1309. X   LogAndSign LogDeterminant() const;
  1310. X   void ReDimension(int n,int lbw)             // change dimensions
  1311. X      { BandMatrix::ReDimension(n,lbw,0); }
  1312. X   Real& operator()(int, int);
  1313. X   Real& element(int, int);
  1314. X#ifdef SETUP_C_SUBSCRIPTS
  1315. X   Real* operator[](int m) { return store+lower*(m+1); }
  1316. X   const Real* operator[](int m) const { return store+lower*(m+1); }
  1317. X#endif
  1318. X   NEW_DELETE(LowerBandMatrix)
  1319. X};
  1320. X
  1321. Xclass SymmetricBandMatrix : public GeneralMatrix
  1322. X{
  1323. X   GeneralMatrix* Image() const;                // copy of matrix
  1324. X   void CornerClear() const;                    // set unused elements to zero
  1325. Xpublic:
  1326. X   int lower;                                   // lower band width
  1327. X   SymmetricBandMatrix() { lower=0; CornerClear(); }
  1328. X   SymmetricBandMatrix(int n, int lb) { ReDimension(n,lb); CornerClear(); }
  1329. X   SymmetricBandMatrix(const BaseMatrix&);
  1330. X   void operator=(const BaseMatrix&);
  1331. X   void operator=(Real f) { GeneralMatrix::operator=(f); }
  1332. X   void operator=(const SymmetricBandMatrix& m)
  1333. X      { operator=((const BaseMatrix&)m); }
  1334. X   Real& operator()(int, int);                  // access element
  1335. X   Real& element(int, int);                     // access element
  1336. X#ifdef SETUP_C_SUBSCRIPTS
  1337. X   Real* operator[](int m) { return store+lower*(m+1); }
  1338. X   const Real* operator[](int m) const { return store+lower*(m+1); }
  1339. X#endif
  1340. X   MatrixType Type() const;
  1341. X   SymmetricBandMatrix(const SymmetricBandMatrix& gm) { GetMatrix(&gm); }
  1342. X#ifndef __ZTC__
  1343. X   Real operator()(int, int) const;             // access element
  1344. X#endif
  1345. X   GeneralMatrix* MakeSolver();
  1346. X   Real SumSquare() const;
  1347. X   Real SumAbsoluteValue() const;
  1348. X   Real MaximumAbsoluteValue() const
  1349. X      { CornerClear(); return GeneralMatrix::MaximumAbsoluteValue(); }
  1350. X   Real Trace() const;
  1351. X   LogAndSign LogDeterminant() const;
  1352. X   void GetRow(MatrixRowCol&);
  1353. X   void GetCol(MatrixRowCol&);
  1354. X   GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
  1355. X   void ReDimension(int,int);                       // change dimensions
  1356. X   MatrixBandWidth BandWidth() const;
  1357. X   void SetParameters(const GeneralMatrix*);
  1358. X   NEW_DELETE(SymmetricBandMatrix)
  1359. X};
  1360. X
  1361. Xclass BandLUMatrix : public GeneralMatrix
  1362. X// for LU decomposition of band matrix
  1363. X{
  1364. X   int* indx;
  1365. X   Boolean d;
  1366. X   Boolean sing;                                   // TRUE if singular
  1367. X   Real* store2;
  1368. X   int storage2;
  1369. X   void ludcmp();
  1370. X   int m1,m2;                                   // lower and upper
  1371. Xpublic:
  1372. X   BandLUMatrix(const BaseMatrix&);
  1373. X   MatrixType Type() const;
  1374. X   void lubksb(Real*, int=0);
  1375. X   ~BandLUMatrix();
  1376. X   GeneralMatrix* MakeSolver() { return this; } // for solving
  1377. X   LogAndSign LogDeterminant() const;
  1378. X   void Solver(MatrixRowCol&, const MatrixRowCol&);
  1379. X   void GetRow(MatrixRowCol&);
  1380. X   void GetCol(MatrixRowCol&);
  1381. X   void operator=(const BaseMatrix&);
  1382. X   void operator=(const BandLUMatrix& m) { operator=((const BaseMatrix&)m); }
  1383. X   NEW_DELETE(BandLUMatrix)
  1384. X   void CleanUp();                                // to clear store
  1385. X};
  1386. X
  1387. X
  1388. X/***************************** temporary classes *************************/
  1389. X
  1390. Xclass MultipliedMatrix : public BaseMatrix
  1391. X{
  1392. Xprotected:
  1393. X   union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
  1394. X                          // pointers to summands
  1395. X   union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
  1396. X   MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1397. X      : bm1(bm1x),bm2(bm2x) {}
  1398. X   int search(const BaseMatrix*) const;
  1399. X   friend class BaseMatrix;
  1400. Xpublic:
  1401. X   GeneralMatrix* Evaluate(MatrixType);
  1402. X   MatrixBandWidth BandWidth() const;
  1403. X   NEW_DELETE(MultipliedMatrix)
  1404. X};
  1405. X
  1406. Xclass AddedMatrix : public MultipliedMatrix
  1407. X{
  1408. Xprotected:
  1409. X   AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1410. X      : MultipliedMatrix(bm1x,bm2x) {}
  1411. X
  1412. X   friend class BaseMatrix;
  1413. Xpublic:
  1414. X   GeneralMatrix* Evaluate(MatrixType);
  1415. X   MatrixBandWidth BandWidth() const;
  1416. X#ifdef __GNUG__
  1417. X   void SelectVersion(MatrixType, int&, int&) const;
  1418. X#else
  1419. X   void SelectVersion(MatrixType, Boolean&, Boolean&) const;
  1420. X#endif
  1421. X   NEW_DELETE(AddedMatrix)
  1422. X};
  1423. X
  1424. Xclass SolvedMatrix : public MultipliedMatrix
  1425. X{
  1426. X   SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1427. X      : MultipliedMatrix(bm1x,bm2x) {}
  1428. X   friend class BaseMatrix;
  1429. X   friend class InvertedMatrix;                        // for operator*
  1430. Xpublic:
  1431. X   GeneralMatrix* Evaluate(MatrixType);
  1432. X   MatrixBandWidth BandWidth() const;
  1433. X   NEW_DELETE(SolvedMatrix)
  1434. X};
  1435. X
  1436. Xclass SubtractedMatrix : public AddedMatrix
  1437. X{
  1438. X   SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
  1439. X      : AddedMatrix(bm1x,bm2x) {}
  1440. X   friend class BaseMatrix;
  1441. Xpublic:
  1442. X   GeneralMatrix* Evaluate(MatrixType);
  1443. X   NEW_DELETE(SubtractedMatrix)
  1444. X};
  1445. X
  1446. Xclass ShiftedMatrix : public BaseMatrix
  1447. X{
  1448. Xprotected:
  1449. X   Real f;
  1450. X   union { const BaseMatrix* bm; GeneralMatrix* gm; };
  1451. X   ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx),f(fx) {}
  1452. X   int search(const BaseMatrix*) const;
  1453. X   friend class BaseMatrix;
  1454. Xpublic:
  1455. X   GeneralMatrix* Evaluate(MatrixType);
  1456. X   NEW_DELETE(ShiftedMatrix)
  1457. X};
  1458. X
  1459. Xclass ScaledMatrix : public ShiftedMatrix
  1460. X{
  1461. X   ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx,fx) {}
  1462. X   friend class BaseMatrix;
  1463. Xpublic:
  1464. X   GeneralMatrix* Evaluate(MatrixType);
  1465. X   MatrixBandWidth BandWidth() const;
  1466. X   NEW_DELETE(ScaledMatrix)
  1467. X};
  1468. X
  1469. Xclass NegatedMatrix : public BaseMatrix
  1470. X{
  1471. Xprotected:
  1472. X   union { const BaseMatrix* bm; GeneralMatrix* gm; };
  1473. X   NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
  1474. X   int search(const BaseMatrix*) const;
  1475. Xprivate:
  1476. X   friend class BaseMatrix;
  1477. Xpublic:
  1478. X   GeneralMatrix* Evaluate(MatrixType);
  1479. X   MatrixBandWidth BandWidth() const;
  1480. X   NEW_DELETE(NegatedMatrix)
  1481. X};
  1482. X
  1483. Xclass TransposedMatrix : public NegatedMatrix
  1484. X{
  1485. X   TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1486. X   friend class BaseMatrix;
  1487. Xpublic:
  1488. X   GeneralMatrix* Evaluate(MatrixType);
  1489. X   MatrixBandWidth BandWidth() const;
  1490. X   NEW_DELETE(TransposedMatrix)
  1491. X};
  1492. X
  1493. Xclass InvertedMatrix : public NegatedMatrix
  1494. X{
  1495. X   InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1496. Xpublic:
  1497. X#ifndef TEMPS_DESTROYED_QUICKLY
  1498. X   SolvedMatrix operator*(const BaseMatrix&) const;  // inverse(A) * B
  1499. X#else
  1500. X   SolvedMatrix& operator*(const BaseMatrix&);       // inverse(A) * B
  1501. X#endif
  1502. X   friend class BaseMatrix;
  1503. X   GeneralMatrix* Evaluate(MatrixType);
  1504. X   MatrixBandWidth BandWidth() const;
  1505. X   NEW_DELETE(InvertedMatrix)
  1506. X};
  1507. X
  1508. Xclass RowedMatrix : public NegatedMatrix
  1509. X{
  1510. X   RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1511. X   friend class BaseMatrix;
  1512. Xpublic:
  1513. X   GeneralMatrix* Evaluate(MatrixType);
  1514. X   MatrixBandWidth BandWidth() const;
  1515. X   NEW_DELETE(RowedMatrix)
  1516. X};
  1517. X
  1518. Xclass ColedMatrix : public NegatedMatrix
  1519. X{
  1520. X   ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1521. X   friend class BaseMatrix;
  1522. Xpublic:
  1523. X   GeneralMatrix* Evaluate(MatrixType);
  1524. X   MatrixBandWidth BandWidth() const;
  1525. X   NEW_DELETE(ColedMatrix)
  1526. X};
  1527. X
  1528. Xclass DiagedMatrix : public NegatedMatrix
  1529. X{
  1530. X   DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
  1531. X   friend class BaseMatrix;
  1532. Xpublic:
  1533. X   GeneralMatrix* Evaluate(MatrixType);
  1534. X   MatrixBandWidth BandWidth() const;
  1535. X   NEW_DELETE(DiagedMatrix)
  1536. X};
  1537. X
  1538. Xclass MatedMatrix : public NegatedMatrix
  1539. X{
  1540. X   int nr, nc;
  1541. X   MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
  1542. X      : NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
  1543. X   friend class BaseMatrix;
  1544. Xpublic:
  1545. X   GeneralMatrix* Evaluate(MatrixType);
  1546. X   MatrixBandWidth BandWidth() const;
  1547. X   NEW_DELETE(MatedMatrix)
  1548. X};
  1549. X
  1550. Xclass ConstMatrix : public BaseMatrix
  1551. X{
  1552. X   const GeneralMatrix* cgm;
  1553. X   int search(const BaseMatrix*) const;
  1554. X   ConstMatrix(const GeneralMatrix* cgmx) : cgm(cgmx) {}
  1555. X   friend class BaseMatrix;
  1556. X   friend class GeneralMatrix;
  1557. Xpublic:
  1558. X   GeneralMatrix* Evaluate(MatrixType);
  1559. X   MatrixBandWidth BandWidth() const;
  1560. X   NEW_DELETE(ConstMatrix)
  1561. X};
  1562. X
  1563. Xclass ReturnMatrixX : public BaseMatrix    // for matrix return
  1564. X{
  1565. X   GeneralMatrix* gm;
  1566. X   int search(const BaseMatrix*) const;
  1567. Xpublic:
  1568. X   GeneralMatrix* Evaluate(MatrixType);
  1569. X   friend class BaseMatrix;
  1570. X#ifdef TEMPS_DESTROYED_QUICKLY
  1571. X   ReturnMatrixX(const ReturnMatrixX& tm);
  1572. X#else
  1573. X   ReturnMatrixX(const ReturnMatrixX& tm) : gm(tm.gm) {}
  1574. X#endif
  1575. X   ReturnMatrixX(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
  1576. X//   ReturnMatrixX(GeneralMatrix&);
  1577. X   MatrixBandWidth BandWidth() const;
  1578. X   NEW_DELETE(ReturnMatrixX)
  1579. X};
  1580. X
  1581. X
  1582. X/**************************** submatrices ******************************/
  1583. X
  1584. Xclass GetSubMatrix : public NegatedMatrix
  1585. X{
  1586. X   int row_skip;
  1587. X   int row_number;
  1588. X   int col_skip;
  1589. X   int col_number;
  1590. X   Boolean IsSym;
  1591. X
  1592. X   GetSubMatrix
  1593. X      (const BaseMatrix* bmx, int rs, int rn, int cs, int cn, Boolean is)
  1594. X      : NegatedMatrix(bmx),
  1595. X      row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
  1596. X   GetSubMatrix(const GetSubMatrix& g)
  1597. X      : NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
  1598. X      col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
  1599. X   void SetUpLHS();
  1600. X   friend class BaseMatrix;
  1601. Xpublic:
  1602. X   GeneralMatrix* Evaluate(MatrixType);
  1603. X   void operator=(const BaseMatrix&);
  1604. X   void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
  1605. X   void operator<<(const BaseMatrix&);
  1606. X   void operator<<(const Real*);                // copy from array
  1607. X   void operator=(Real);                        // copy from constant
  1608. X   void Inject(const GeneralMatrix&);           // copy stored els only
  1609. X   MatrixBandWidth BandWidth() const;
  1610. X   NEW_DELETE(GetSubMatrix)
  1611. X};
  1612. X
  1613. X/**************************** matrix input *******************************/
  1614. X
  1615. Xclass MatrixInput          // for reading a list of values into a matrix
  1616. X                           // the difficult part is detecting a mismatch
  1617. X                           // in the number of elements
  1618. X{
  1619. X   static int n;           // number values still to be read
  1620. X   static Real* r;         // pointer to last thing read
  1621. X   static int depth;       // number of objects of this class in existence
  1622. Xpublic:
  1623. X   MatrixInput() { depth++; }
  1624. X   MatrixInput(const MatrixInput&) { depth++; }
  1625. X   ~MatrixInput();
  1626. X   MatrixInput operator<<(Real);
  1627. X                           // could return a reference if we don't have
  1628. X                           // TEMPS_DESTROYED_QUICKLY
  1629. X   friend GeneralMatrix;                             
  1630. X};
  1631. X
  1632. X/***************************** exceptions ********************************/
  1633. X
  1634. Xclass MatrixDetails
  1635. X{
  1636. X   MatrixType type;
  1637. X   int nrows;
  1638. X   int ncols;
  1639. X   int ubw;
  1640. X   int lbw;
  1641. Xpublic:
  1642. X   MatrixDetails(const GeneralMatrix& A);
  1643. X   void PrintOut();
  1644. X};
  1645. X   
  1646. X
  1647. X
  1648. Xclass SpaceException : public Exception
  1649. X{
  1650. Xpublic:
  1651. X   static long st_type() { return 2; }
  1652. X   long type() const { return 2; }
  1653. X   static int action;
  1654. X   SpaceException();
  1655. X   static void SetAction(int a) { action=a; }
  1656. X};
  1657. X
  1658. X
  1659. Xclass MatrixException : public Exception
  1660. X{
  1661. Xpublic:
  1662. X   static long st_type() { return 3; }
  1663. X   long type() const { return 3; }
  1664. X   MatrixException(int);
  1665. X   MatrixException(int, const GeneralMatrix&);
  1666. X   MatrixException(int, const GeneralMatrix&, const GeneralMatrix&);
  1667. X};
  1668. X
  1669. Xclass DataException : public MatrixException
  1670. X{
  1671. Xpublic:
  1672. X   static long st_type() { return 3*53; }
  1673. X   long type() const { return 3*53; }
  1674. X   static int action;
  1675. X   DataException(const GeneralMatrix& A);
  1676. X   static void SetAction(int a) { action=a; }
  1677. X};
  1678. X
  1679. Xclass SingularException : public DataException
  1680. X{
  1681. Xpublic:
  1682. X   static long st_type() { return 3*53*109; }
  1683. X   long type() const { return 3*53*109; }
  1684. X   SingularException(const GeneralMatrix& A);
  1685. X};
  1686. X
  1687. Xclass NPDException : public DataException     // Not positive definite
  1688. X{
  1689. Xpublic:
  1690. X   static long st_type() { return 3*53*113; }
  1691. X   long type() const { return 3*53*113; }
  1692. X   NPDException(const GeneralMatrix&);
  1693. X};
  1694. X
  1695. Xclass ConvergenceException : public MatrixException
  1696. X{
  1697. Xpublic:
  1698. X   static long st_type() { return 3*59; }
  1699. X   long type() const { return 3*59; }
  1700. X   static int action;
  1701. X   ConvergenceException(const GeneralMatrix& A);
  1702. X   static void SetAction(int a) { action=a; }
  1703. X};
  1704. X
  1705. Xclass ProgramException : public MatrixException
  1706. X{
  1707. Xpublic:
  1708. X   static long st_type() { return 3*61; }
  1709. X   long type() const { return 3*61; }
  1710. X   static int action;
  1711. X   ProgramException(char* c);
  1712. X   ProgramException(char* c, const GeneralMatrix&);
  1713. X   ProgramException(char* c, const GeneralMatrix&, const GeneralMatrix&);
  1714. X   ProgramException();
  1715. X   ProgramException(const GeneralMatrix&);
  1716. X   static void SetAction(int a) { action=a; }
  1717. X};
  1718. X
  1719. Xclass IndexException : public ProgramException
  1720. X{
  1721. Xpublic:
  1722. X   static long st_type() { return 3*61*101; }
  1723. X   long type() const { return 3*61*101; }
  1724. X   IndexException(int i, const GeneralMatrix& A);
  1725. X   IndexException(int i, int j, const GeneralMatrix& A);
  1726. X   // next two are for access via element function
  1727. X   IndexException(int i, const GeneralMatrix& A, Boolean);
  1728. X   IndexException(int i, int j, const GeneralMatrix& A, Boolean);
  1729. X};
  1730. X
  1731. Xclass VectorException : public ProgramException    // can't convert to vector
  1732. X{
  1733. Xpublic:
  1734. X   static long st_type() { return 3*61*107; }
  1735. X   long type() const { return 3*61*107; }
  1736. X   VectorException();
  1737. X   VectorException(const GeneralMatrix& A);
  1738. X};
  1739. X
  1740. Xclass NotSquareException : public ProgramException
  1741. X{
  1742. Xpublic:
  1743. X   static long st_type() { return 3*61*109; }
  1744. X   long type() const { return 3*61*109; }
  1745. X   NotSquareException(const GeneralMatrix& A);
  1746. X};
  1747. X
  1748. Xclass SubMatrixDimensionException : public ProgramException
  1749. X{
  1750. Xpublic:
  1751. X   static long st_type() { return 3*61*113; }
  1752. X   long type() const { return 3*61*113; }
  1753. X   SubMatrixDimensionException();
  1754. X};
  1755. X
  1756. Xclass IncompatibleDimensionsException : public ProgramException
  1757. X{
  1758. Xpublic:
  1759. X   static long st_type() { return 3*61*127; }
  1760. X   long type() const { return 3*61*127; }
  1761. X   IncompatibleDimensionsException();
  1762. X};
  1763. X
  1764. Xclass NotDefinedException : public ProgramException
  1765. X{
  1766. Xpublic:
  1767. X   static long st_type() { return 3*61*131; }
  1768. X   long type() const { return 3*61*131; }
  1769. X   NotDefinedException(char* op, char* matrix);
  1770. X};
  1771. X
  1772. Xclass CannotBuildException : public ProgramException
  1773. X{
  1774. Xpublic:
  1775. X   static long st_type() { return 3*61*137; }
  1776. X   long type() const { return 3*61*137; }
  1777. X   CannotBuildException(char* matrix);
  1778. X};
  1779. X
  1780. X
  1781. Xclass InternalException : public MatrixException
  1782. X{
  1783. Xpublic:
  1784. X   static long st_type() { return 3*67; }
  1785. X   long type() const { return 3*67; }
  1786. X   static int action;
  1787. X   InternalException(char* c);
  1788. X   static void SetAction(int a) { action=a; }
  1789. X};
  1790. X
  1791. X
  1792. X/***************************** functions ***********************************/
  1793. X
  1794. X
  1795. Xinline LogAndSign LogDeterminant(const BaseMatrix& B)
  1796. X   { return B.LogDeterminant(); }
  1797. Xinline Real SumSquare(const BaseMatrix& B) { return B.SumSquare(); }
  1798. Xinline Real Trace(const BaseMatrix& B) { return B.Trace(); }
  1799. Xinline Real SumAbsoluteValue(const BaseMatrix& B)
  1800. X   { return B.SumAbsoluteValue(); }
  1801. Xinline Real MaximumAbsoluteValue(const BaseMatrix& B)
  1802. X   { return B.MaximumAbsoluteValue(); }
  1803. Xinline Real Norm1(const BaseMatrix& B) { return B.Norm1(); }
  1804. Xinline Real Norm1(RowVector& RV) { return RV.MaximumAbsoluteValue(); }
  1805. Xinline Real NormInfinity(const BaseMatrix& B) { return B.NormInfinity(); }
  1806. Xinline Real NormInfinity(ColumnVector& CV)
  1807. X   { return CV.MaximumAbsoluteValue(); } 
  1808. X
  1809. X#endif
  1810. END_OF_FILE
  1811. if test 47360 -ne `wc -c <'newmat.h'`; then
  1812.     echo shar: \"'newmat.h'\" unpacked with wrong size!
  1813. fi
  1814. # end of 'newmat.h'
  1815. fi
  1816. if test -f 'newmatc.txt' -a "${1}" != "-c" ; then 
  1817.   echo shar: Will not clobber existing file \"'newmatc.txt'\"
  1818. else
  1819. echo shar: Extracting \"'newmatc.txt'\" \(1948 characters\)
  1820. sed "s/^X//" >'newmatc.txt' <<'END_OF_FILE'
  1821. X//$$ newmatc.txt                                      Testing
  1822. X
  1823. X                    Testing newmat on your compiler
  1824. X                    ===============================
  1825. X
  1826. XThere are a series of files of the form tmt?.cxx included with some
  1827. Xversions of the package. These are used to check that the package is
  1828. Xperforming correctly. They consist of a large number of matrix formulae
  1829. Xall of which evaluate to zero (except the first one which is used to
  1830. Xcheck that we are detecting non-zero matrices). The printout should
  1831. Xstate that it has found one non-zero matrix.
  1832. X
  1833. XVarious versions of the make file (extension .mak) are included with the
  1834. Xpackage.
  1835. X
  1836. XThere is an optional #define DO_FREE_CHECK in include.h. If this is
  1837. Xactivated the program will check that the new-s and delete-s are
  1838. Xbalanced. There should be no unbalanced new-s or delete-s. You need the
  1839. XANSI version of the preprocessor to run DO_FREE_CHECK.
  1840. X
  1841. XThe program also allocates and deletes a large block and small block of
  1842. Xmemory before it starts the main testing and then at the end of the
  1843. Xtest. It then checks that the blocks of memory were allocated in the
  1844. Xsame place. If not then one suspects that there has been a memory leak.
  1845. Xi.e. a piece of memory has been allocated and not deleted.
  1846. X
  1847. XThis is not completely foolproof. Programs may allocate extra print
  1848. Xbuffers while the program is running. I have tried to overcome this by
  1849. Xdoing a print before I allocate the first memory block. Programs may
  1850. Xallocate memory for different sized items in different places, or might
  1851. Xnot allocate items consecutively. Or they might mix the items with memory
  1852. Xblocks from other programs. Nevertheless, I seem to get consistent
  1853. Xanswers from most of the compilers I am working with, so I think this is
  1854. Xa worthwhile test.
  1855. X
  1856. XGnu 2.2 and Zortech do not pass this test. There may be good reasons for
  1857. Xthis, but I think it would be a good idea if the authors of these
  1858. Xpackages made sure they knew what was happening.
  1859. X
  1860. END_OF_FILE
  1861. if test 1948 -ne `wc -c <'newmatc.txt'`; then
  1862.     echo shar: \"'newmatc.txt'\" unpacked with wrong size!
  1863. fi
  1864. # end of 'newmatc.txt'
  1865. fi
  1866. echo shar: End of archive 3 \(of 8\).
  1867. cp /dev/null ark3isdone
  1868. MISSING=""
  1869. for I in 1 2 3 4 5 6 7 8 ; do
  1870.     if test ! -f ark${I}isdone ; then
  1871.     MISSING="${MISSING} ${I}"
  1872.     fi
  1873. done
  1874. if test "${MISSING}" = "" ; then
  1875.     echo You have unpacked all 8 archives.
  1876.     rm -f ark[1-9]isdone
  1877. else
  1878.     echo You still need to unpack the following archives:
  1879.     echo "        " ${MISSING}
  1880. fi
  1881. ##  End of shell archive.
  1882. exit 0
  1883.  
  1884. exit 0 # Just in case...
  1885.