home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l391 / 2.ddi / MATH.BA$ / MATH.bin
Encoding:
Text File  |  1992-08-19  |  46.7 KB  |  1,203 lines
  1. ' ------------------------------------------------------------------------
  2. ' Visual Basic for MS-DOS Matrix Math ToolKit
  3. '
  4. ' The Matrix Math ToolKit (MATH.BAS and MATHA.BAS)
  5. ' contains routines which perform elementary operations
  6. ' on systems of linear equations represented as
  7. ' matrices.  The functions return integer error codes
  8. ' in the name and results in the parameter list.  The
  9. ' functions matbs?% and matlu?% found in this module
  10. ' are intended for internal use only.
  11. '
  12. ' Error codes returned:
  13. '     0  no error                     -1  matrix not invertible
  14. '    -2  matrix not square            -3  inner dimensions different
  15. '    -4  matrix dimensions different  -5  result matrix dimensioned incorrectly
  16. '    any other codes returned are standard BASIC errors
  17. '
  18. ' MATH.BAS is an 80x87 or emulator mathpack version of the
  19. ' toolkit and MATHA.BAS is an Alternate mathpack version.
  20. ' To use the Matrix Math routines in your program
  21. ' either include the appropriate source file in your
  22. ' program or use the supplied library (MATH.LIB, MATHA.LIB)
  23. ' and Quick Library (MATH.QLB) and call the appropriate
  24. ' procedures.
  25. '
  26. ' Copyright (C) 1982-1992 Microsoft Corporation
  27. '
  28. ' You have a royalty-free right to use, modify, reproduce
  29. ' and distribute the sample applications and toolkits provided with
  30. ' Visual Basic for MS-DOS (and/or any modified version)
  31. ' in any way you find useful, provided that you agree that
  32. ' Microsoft has no warranty, obligations or liability for
  33. ' any of the sample applications or toolkits.
  34. ' ------------------------------------------------------------------------
  35.  
  36. '===================================================================
  37. 'MatDet, MatSEqn, and MatInv all use LU-decomposition to implement Gaussian
  38. 'elimination.  A brief explanation of what is meant by an LU matrix is given
  39. 'below, followed by simplified versions of the two internal routines used to
  40. 'do all elimination.
  41. '
  42. 'What is meant by an LU matrix:
  43. 'An upper triangle matrix (one with all nonzero entries on or above the main
  44. 'diagonal) can be solved immediately.  The goal of Gaussian elimination is to
  45. 'transform a non upper triangle system into an equivalent triangular one.
  46. '
  47. 'Given a system of equations represented in matrix form by Ax=b, we need a
  48. 'linear transformation L such that LA=U where U is and upper triangular matrix.
  49. 'Then Ux=LAx=Lb and Ux=Lb is an upper triangular system.
  50. '
  51. 'This library explicitly calculates U, but L is never saved in its own array.
  52. 'When we do a row operation to create a zero below the main diagonal, we no
  53. 'longer need to save that value because we know it is zero.  This leaves the
  54. 'space available to save the multiplier used in the row operation.  When
  55. 'elimination is completed (ie, when the matrix is upper triangular), these
  56. 'multipliers give us a complete record of what we did to A to make it upper
  57. 'triangular.  This is equivalent to saying the multipliers represent L.  We now
  58. 'have a U and an L stored in the same matrix!  This type of matrix will be
  59. 'referred to as an LU matrix, or just LU.
  60. '
  61. 'The following code fragments get LU and backsolve Ux=Lb.  The actual routines
  62. 'used in the toolbox are much more involved because they implement total
  63. 'pivoting and implicit row scaling to reduce round off errors.  However, all the
  64. 'extras (pivoting, scaling, error checking) are extraneous to the main routines,
  65. 'which total only 20 lines.  If you are unfamilar with this type of matrix math,
  66. 'gaining an understanding of these 20 lines is a very good introduction.  Try
  67. 'working through a 2x2 or 3x3 example by hand to see what is happening.  The
  68. 'numerical techniques used to reduce round off error will not be discussed.
  69. '
  70. '-------------------------------------------------------------------
  71. 'Given the coefficient matrix A(1 TO N, 1 TO N) and the vector b(1 TO N),
  72. 'the following fragments will find x(1 TO N) satisfying Ax=b using Gaussian
  73. 'elimination.
  74. '
  75. 'matlu:
  76. 'Perform row operations to get all zeroes below the main diagonal.
  77. 'Define Rj(1 TO N) to be the vector corresponding to the jth row of A.
  78. 'Let Rrow = Rrow + m*Rpvt where m = -Rrow(pvt)/Rpvt(pvt).
  79. 'Then A(row, pvt)=0.
  80. '
  81. '** FOR pvt = 1 TO (N - 1)
  82. '**    FOR row = (pvt + 1) TO N
  83. '**       'Save m for later use in the space just made 0.
  84. '**       A(row, pvt) = -A(row, pvt) / A(pvt, pvt)
  85. '**       'Do the row operation.
  86. '**       FOR col = (pvt + 1) TO N
  87. '**          A(row, col) = A(row, col) + A(row, pvt) * A(pvt, col)
  88. '**       NEXT col
  89. '**    NEXT row
  90. '** NEXT pvt
  91. '
  92. 'matbs:
  93. 'Do the same row operations on b using the multipliers saved in A.
  94. '
  95. '** FOR pvt = 1 TO (N - 1)
  96. '**    FOR row = (pvt + 1) TO N
  97. '**       b(row) = b(row) + A(row, pvt) * b(pvt)
  98. '**    NEXT row
  99. '** NEXT pvt
  100. '
  101. 'Backsolve Ux=Lb to find x.
  102. '                               N
  103. 'For r = N to 1, x(r) = [b(r) - Σ (A(r,c)*x(c))]/A(r,r)
  104. '                              c=r+1
  105. '** FOR row = N TO 1 STEP -1
  106. '**    x(row) = b(row)
  107. '**    FOR col = (row + 1) TO N
  108. '**       x(row) = x(row) - A(row, col) * x(col)
  109. '**    NEXT col
  110. '**    x(row) = x(row) / A(row, row)
  111. '** NEXT row
  112. '
  113. '===================================================================
  114.  
  115. '$INCLUDE: 'math.bi'
  116. DECLARE FUNCTION matbsD% (A() AS DOUBLE, b() AS DOUBLE, x() AS DOUBLE)
  117. DECLARE FUNCTION matbsS% (A() AS SINGLE, b() AS SINGLE, x() AS SINGLE)
  118. DECLARE FUNCTION matluD% (A() AS DOUBLE)
  119. DECLARE FUNCTION matluS% (A() AS SINGLE)
  120. DIM SHARED lo AS INTEGER, up AS INTEGER
  121. DIM SHARED continue AS INTEGER, count AS INTEGER
  122. DIM SHARED rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  123. END
  124.  
  125. '=======================MatAddC%====================================
  126. 'MatAddC% adds two currency type matrices and places the sum in
  127. 'the first.
  128. '
  129. 'Parameters: matrices Alpha,Beta
  130. '
  131. 'Returns: Alpha() = Alpha() + Beta()
  132. '===================================================================
  133. FUNCTION MatAddC% (Alpha() AS CURRENCY, Beta() AS CURRENCY)
  134. ON LOCAL ERROR GOTO cadderr: MatAddC% = 0
  135. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  136. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  137. 'loop through and add elements
  138. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  139.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  140.       Alpha(row%, col%) = Alpha(row%, col%) + Beta(row%, col%)
  141.    NEXT col%
  142. NEXT row%
  143. caddexit:
  144. EXIT FUNCTION
  145. cadderr:
  146.    MatAddC% = (ERR + 5) MOD 200 - 5
  147.    RESUME caddexit
  148. END FUNCTION
  149.  
  150. '=======================MatAddD%====================================
  151. 'MatAddD% adds two double precision matrices and places the sum in
  152. 'the first.
  153. '
  154. 'Parameters: matrices Alpha,Beta
  155. '
  156. 'Returns: Alpha() = Alpha() + Beta()
  157. '===================================================================
  158. FUNCTION MatAddD% (Alpha() AS DOUBLE, Beta() AS DOUBLE)
  159. ON LOCAL ERROR GOTO dadderr: MatAddD% = 0
  160. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  161. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  162. 'loop through and add elements
  163. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  164.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  165.       Alpha(row%, col%) = Alpha(row%, col%) + Beta(row%, col%)
  166.    NEXT col%
  167. NEXT row%
  168. daddexit:
  169. EXIT FUNCTION
  170. dadderr:
  171.    MatAddD% = (ERR + 5) MOD 200 - 5
  172.    RESUME daddexit
  173. END FUNCTION
  174.  
  175. '=======================MatAddI%====================================
  176. 'MatAddI% adds two integer matrices and places the sum in
  177. 'the first.
  178. '
  179. 'Parameters: matrices Alpha,Beta
  180. '
  181. 'Returns: Alpha() = Alpha() + Beta()
  182. '===================================================================
  183. FUNCTION MatAddI% (Alpha() AS INTEGER, Beta() AS INTEGER)
  184. ON LOCAL ERROR GOTO iadderr: MatAddI% = 0
  185. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  186. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  187. 'loop through and add elements
  188. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  189.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  190.       Alpha(row%, col%) = Alpha(row%, col%) + Beta(row%, col%)
  191.    NEXT col%
  192. NEXT row%
  193. iaddexit:
  194. EXIT FUNCTION
  195. iadderr:
  196.    MatAddI% = (ERR + 5) MOD 200 - 5
  197.    RESUME iaddexit
  198. END FUNCTION
  199.  
  200. '=======================MatAddL%====================================
  201. 'MatAddL% adds two long integer matrices and places the sum in
  202. 'the first.
  203. '
  204. 'Parameters: matrices Alpha,Beta
  205. '
  206. 'Returns: Alpha() = Alpha() + Beta()
  207. '===================================================================
  208. FUNCTION MatAddL% (Alpha() AS LONG, Beta() AS LONG)
  209. ON LOCAL ERROR GOTO ladderr: MatAddL% = 0
  210. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  211. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  212. 'loop through and add elements
  213. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  214.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  215.       Alpha(row%, col%) = Alpha(row%, col%) + Beta(row%, col%)
  216.    NEXT col%
  217. NEXT row%
  218. laddexit:
  219. EXIT FUNCTION
  220. ladderr:
  221.    MatAddL% = (ERR + 5) MOD 200 - 5
  222.    RESUME laddexit
  223. END FUNCTION
  224.  
  225. '=======================MatAddS%====================================
  226. 'MatAddS% adds two single precision matrices and places the sum in
  227. 'the first.
  228. '
  229. 'Parameters: matrices Alpha,Beta
  230. '
  231. 'Returns: Alpha() = Alpha() + Beta()
  232. '===================================================================
  233. FUNCTION MatAddS% (Alpha() AS SINGLE, Beta() AS SINGLE)
  234. ON LOCAL ERROR GOTO sadderr: MatAddS% = 0
  235. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  236. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  237. 'loop through and add elements
  238. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  239.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  240.       Alpha(row%, col%) = Alpha(row%, col%) + Beta(row%, col%)
  241.    NEXT col%
  242. NEXT row%
  243. saddexit:
  244. EXIT FUNCTION
  245. sadderr:
  246.    MatAddS% = (ERR + 5) MOD 200 - 5
  247.    RESUME saddexit
  248. END FUNCTION
  249.  
  250. '========================matbsD=====================================
  251. 'matbsD% takes a matrix in LU form, found by matluD%, and a vector b
  252. 'and solves the system Ux=Lb for x. matrices A,b,x are double precision.
  253. '
  254. 'Parameters: LU matrix in A, corresponding pivot vectors in rpvt and cpvt,
  255. '            right side in b
  256. '
  257. 'Returns: solution in x, b is modified, rest unchanged
  258. '===================================================================
  259. FUNCTION matbsD% (A() AS DOUBLE, b() AS DOUBLE, x() AS DOUBLE)
  260. ON LOCAL ERROR GOTO dbserr: matbsD% = 0
  261. 'do row operations on b using the multipliers in L to find Lb
  262. FOR pvt% = lo TO (up - 1)
  263.    c% = cpvt(pvt%)
  264.    FOR row% = (pvt% + 1) TO up
  265.       r% = rpvt(row%)
  266.       b(r%) = b(r%) + A(r%, c%) * b(rpvt(pvt%))
  267.    NEXT row%
  268. NEXT pvt%
  269. 'backsolve Ux=Lb to find x
  270. FOR row% = up TO lo STEP -1
  271.    c% = cpvt(row%)
  272.    r% = rpvt(row%)
  273.    x(c%) = b(r%)
  274.    FOR col% = (row% + 1) TO up
  275.       x(c%) = x(c%) - A(r%, cpvt(col%)) * x(cpvt(col%))
  276.    NEXT col%
  277.    x(c%) = x(c%) / A(r%, c%)
  278. NEXT row%
  279. dbsexit:
  280. EXIT FUNCTION
  281. dbserr:
  282.    matbsD% = ERR
  283.    RESUME dbsexit
  284. END FUNCTION
  285.  
  286. '========================matbsS=====================================
  287. 'matbsS% takes a matrix in LU form, found by matluS%, and a vector b
  288. 'and solves the system Ux=Lb for x. matrices A,b,x are single precision.
  289. '
  290. 'Parameters: LU matrix in A, corresponding pivot vectors in rpvt and cpvt,
  291. '            right side in b
  292. '
  293. 'Returns: solution in x, b is modified, rest unchanged
  294. '===================================================================
  295. FUNCTION matbsS% (A() AS SINGLE, b() AS SINGLE, x() AS SINGLE)
  296. ON LOCAL ERROR GOTO sbserr: matbsS% = 0
  297. 'do row operations on b using the multipliers in L to find Lb
  298. FOR pvt% = lo TO (up - 1)
  299.    c% = cpvt(pvt%)
  300.    FOR row% = (pvt% + 1) TO up
  301.       r% = rpvt(row%)
  302.       b(r%) = b(r%) + A(r%, c%) * b(rpvt(pvt%))
  303.    NEXT row%
  304. NEXT pvt%
  305. 'backsolve Ux=Lb to find x
  306. FOR row% = up TO lo STEP -1
  307.    c% = cpvt(row%)
  308.    r% = rpvt(row%)
  309.    x(c%) = b(r%)
  310.    FOR col% = (row% + 1) TO up
  311.       x(c%) = x(c%) - A(r%, cpvt(col%)) * x(cpvt(col%))
  312.    NEXT col%
  313.    x(c%) = x(c%) / A(r%, c%)
  314. NEXT row%
  315. sbsexit:
  316. EXIT FUNCTION
  317. sbserr:
  318.    matbsS% = ERR
  319.    RESUME sbsexit
  320. END FUNCTION
  321.  
  322. '========================MatDetC%===================================
  323. 'MatDetC% finds the determinant of a square, currency type matrix
  324. '
  325. 'Parameters: A(n x n) matrix, det@ to return the determinant
  326. '
  327. 'Returns: matrix A in LU form, determinant
  328. '===================================================================
  329. FUNCTION MatDetC% (A() AS CURRENCY, det@)
  330. ON LOCAL ERROR GOTO cdeterr: errcode% = 0
  331. lo = LBOUND(A, 1)
  332. up = UBOUND(A, 1)
  333. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  334. 'make temporary double precision matrix to find pivots
  335. DIM Tmp(lo TO up, LBOUND(A, 2) TO UBOUND(A, 2)) AS DOUBLE
  336. FOR row% = lo TO up
  337.    FOR col% = LBOUND(A, 2) TO UBOUND(A, 2)
  338.       Tmp(row%, col%) = CDBL(A(row%, col%))
  339.    NEXT col%
  340. NEXT row%
  341. errcode% = matluD%(Tmp())              'Get LU matrix
  342. IF NOT continue THEN
  343.    IF errcode% = 199 THEN det@ = 0@
  344.    ERROR errcode%
  345. ELSE
  346.    detD# = 1#                          '+/- determinant = product of the pivots
  347.    FOR pvt% = lo TO up
  348.       detD# = detD# * Tmp(rpvt(pvt%), cpvt(pvt%))
  349.    NEXT pvt%                           'count contains the total number of row
  350.    det@ = (-1@) ^ count * CCUR(detD#)  'and column switches due to pivoting.
  351.    IF errcode% THEN ERROR errcode%     'multiply the determinant by -1 for
  352. END IF                                 'each switch.
  353. cdetexit:
  354. ERASE rpvt, cpvt, Tmp
  355. MatDetC% = errcode%
  356. EXIT FUNCTION
  357. cdeterr:
  358.    errcode% = (ERR + 5) MOD 200 - 5
  359.    RESUME cdetexit
  360. END FUNCTION
  361.  
  362. '========================MatDetD%===================================
  363. 'MatDetD% finds the determinant of a square, double precision matrix
  364. '
  365. 'Parameters: A(n x n) matrix, det# to return the determinant
  366. '
  367. 'Returns: matrix A in LU form, determinant
  368. '===================================================================
  369. FUNCTION MatDetD% (A() AS DOUBLE, det#)
  370. ON LOCAL ERROR GOTO ddeterr: errcode% = 0
  371. lo = LBOUND(A, 1)
  372. up = UBOUND(A, 1)
  373. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  374. errcode% = matluD%(A())             'Get LU matrix
  375. IF NOT continue THEN
  376.    IF errcode% = 199 THEN det# = 0#
  377.    ERROR errcode%
  378. ELSE
  379.    det# = 1#                        '+/- determinant = product of the pivots
  380.    FOR pvt% = lo TO up
  381.       det# = det# * A(rpvt(pvt%), cpvt(pvt%))
  382.    NEXT pvt%                         'count contains the total number of row
  383.    det# = (-1) ^ count * det#        'and column switches due to pivoting.
  384.    IF errcode% THEN ERROR errcode%   'multiply the determinant by -1 for
  385. END IF                               'each switch
  386. ddetexit:           
  387. ERASE rpvt, cpvt
  388. MatDetD% = errcode%
  389. EXIT FUNCTION
  390. ddeterr:
  391.    errcode% = (ERR + 5) MOD 200 - 5
  392.    RESUME ddetexit
  393. END FUNCTION
  394.  
  395. '========================MatDetI%===================================
  396. 'MatDetI% finds the determinant of a square, integer matrix
  397. '
  398. 'Parameters: A(n x n) matrix, det% to return the determinant
  399. '
  400. 'Returns: matrix A unchanged, determinant
  401. '===================================================================
  402. FUNCTION MatDetI% (A() AS INTEGER, det%)
  403. ON LOCAL ERROR GOTO ideterr: errcode% = 0
  404. lo = LBOUND(A, 1)
  405. up = UBOUND(A, 1)
  406. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  407. 'make temporary single precision matrix to find pivots
  408. DIM Tmp(lo TO up, LBOUND(A, 2) TO UBOUND(A, 2)) AS SINGLE
  409. FOR row% = lo TO up
  410.    FOR col% = LBOUND(A, 2) TO UBOUND(A, 2)
  411.       Tmp(row%, col%) = CSNG(A(row%, col%))
  412.    NEXT col%
  413. NEXT row%
  414. errcode% = matluS%(Tmp())              'Get LU matrix
  415. IF NOT continue THEN
  416.    IF errcode% = 199 THEN det% = 0
  417.    ERROR errcode%
  418. ELSE
  419.    detS! = 1!                          '+/- determinant = product of the pivots
  420.    FOR pvt% = lo TO up
  421.       detS! = detS! * Tmp(rpvt(pvt%), cpvt(pvt%))
  422.    NEXT pvt%                           'count contains the total number of row
  423.    det% = (-1) ^ count * CINT(detS!)   'and column switches due to pivoting.
  424.    IF errcode% THEN ERROR errcode%     'multiply the determinant by -1 for
  425. END IF                                 'each switch
  426. idetexit:
  427. ERASE rpvt, cpvt, Tmp
  428. MatDetI% = errcode%
  429. EXIT FUNCTION
  430. ideterr:
  431.    errcode% = (ERR + 5) MOD 200 - 5
  432.    RESUME idetexit
  433. END FUNCTION
  434.  
  435. '========================MatDetL%===================================
  436. 'MatDetL% finds the determinant of a square, long integer matrix
  437. '
  438. 'Parameters: A(n x n) matrix, det& to return the determinant
  439. '
  440. 'Returns: matrix A unchanged, determinant
  441. '===================================================================
  442. FUNCTION MatDetL% (A() AS LONG, det&)
  443. ON LOCAL ERROR GOTO ldeterr: errcode% = 0
  444. lo = LBOUND(A, 1)
  445. up = UBOUND(A, 1)
  446. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  447. 'make temporary double precision matrix to find pivots
  448. DIM Tmp(lo TO up, LBOUND(A, 2) TO UBOUND(A, 2)) AS DOUBLE
  449. FOR row% = lo TO up
  450.    FOR col% = LBOUND(A, 2) TO UBOUND(A, 2)
  451.       Tmp(row%, col%) = CDBL(A(row%, col%))
  452.    NEXT col%
  453. NEXT row%
  454. errcode% = matluD%(Tmp())              'Get LU matrix
  455. IF NOT continue THEN
  456.    IF errcode% = 199 THEN det& = 0&
  457.    ERROR errcode%
  458. ELSE
  459.    detD# = 1#                          '+/- determinant = product of the pivots
  460.    FOR pvt% = lo TO up
  461.       detD# = detD# * Tmp(rpvt(pvt%), cpvt(pvt%))
  462.    NEXT pvt%                           'count contains the total number of row
  463.    det& = (-1&) ^ count * CLNG(detD#)  'and column switches due to pivoting.
  464.    IF errcode% THEN ERROR errcode%     'multiply the determinant by -1 for
  465. END IF                                 'each switch
  466. ldetexit:
  467. ERASE rpvt, cpvt, Tmp
  468. MatDetL% = errcode%
  469. EXIT FUNCTION
  470. ldeterr:
  471.    errcode% = (ERR + 5) MOD 200 - 5
  472.    RESUME ldetexit
  473. END FUNCTION
  474.  
  475. '========================MatDetS%===================================
  476. 'MatDetS% finds the determinant of a square, single precision matrix
  477. '
  478. 'Parameters: A(n x n) matrix, det! to return the determinant
  479. '
  480. 'Returns: matrix A in LU form, determinant
  481. '===================================================================
  482. FUNCTION MatDetS% (A() AS SINGLE, det!)
  483. ON LOCAL ERROR GOTO sdeterr: errcode% = 0
  484. lo = LBOUND(A, 1)
  485. up = UBOUND(A, 1)
  486. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  487. errcode% = matluS%(A())                'Get LU matrix
  488. IF NOT continue THEN
  489.    IF errcode% = 199 THEN det! = 0!
  490.    ERROR errcode%
  491. ELSE
  492.    det! = 1!                           '+/- determinant = product of the pivots
  493.    FOR pvt% = lo TO up
  494.       det! = det! * A(rpvt(pvt%), cpvt(pvt%))
  495.    NEXT pvt%                           'count contains the total number of row
  496.    det! = (-1) ^ count * det!          'and column switches due to pivoting.
  497.    IF errcode% THEN ERROR errcode%     'multiply the determinant by -1 for
  498. END IF                                 'each switch
  499. sdetexit:
  500. ERASE rpvt, cpvt
  501. MatDetS% = errcode%
  502. EXIT FUNCTION
  503. sdeterr:
  504.    errcode% = (ERR + 5) MOD 200 - 5
  505.    RESUME sdetexit
  506. END FUNCTION
  507.  
  508. '========================MatInvC%===================================
  509. 'MatInvC% uses the matluD% and matbsD procedures to invert a square, currency
  510. 'type matrix.  Let e(N) contain all zeroes except for the jth position, which
  511. 'is 1.  Then the jth column of A^-1 is x, where Ax=e.
  512. '
  513. 'Parameters: A(n x n) matrix
  514. '
  515. 'Returns: A^-1
  516. '===================================================================
  517. FUNCTION MatInvC% (A() AS CURRENCY)
  518. ON LOCAL ERROR GOTO cinverr: errcode% = 0
  519. lo = LBOUND(A, 1)
  520. up = UBOUND(A, 1)
  521. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  522. 'duplicate A() in a double precision work matrix, Tmp()
  523. DIM Tmp(lo TO up, LBOUND(A, 2) TO UBOUND(A, 2)) AS DOUBLE
  524. DIM e(lo TO up) AS DOUBLE, x(lo TO up) AS DOUBLE
  525. FOR row% = lo TO up
  526.    FOR col% = LBOUND(A, 2) TO UBOUND(A, 2)
  527.       Tmp(row%, col%) = CDBL(A(row%, col%))
  528.    NEXT col%
  529. NEXT row%
  530. errcode% = matluD%(Tmp())                    'Put LU in Tmp
  531. IF NOT continue THEN ERROR errcode%
  532. FOR col% = lo TO up                          'Find A^-1 one column at a time
  533.    e(col%) = 1#
  534.    bserrcode% = matbsD%(Tmp(), e(), x())
  535.    IF bserrcode% THEN ERROR bserrcode%
  536.    FOR row% = lo TO up
  537.       A(row%, col%) = CCUR(x(row%))          'Put the column into A
  538.       e(row%) = 0#
  539.    NEXT row%
  540. NEXT col%
  541. IF errcode% THEN ERROR errcode%
  542. cinvexit:
  543. ERASE Tmp, e, x, rpvt, cpvt
  544. MatInvC% = errcode%
  545. EXIT FUNCTION
  546. cinverr:
  547.    errcode% = (ERR + 5) MOD 200 - 5
  548.    RESUME cinvexit
  549. END FUNCTION
  550.  
  551. '========================MatInvD%===================================
  552. 'MatInvD% uses the matluD% and matbsD procedures to invert a square, double
  553. 'precision matrix.  Let e(N) contain all zeroes except for the jth position,
  554. 'which is 1.  Then the jth column of A^-1 is x, where Ax=e.
  555. '
  556. 'Parameters: A(n x n) matrix
  557. '
  558. 'Returns: A^-1
  559. '===================================================================
  560. FUNCTION MatInvD% (A() AS DOUBLE)
  561. ON LOCAL ERROR GOTO dinverr: errcode% = 0
  562. lo = LBOUND(A, 1)
  563. up = UBOUND(A, 1)
  564. DIM Ain(lo TO up, lo TO up) AS DOUBLE
  565. DIM e(lo TO up) AS DOUBLE, x(lo TO up) AS DOUBLE
  566. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  567. errcode% = matluD%(A())                     'Get LU matrix
  568. IF NOT continue THEN ERROR errcode%
  569. FOR col% = lo TO up                         'Find A^-1 one column at a time
  570.    e(col%) = 1#
  571.    bserrcode% = matbsD%(A(), e(), x())
  572.    IF bserrcode% THEN ERROR bserrcode%
  573.    FOR row% = lo TO up
  574.       Ain(row%, col%) = x(row%)
  575.       e(row%) = 0#
  576.    NEXT row%
  577. NEXT col%
  578. FOR col% = lo TO up                         'Put A^-1 in A
  579.    FOR row% = lo TO up
  580.       A(row%, col%) = Ain(row%, col%)
  581.    NEXT row%
  582. NEXT col%
  583. IF errcode% THEN ERROR errcode%
  584. dinvexit:
  585. ERASE e, x, Ain, rpvt, cpvt
  586. MatInvD% = errcode%
  587. EXIT FUNCTION
  588. dinverr:
  589.    errcode% = (ERR + 5) MOD 200 - 5
  590.    RESUME dinvexit
  591. END FUNCTION
  592.  
  593. '========================MatInvS%===================================
  594. 'MatInvS% uses the matluS% and matbsS procedures to invert a square, single
  595. 'precision matrix.  Let e(N) contain all zeroes except for the jth position,
  596. 'which is 1. Then the jth column of A^-1 is x, where Ax=e.
  597. '
  598. 'Parameters: A(n x n) matrix
  599. '
  600. 'Returns: A^-1
  601. '===================================================================
  602. FUNCTION MatInvS% (A() AS SINGLE)
  603. ON LOCAL ERROR GOTO sinverr: errcode% = 0
  604. lo = LBOUND(A, 1)
  605. up = UBOUND(A, 1)
  606. DIM Ain(lo TO up, lo TO up) AS SINGLE
  607. DIM e(lo TO up) AS SINGLE, x(lo TO up) AS SINGLE
  608. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  609. errcode% = matluS%(A())                     'Get LU matrix
  610. IF NOT continue THEN ERROR errcode%
  611. FOR col% = lo TO up                         'find A^-1 one column at a time
  612.    e(col%) = 1!
  613.    bserrcode% = matbsS%(A(), e(), x())
  614.    IF bserrcode% THEN ERROR bserrcode%
  615.    FOR row% = lo TO up
  616.       Ain(row%, col%) = x(row%)
  617.       e(row%) = 0!
  618.    NEXT row%
  619. NEXT col%
  620. FOR col% = lo TO up                         'put A^-1 in A
  621.    FOR row% = lo TO up
  622.       A(row%, col%) = Ain(row%, col%)
  623.    NEXT row%
  624. NEXT col%
  625. IF errcode% THEN ERROR errcode%
  626. sinvexit:
  627. ERASE e, x, Ain, rpvt, cpvt
  628. MatInvS% = errcode%
  629. EXIT FUNCTION
  630. sinverr:
  631.    errcode% = (ERR + 5) MOD 200 - 5
  632.    RESUME sinvexit
  633. END FUNCTION
  634.  
  635. '========================matluD%====================================
  636. 'matluD% does Gaussian elimination with total pivoting to put a square, double
  637. 'precision matrix in LU form. The multipliers used in the row operations to
  638. 'create zeroes below the main diagonal are saved in the zero spaces.
  639. '
  640. 'Parameters: A(n x n) matrix, rpvt(n) and cpvt(n) permutation vectors
  641. '            used to index the row and column pivots
  642. '
  643. 'Returns: A in LU form with corresponding pivot vectors; the total number of
  644. '         pivots in count, which is used to find the sign of the determinant.
  645. '===================================================================
  646. FUNCTION matluD% (A() AS DOUBLE)
  647. ON LOCAL ERROR GOTO dluerr: errcode% = 0
  648. 'Checks if A is square, returns error code if not
  649. IF NOT (lo = LBOUND(A, 2) AND up = UBOUND(A, 2)) THEN ERROR 198
  650. DIM rownorm(lo TO up) AS DOUBLE
  651. count = 0                            'initialize count, continue
  652. continue = -1
  653. FOR row% = lo TO up                  'initialize rpvt and cpvt
  654.    rpvt(row%) = row%
  655.    cpvt(row%) = row%
  656.    rownorm(row%) = 0#                'find the row norms of A()
  657.    FOR col% = lo TO up
  658.       rownorm(row%) = rownorm(row%) + ABS(A(row%, col%))
  659.    NEXT col%
  660.    IF rownorm(row%) = 0# THEN        'if any rownorm is zero, the matrix
  661.       continue = 0                   'is singular, set error, exit and
  662.       ERROR 199                      'do not continue
  663.    END IF
  664. NEXT row%
  665. FOR pvt% = lo TO (up - 1)
  666. 'Find best available pivot
  667.    max# = 0#                         'checks all values in rows and columns not
  668.    FOR row% = pvt% TO up             'already used for pivoting and saves the
  669.       r% = rpvt(row%)                'largest absolute number and its position
  670.       FOR col% = pvt% TO up
  671.          c% = cpvt(col%)
  672.          temp# = ABS(A(r%, c%)) / rownorm(r%)
  673.          IF temp# > max# THEN
  674.             max# = temp#
  675.             bestrow% = row%          'save the position of new max#
  676.             bestcol% = col%
  677.          END IF
  678.       NEXT col%
  679.    NEXT row%
  680.    IF max# = 0# THEN                 'if no nonzero number is found, A is
  681.       continue = 0                   'singular, send back error, do not continue
  682.       ERROR 199
  683.    ELSEIF pvt% > 1 THEN              'check if drop in pivots is too much
  684.       IF max# < (deps# * oldmax#) THEN errcode% = 199
  685.    END IF
  686.    oldmax# = max#
  687.    IF rpvt(pvt%) <> rpvt(bestrow%) THEN
  688.       count = count + 1                    'if a row or column pivot is
  689.       SWAP rpvt(pvt%), rpvt(bestrow%)      'necessary, count it and permute
  690.    END IF                                  'rpvt or cpvt. Note: the rows and
  691.    IF cpvt(pvt%) <> cpvt(bestcol%) THEN    'columns are not actually switched,
  692.       count = count + 1                    'only the order in which they are
  693.       SWAP cpvt(pvt%), cpvt(bestcol%)      'used.
  694.    END IF
  695. 'Eliminate all values below the pivot
  696.    rp% = rpvt(pvt%)
  697.    cp% = cpvt(pvt%)
  698.    FOR row% = (pvt% + 1) TO up
  699.       r% = rpvt(row%)
  700.       A(r%, cp%) = -A(r%, cp%) / A(rp%, cp%)  'save multipliers
  701.       FOR col% = (pvt% + 1) TO up
  702.          c% = cpvt(col%)                      'complete row operations
  703.          A(r%, c%) = A(r%, c%) + A(r%, cp%) * A(rp%, c%)
  704.       NEXT col%
  705.    NEXT row%
  706. NEXT pvt%
  707. IF A(rpvt(up), cpvt(up)) = 0# THEN
  708.    continue = 0                      'if last pivot is zero or pivot drop is
  709.    ERROR 199                         'too large, A is singular, send back error
  710. ELSEIF (ABS(A(rpvt(up), cpvt(up))) / rownorm(rpvt(up))) < (deps# * oldmax#) THEN
  711.    errcode% = 199                    'if pivot is not identically zero then
  712. END IF                               'continue remains TRUE
  713. IF errcode% THEN ERROR errcode%
  714. dluexit:
  715. matluD% = errcode%
  716. EXIT FUNCTION
  717. dluerr:
  718.    IF errcode% < 199 THEN continue = 0
  719.    errcode% = ERR
  720.    RESUME dluexit
  721. END FUNCTION
  722.  
  723. '========================matluS%====================================
  724. 'matluS% does Gaussian elimination with total pivoting to put a square, single
  725. 'precision matrix in LU form. The multipliers used in the row operations to
  726. 'create zeroes below the main diagonal are saved in the zero spaces.
  727. '
  728. 'Parameters: A(n x n) matrix, rpvt(n) and cpvt(n) permutation vectors
  729. '            used to index the row and column pivots
  730. '
  731. 'Returns: A in LU form with corresponding pivot vectors; the total number of
  732. '         pivots in count, which is used to find the sign of the determinant.
  733. '===================================================================
  734. FUNCTION matluS% (A() AS SINGLE)
  735. ON LOCAL ERROR GOTO sluerr: errcode% = 0
  736. 'Checks if A is square, returns error code if not
  737. IF NOT (lo = LBOUND(A, 2) AND up = UBOUND(A, 2)) THEN ERROR 198
  738. DIM rownorm(lo TO up) AS SINGLE
  739. count = 0                            'initialize count, continue
  740. continue = -1
  741. FOR row% = lo TO up                  'initialize rpvt and cpvt
  742.    rpvt(row%) = row%
  743.    cpvt(row%) = row%
  744.    rownorm(row%) = 0!                'find the row norms of A()
  745.    FOR col% = lo TO up
  746.       rownorm(row%) = rownorm(row%) + ABS(A(row%, col%))
  747.    NEXT col%
  748.    IF rownorm(row%) = 0! THEN        'if any rownorm is zero, the matrix
  749.       continue = 0                   'is singular, set error, exit and do
  750.       ERROR 199                      'not continue
  751.    END IF
  752. NEXT row%
  753. FOR pvt% = lo TO (up - 1)
  754. 'Find best available pivot
  755.    max! = 0!                         'checks all values in rows and columns not
  756.    FOR row% = pvt% TO up             'already used for pivoting and finds the
  757.       r% = rpvt(row%)                'number largest in absolute value relative
  758.       FOR col% = pvt% TO up          'to its row norm
  759.          c% = cpvt(col%)
  760.          temp! = ABS(A(r%, c%)) / rownorm(r%)
  761.          IF temp! > max! THEN
  762.             max! = temp!
  763.             bestrow% = row%          'save the position of new max!
  764.             bestcol% = col%
  765.          END IF
  766.       NEXT col%
  767.    NEXT row%
  768.    IF max! = 0! THEN                 'if no nonzero number is found, A is
  769.       continue = 0                   'singular, send back error, do not continue
  770.       ERROR 199
  771.    ELSEIF pvt% > 1 THEN              'check if drop in pivots is too much
  772.       IF max! < (seps! * oldmax!) THEN errcode% = 199
  773.    END IF
  774.    oldmax! = max!
  775.    IF rpvt(pvt%) <> rpvt(bestrow%) THEN
  776.       count = count + 1                    'if a row or column pivot is
  777.       SWAP rpvt(pvt%), rpvt(bestrow%)      'necessary, count it and permute
  778.    END IF                                  'rpvt or cpvt. Note: the rows and
  779.    IF cpvt(pvt%) <> cpvt(bestcol%) THEN    'columns are not actually switched,
  780.       count = count + 1                    'only the order in which they are
  781.       SWAP cpvt(pvt%), cpvt(bestcol%)      'used.
  782.    END IF
  783. 'Eliminate all values below the pivot
  784.    rp% = rpvt(pvt%)
  785.    cp% = cpvt(pvt%)
  786.    FOR row% = (pvt% + 1) TO up
  787.       r% = rpvt(row%)
  788.       A(r%, cp%) = -A(r%, cp%) / A(rp%, cp%)  'save multipliers
  789.       FOR col% = (pvt% + 1) TO up
  790.          c% = cpvt(col%)                      'complete row operations
  791.          A(r%, c%) = A(r%, c%) + A(r%, cp%) * A(rp%, c%)
  792.       NEXT col%
  793.    NEXT row%
  794. NEXT pvt%
  795. IF A(rpvt(up), cpvt(up)) = 0! THEN
  796.    continue = 0                      'if last pivot is zero or pivot drop is
  797.    ERROR 199                         'too large, A is singular, send back error
  798. ELSEIF (ABS(A(rpvt(up), cpvt(up))) / rownorm(rpvt(up))) < (seps! * oldmax!) THEN
  799.    errcode% = 199                    'if pivot is not identically zero then
  800. END IF                               'continue remains TRUE
  801. IF errcode% THEN ERROR errcode%
  802. sluexit:
  803. matluS% = errcode%
  804. EXIT FUNCTION
  805. sluerr:
  806.    errcode% = ERR
  807.    IF errcode% < 199 THEN continue = 0
  808.    RESUME sluexit
  809. END FUNCTION
  810.  
  811. '=======================MatMultC%===================================
  812. 'MatMultC% multiplies two currency type matrices and places the
  813. 'product in a result matrix
  814. '
  815. 'Parameters: matrices Alpha,Beta,Gamma
  816. '
  817. 'Returns: Gamma() = Alpha() * Beta()
  818. '===================================================================
  819. FUNCTION MatMultC% (Alpha() AS CURRENCY, Beta() AS CURRENCY, Gamma() AS CURRENCY)
  820. ON LOCAL ERROR GOTO cmulterr: MatMultC% = 0
  821. IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN
  822.    ERROR 197                   'check inside dimensions
  823. ELSEIF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) OR (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) OR (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN
  824.    ERROR 195                   'check dimensions of result matrix
  825. END IF
  826. 'loop through, Gamma(row,col)=inner product of Alpha(row,*) and Beta(*,col)
  827. FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
  828.    FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
  829.       Gamma(row%, col%) = 0@
  830.       FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  831.          Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
  832.       NEXT inside%
  833.    NEXT col%
  834. NEXT row%
  835. cmultexit:
  836. EXIT FUNCTION
  837. cmulterr:
  838.    MatMultC% = (ERR + 5) MOD 200 - 5
  839.    RESUME cmultexit
  840. END FUNCTION
  841.  
  842. '=======================MatMultD%===================================
  843. 'MatMultD% multiplies two double precision matrices and places the
  844. 'product in a result matrix
  845. '
  846. 'Parameters: matrices Alpha,Beta,Gamma
  847. '
  848. 'Returns: Gamma() = Alpha() * Beta()
  849. '===================================================================
  850. FUNCTION MatMultD% (Alpha() AS DOUBLE, Beta() AS DOUBLE, Gamma() AS DOUBLE)
  851. ON LOCAL ERROR GOTO dmulterr: MatMultD% = 0
  852. IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN
  853.    ERROR 197                   'check inside dimensions
  854. ELSEIF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) OR (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) OR (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN
  855.    ERROR 195                   'check dimensions of result matrix
  856. END IF
  857. 'loop through, Gamma(row,col)=inner product of Alpha(row,*) and Beta(*,col)
  858. FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
  859.    FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
  860.       Gamma(row%, col%) = 0#
  861.       FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  862.          Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
  863.       NEXT inside%
  864.    NEXT col%
  865. NEXT row%
  866. dmultexit:
  867. EXIT FUNCTION
  868. dmulterr:
  869.    MatMultD% = (ERR + 5) MOD 200 - 5
  870.    RESUME dmultexit
  871. END FUNCTION
  872.  
  873. '=======================MatMultI%===================================
  874. 'MatMultI% multiplies two integer matrices and places the product in
  875. 'a result matrix
  876. '
  877. 'Parameters: matrices Alpha,Beta,Gamma
  878. '
  879. 'Returns: Gamma() = Alpha() * Beta()
  880. '===================================================================
  881. FUNCTION MatMultI% (Alpha() AS INTEGER, Beta() AS INTEGER, Gamma() AS INTEGER)
  882. ON LOCAL ERROR GOTO imulterr: MatMultI% = 0
  883. IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN
  884.    ERROR 197                   'check inside dimensions
  885. ELSEIF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) OR (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) OR (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN
  886.    ERROR 195                   'check dimensions of result matrix
  887. END IF
  888. 'loop through, Gamma(row,col)=inner product of Alpha(row,*) and Beta(*,col)
  889. FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
  890.    FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
  891.       Gamma(row%, col%) = 0
  892.       FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  893.          Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
  894.       NEXT inside%
  895.    NEXT col%
  896. NEXT row%
  897. imultexit:
  898. EXIT FUNCTION
  899. imulterr:
  900.    MatMultI% = (ERR + 5) MOD 200 - 5
  901.    RESUME imultexit
  902. END FUNCTION
  903.  
  904. '=======================MatMultL%===================================
  905. 'MatMultL% multiplies two long integer matrices and places the product
  906. 'in a result matrix
  907. '
  908. 'Parameters: matrices Alpha,Beta,Gamma
  909. '
  910. 'Returns: Gamma() = Alpha() * Beta()
  911. '===================================================================
  912. FUNCTION MatMultL% (Alpha() AS LONG, Beta() AS LONG, Gamma() AS LONG)
  913. ON LOCAL ERROR GOTO lmulterr: MatMultL% = 0
  914. IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN
  915.    ERROR 197                   'check inside dimensions
  916. ELSEIF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) OR (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) OR (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN
  917.    ERROR 195                   'check dimensions of result matrix
  918. END IF
  919. 'loop through, Gamma(row,col)=inner product of Alpha(row,*) and Beta(*,col)
  920. FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
  921.    FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
  922.       Gamma(row%, col%) = 0&
  923.       FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  924.          Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
  925.       NEXT inside%
  926.    NEXT col%
  927. NEXT row%
  928. lmultexit:
  929. EXIT FUNCTION
  930. lmulterr:
  931.    MatMultL% = (ERR + 5) MOD 200 - 5
  932.    RESUME lmultexit
  933. END FUNCTION
  934.  
  935. '=======================MatMultS%===================================
  936. 'MatMultS% multiplies two single precision matrices and places the
  937. 'product in a result matrix
  938. '
  939. 'Parameters: matrices Alpha,Beta,Gamma
  940. '
  941. 'Returns: Gamma() = Alpha() * Beta()
  942. '===================================================================
  943. FUNCTION MatMultS% (Alpha() AS SINGLE, Beta() AS SINGLE, Gamma() AS SINGLE)
  944. ON LOCAL ERROR GOTO smulterr: MatMultS% = 0
  945. IF (LBOUND(Alpha, 2) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 1)) THEN
  946.    ERROR 197                   'check inside dimensions
  947. ELSEIF (LBOUND(Alpha, 1) <> LBOUND(Gamma, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Gamma, 1)) OR (LBOUND(Beta, 2) <> LBOUND(Gamma, 2)) OR (UBOUND(Beta, 2) <> UBOUND(Gamma, 2)) THEN
  948.    ERROR 195                   'check dimensions of result matrix
  949. END IF
  950. 'loop through, Gamma(row,col)=inner product of Alpha(row,*) and Beta(*,col)
  951. FOR row% = LBOUND(Gamma, 1) TO UBOUND(Gamma, 1)
  952.    FOR col% = LBOUND(Gamma, 2) TO UBOUND(Gamma, 2)
  953.       Gamma(row%, col%) = 0!
  954.       FOR inside% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  955.          Gamma(row%, col%) = Gamma(row%, col%) + Alpha(row%, inside%) * Beta(inside%, col%)
  956.       NEXT inside%
  957.    NEXT col%
  958. NEXT row%
  959. smultexit:
  960. EXIT FUNCTION
  961. smulterr:
  962.    MatMultS% = (ERR + 5) MOD 200 - 5
  963.    RESUME smultexit
  964. END FUNCTION
  965.  
  966. '========================MatSEqnC%==================================
  967. 'MatSEqnC% solves a system of n linear equations, Ax=b, and puts the
  968. 'answer in b. A is first put in LU form by matluC%, then matbsC is called
  969. 'to solve the system.  matrices A,b are currency type.
  970. '
  971. 'Parameters: A(n x n) contains coefficient matrix, b(N) contains the right side
  972. '
  973. 'Returns: A in LU form, solution in b
  974. '===================================================================
  975. FUNCTION MatSEqnC% (A() AS CURRENCY, b() AS CURRENCY)
  976. ON LOCAL ERROR GOTO cseqnerr: errcode% = 0
  977. lo = LBOUND(A, 1)
  978. up = UBOUND(A, 1)
  979. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  980. 'duplicate A(), b() in temporary double precision matrices Tmp(), btmp()
  981. DIM Tmp(lo TO up, LBOUND(A, 2) TO UBOUND(A, 2)) AS DOUBLE
  982. DIM x(lo TO up) AS DOUBLE, btmp(lo TO up) AS DOUBLE
  983. FOR row% = lo TO up
  984.    FOR col% = LBOUND(A, 2) TO UBOUND(A, 2)
  985.       Tmp(row%, col%) = CDBL(A(row%, col%))
  986.    NEXT col%
  987. NEXT row%
  988. errcode% = matluD%(Tmp())                   'Get LU matrix
  989. IF NOT continue THEN ERROR errcode%
  990. 'check dimensions of b, make double precision copy if ok.
  991. IF (lo <> LBOUND(b)) OR (up <> UBOUND(b)) THEN ERROR 197
  992. FOR row% = lo TO up
  993.    btmp(row%) = CDBL(b(row%))
  994. NEXT row%
  995. bserrcode% = matbsD%(Tmp(), btmp(), x())    'Backsolve system
  996. IF bserrcode% THEN ERROR bserrcode%
  997. FOR row% = lo TO up
  998.    b(row%) = CCUR(x(row%))                  'Put solution in b for return
  999. NEXT row%
  1000. IF errcode% THEN ERROR errcode%
  1001. cseqnexit:
  1002. ERASE Tmp, btmp, x, rpvt, cpvt
  1003. MatSEqnC% = errcode%
  1004. EXIT FUNCTION
  1005. cseqnerr:
  1006.    errcode% = (ERR + 5) MOD 200 - 5
  1007.    RESUME cseqnexit
  1008. END FUNCTION
  1009.  
  1010. '========================MatSEqnD%==================================
  1011. 'MatSEqnD% solves a system of n linear equations, Ax=b, and puts the
  1012. 'answer in b. A is first put in LU form by matluD%, then matbsD is called
  1013. 'to solve the system.  matrices A,b are double precision.
  1014. '
  1015. 'Parameters: A(n x n) contains coefficient matrix, b(N) contains the right side
  1016. '
  1017. 'Returns: A in LU form, solution in b
  1018. '===================================================================
  1019. FUNCTION MatSEqnD% (A() AS DOUBLE, b() AS DOUBLE)
  1020. ON LOCAL ERROR GOTO dseqnerr: errcode% = 0
  1021. lo = LBOUND(A, 1)
  1022. up = UBOUND(A, 1)
  1023. DIM x(lo TO up) AS DOUBLE
  1024. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  1025. errcode% = matluD%(A())                      'Get LU matrix
  1026. IF NOT continue THEN ERROR errcode%
  1027. 'check dimensions of b
  1028. IF (lo <> LBOUND(b)) OR (up <> UBOUND(b)) THEN ERROR 197
  1029. bserrcode% = matbsD%(A(), b(), x())          'Backsolve system
  1030. IF bserrcode% THEN ERROR bserrcode%
  1031. FOR row% = lo TO up
  1032.    b(row%) = x(row%)                         'Put solution in b for return
  1033. NEXT row%
  1034. IF errcode% THEN ERROR errcode%
  1035. dseqnexit:
  1036. ERASE x, rpvt, cpvt
  1037. MatSEqnD% = errcode%
  1038. EXIT FUNCTION
  1039. dseqnerr:
  1040.    errcode% = (ERR + 5) MOD 200 - 5
  1041.    RESUME dseqnexit
  1042. END FUNCTION
  1043.  
  1044. '========================MatSEqnS%==================================
  1045. 'MatSEqnS% solves a system of n linear equations, Ax=b, and puts the
  1046. 'answer in b. A is first put in LU form by matluS%, then matbsS is called
  1047. 'to solve the system.  matrices A,b are single precision.
  1048. '
  1049. 'Parameters: A(n x n) contains coefficient matrix, b(N) contains the right side
  1050. '
  1051. 'Returns: A in LU form, solution in b
  1052. '===================================================================
  1053. FUNCTION MatSEqnS% (A() AS SINGLE, b() AS SINGLE)
  1054. ON LOCAL ERROR GOTO sseqnerr: errcode% = 0
  1055. lo = LBOUND(A, 1)
  1056. up = UBOUND(A, 1)
  1057. DIM x(lo TO up) AS SINGLE
  1058. REDIM rpvt(lo TO up) AS INTEGER, cpvt(lo TO up) AS INTEGER
  1059. errcode% = matluS%(A())                      'Get LU matrix
  1060. IF NOT continue THEN ERROR errcode%
  1061. 'check dimensions of b
  1062. IF (lo <> LBOUND(b)) OR (up <> UBOUND(b)) THEN ERROR 197
  1063. bserrcode% = matbsS%(A(), b(), x())          'Backsolve system
  1064. IF bserrcode% THEN ERROR bserrcode%
  1065. FOR row% = lo TO up
  1066.    b(row%) = x(row%)                         'Put solution in b for return
  1067. NEXT row%
  1068. IF errcode% THEN ERROR errcode%
  1069. sseqnexit:
  1070. ERASE x, rpvt, cpvt
  1071. MatSEqnS% = errcode%
  1072. EXIT FUNCTION
  1073. sseqnerr:
  1074.    errcode% = (ERR + 5) MOD 200 - 5
  1075.    RESUME sseqnexit
  1076. END FUNCTION
  1077.  
  1078. '=======================MatSubC%====================================
  1079. 'MatSubC% takes the difference of two currency type matrices and
  1080. 'places the result in the first.
  1081. '
  1082. 'Params: matrices Alpha,Beta
  1083. '
  1084. 'Returns: Alpha=Alpha-Beta
  1085. '===================================================================
  1086. FUNCTION MatSubC% (Alpha() AS CURRENCY, Beta() AS CURRENCY)
  1087. ON LOCAL ERROR GOTO csuberr: MatSubC% = 0
  1088. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  1089. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  1090. 'loop through and subtract elements
  1091. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  1092.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  1093.       Alpha(row%, col%) = Alpha(row%, col%) - Beta(row%, col%)
  1094.    NEXT col%
  1095. NEXT row%
  1096. csubexit:
  1097. EXIT FUNCTION
  1098. csuberr:
  1099.    MatSubC% = (ERR + 5) MOD 200 - 5
  1100.    RESUME csubexit:
  1101. END FUNCTION
  1102.  
  1103. '=======================MatSubD%====================================
  1104. 'MatSubD% takes the difference of two double precision matrices and
  1105. 'places the result in the first.
  1106. '
  1107. 'Parameters: matrices Alpha,Beta
  1108. '
  1109. 'Returns: Alpha() = Alpha() - Beta()
  1110. '===================================================================
  1111. FUNCTION MatSubD% (Alpha() AS DOUBLE, Beta() AS DOUBLE)
  1112. ON LOCAL ERROR GOTO dsuberr: MatSubD% = 0
  1113. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  1114. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  1115. 'loop through and subtract elements
  1116. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  1117.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  1118.       Alpha(row%, col%) = Alpha(row%, col%) - Beta(row%, col%)
  1119.    NEXT col%
  1120. NEXT row%
  1121. dsubexit:
  1122. EXIT FUNCTION
  1123. dsuberr:
  1124.    MatSubD% = (ERR + 5) MOD 200 - 5
  1125.    RESUME dsubexit:
  1126. END FUNCTION
  1127.  
  1128. '=======================MatSubI%====================================
  1129. 'MatSubI% takes the difference of two integer matrices and places the
  1130. 'result in the first.
  1131. '
  1132. 'Parameters: matrices Alpha,Beta
  1133. '
  1134. 'Returns: Alpha() = Alpha() - Beta()
  1135. '===================================================================
  1136. FUNCTION MatSubI% (Alpha() AS INTEGER, Beta() AS INTEGER)
  1137. ON LOCAL ERROR GOTO isuberr: MatSubI% = 0
  1138. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  1139. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  1140. 'loop through and subtract elements
  1141. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  1142.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  1143.       Alpha(row%, col%) = Alpha(row%, col%) - Beta(row%, col%)
  1144.    NEXT col%
  1145. NEXT row%
  1146. isubexit:
  1147. EXIT FUNCTION
  1148. isuberr:
  1149.    MatSubI% = (ERR + 5) MOD 200 - 5
  1150.    RESUME isubexit:
  1151. END FUNCTION
  1152.  
  1153. '=======================MatSubL%====================================
  1154. 'MatSubL% takes the difference of two long integer matrices and places
  1155. 'the result in the first.
  1156. '
  1157. 'Parameters: matrices Alpha,Beta
  1158. '
  1159. 'Returns: Alpha() = Alpha() - Beta()
  1160. '===================================================================
  1161. FUNCTION MatSubL% (Alpha() AS LONG, Beta() AS LONG)
  1162. ON LOCAL ERROR GOTO lsuberr: MatSubL% = 0
  1163. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  1164. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  1165. 'loop through and subtract elements
  1166. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  1167.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  1168.       Alpha(row%, col%) = Alpha(row%, col%) - Beta(row%, col%)
  1169.    NEXT col%
  1170. NEXT row%
  1171. lsubexit:
  1172. EXIT FUNCTION
  1173. lsuberr:
  1174.    MatSubL% = (ERR + 5) MOD 200 - 5
  1175.    RESUME lsubexit:
  1176. END FUNCTION
  1177.  
  1178. '=======================MatSubS%====================================
  1179. 'MatSubS% takes the difference of two single precision matrices and
  1180. 'places the result in the first.
  1181. '
  1182. 'Parameters: matrices Alpha,Beta
  1183. '
  1184. 'Returns: Alpha() = Alpha() - Beta()
  1185. '===================================================================
  1186. FUNCTION MatSubS% (Alpha() AS SINGLE, Beta() AS SINGLE)
  1187. ON LOCAL ERROR GOTO ssuberr: MatSubS% = 0
  1188. 'check if Alpha, Beta have same dimensions if not, exit and send back error
  1189. IF (LBOUND(Alpha, 1) <> LBOUND(Beta, 1)) OR (UBOUND(Alpha, 1) <> UBOUND(Beta, 1)) OR (LBOUND(Alpha, 2) <> LBOUND(Beta, 2)) OR (UBOUND(Alpha, 2) <> UBOUND(Beta, 2)) THEN ERROR 196
  1190. 'loop through and subtract elements
  1191. FOR row% = LBOUND(Alpha, 1) TO UBOUND(Alpha, 1)
  1192.    FOR col% = LBOUND(Alpha, 2) TO UBOUND(Alpha, 2)
  1193.       Alpha(row%, col%) = Alpha(row%, col%) - Beta(row%, col%)
  1194.    NEXT col%
  1195. NEXT row%
  1196. ssubexit:
  1197. EXIT FUNCTION
  1198. ssuberr:
  1199.    MatSubS% = (ERR + 5) MOD 200 - 5
  1200.    RESUME ssubexit:
  1201. END FUNCTION
  1202.  
  1203.