home *** CD-ROM | disk | FTP | other *** search
/ Boston 2 / boston-2.iso / DOS / PROGRAM / C / NEWMAT / NEWMAT.TXT < prev    next >
Text File  |  1993-12-01  |  54KB  |  1,481 lines

  1. //$$ newmat.txt            Documentation file
  2.  
  3.  
  4.    Documentation for newmat03, an experimental matrix package in C++.
  5.    ==================================================================
  6.  
  7.  
  8. MATRIX PACKAGE                                           25 November, 1991
  9.  
  10. Copyright (C) 1991: R B Davies and DSIR
  11.  
  12. Permission is granted to use but not to sell.
  13.  
  14.  
  15. Contents
  16. ========
  17.  
  18. General description
  19. Is this the package you need?
  20. Changes
  21. Where you can get a copy of this package
  22. Compiler performance
  23. Detailed documentation
  24.    Customising
  25.    Constructors
  26.    Elements of matrices
  27.    Matrix copy
  28.    Unary operators
  29.    Binary operators
  30.    Combination of a matrix and scalar
  31.    Scalar functions of matrices
  32.    Submatrix operations
  33.    Change dimensions
  34.    Change type
  35.    Multiple matrix solve
  36.    Memory management
  37.    Output
  38.    Accessing matrices of unspecified type
  39.    Cholesky decomposition
  40.    Householder triangularisation
  41.    Singular Value Decomposition
  42.    Eigenvalues
  43.    Sorting
  44.    Fast Fourier Transform
  45.    Interface to Numerical Recipes in C
  46. List of files
  47. Notes on the design of the package
  48.    What this is package for
  49.    What size of matrices?
  50.    Allow matrix expressions?
  51.    Which matrix types?
  52.    What element types?
  53.    Naming convention
  54.    Row and Column index ranges
  55.    Structure of matrix objects
  56.    Data storage - one block or several
  57.    Data storage - by row or by column or other
  58.    Storage of symmetric matrices
  59.    Element access - method and checking
  60.    Use iterators?
  61.    Memory management - reference counting or status variable?
  62.    Evaluation of expressions - use two stage method?
  63.    How to overcome an explosion in number of operations
  64.    Using const
  65.    A calculus of matrix types
  66.    Error handling
  67.    Band and sparse matrices
  68. Problem report form
  69.  
  70.  
  71. ---------------------------------------------------------------------------
  72.  
  73.  
  74. General description
  75. ===================
  76.  
  77. The package is intented for scientists and engineers who need to
  78. manipulate a variety of types of matrices using standard matrix
  79. operations. Emphasis is on the kind of operations needed in statistical
  80. calculations such as least squares, linear equation solve and
  81. eigenvalues.
  82.  
  83. It supports matrix types
  84.  
  85.     Matrix                       (rectangular matrix)
  86.     nricMatrix                   (variant of rectangular matrix)
  87.     UpperTriangularMatrix
  88.     LowerTriangularMatrix
  89.     DiagonalMatrix
  90.     SymmetricMatrix
  91.     RowVector                    (derived from Matrix)
  92.     ColumnVector                 (derived from Matrix).
  93.  
  94. Only one element type (float or double) is supported.
  95.  
  96. The package includes the operations *, +, -, inverse, transpose,
  97. conversion between types, submatrix, determinant, Cholesky
  98. decomposition, Householder triangularisation, singular value
  99. decomposition, eigenvalues of a symmetric matrix, sorting, fast fourier
  100. transform, printing and an interface with "Numerical Recipes in C".
  101.  
  102. It is intended for matrices in the range 4 x 4 to about 90 x 90 (125 x
  103. 125 for triangular matrices). The upper limit is imposed by the maximum
  104. number of elements that can be contained in a single array (8192 doubles
  105. in some machines).
  106.  
  107. A two-stage approach to evaluating matrix expressions is used to improve
  108. efficiency and reduce use of temporary storage.
  109.  
  110. The package is designed for version 2 of C++. It works with Turbo C++,
  111. Borland C++, Glockenspiel C++ (2.00a) on a PC and AT&T C++ (2.0) and Gnu
  112. C++ on a Sun. It works with some problems with Zortech C++ (version 2).
  113.  
  114.  
  115. ---------------------------------------------------------------------------
  116.  
  117.  
  118. Is this the package you need?
  119. =============================
  120.  
  121. Do you
  122.  
  123. 1.   need matrix operators such as * and + defined as operators so you
  124.      can write things like
  125.  
  126.         X  = A * (B + C);
  127.  
  128. 2.   need a variety of types of matrices
  129.  
  130. 3.   need only one element type (float or double)
  131.  
  132. 4.   work with matrices in the range 4x4 to 90x90
  133.  
  134. 5.   tolerate a large and complex package
  135.  
  136.  
  137. Then maybe this is the right package for you. 
  138.  
  139. If you don't need (1) then there may be better options. Likewise if you
  140. don't need (2) there may be better options. If you require "not (5)"
  141. then this is not the package for you.
  142.  
  143.  
  144. If you need (2) and "not (3)" and have some spare money, then maybe you
  145. should look at M++ from Dyad or the Rogue Wave matrix package.
  146.  
  147.  
  148. ---------------------------------------------------------------------------
  149.  
  150.  
  151. Changes
  152. =======
  153.  
  154. Newmat03 - November 1991:
  155.  
  156. Col and Cols become Column and Columns. Added Sort, SVD, Jacobi,
  157. Eigenvalues, FFT, real conversion of 1x1 matrix, "Numerical Recipes in
  158. C" interface, output operations, various scalar functions. Improved
  159. return from functions. Reorganised setting options in "include.hxx".
  160.  
  161.  
  162. Newmat02 - July 1991:
  163.  
  164. Version with matrix row/column operations and numerous additional
  165. functions.
  166.  
  167.  
  168. Matrix - October 1990:
  169.  
  170. Early version of package.
  171.  
  172.  
  173. ---------------------------------------------------------------------------
  174.  
  175.  
  176. How to get a copy of this package
  177. =================================
  178.  
  179. I am putting copies on Compuserve (Borland library, zip format),
  180. SIMTEL20 (MsDos library, zip format), comp.sources.misc on Internet
  181. (shar format), and on the MsDos program library at Victoria University,
  182. Wellington.
  183.  
  184.  
  185. ---------------------------------------------------------------------------
  186.  
  187.  
  188. Compiler performance
  189. ====================
  190.  
  191. I have tested this package on a number of compilers. Here are the
  192. levels of success with this package. In most cases I have chosen code
  193. that works under all the compilers I have access to, but I have had to
  194. include some specific work-arounds for some compilers. For the MsDos
  195. versions, I am using a 386/387sx computer running MsDos 5, except that
  196. Turbo is on an old XT. The unix versions are on a Sun Sparc station.
  197.  
  198. A series of #defines at the beginning of "include.hxx" customises the
  199. package for the compiler you are using. Turbo, Borland and Zortech are
  200. recognised automatically, otherwise you have to set the appropriate
  201. #define statement.
  202.  
  203. The compilers are looking a bit old now. I do intend to test the package
  204. against newer versions as they become available.
  205.  
  206. Borland C++ 2.0: Recently, this has been my main development platform,
  207. so naturally almost everything works with this compiler. The library
  208. manager "tlib" fails but you can use "zorlib" from Zortech instead.
  209. Sometimes Borland crashes during a compiler or mis-compiles. You just
  210. have to reboot and continue the compile.
  211.  
  212. Turbo C++ (? version): Almost works OK. My rather elderly version does
  213. show a problem. Probably not worth tracking down - buy a newer version.
  214. Haven't tried the linker.
  215.  
  216. Zortech C++ 2.12: "const" doesn't work correctly with this compiler, so
  217. the package skips all of the statements Zortech can't handle. If you are
  218. using a later version of Zortech you could probably re-activate these
  219. statements. Zortech leaves rubbish on the heap. I don't know whether
  220. this is my programming error or a Zortech error. It works better when
  221. one doesn't optimise but there still are problems. The nric routines
  222. don't work. Zortech does not support IO manipulators.
  223.  
  224. Glockenspiel C++ (2.00a for MsDos loading into Microsoft C 5.1): I
  225. haven't tested the latest version of my package with Glockenspiel. I had
  226. to #define the matrix names to shorter names to avoid ambiguities and
  227. had quite a bit of difficulty stopping the compiles from running out of
  228. space and not exceeding Microsoft's block nesting limit. A couple of my
  229. test statements produced statements too complex for Microsoft, but
  230. basically the package worked. This was my original development platform
  231. and I still use .cxx and .hxx as my file name extensions. However,
  232. Glockenspiel is no longer competitive for MsDos and I am not updating my
  233. copy of the compiler.
  234.  
  235. Sun AT&T C++ 2.00: This works fine. Except aggregates are not supported.
  236.  
  237. Gnu G++ 1.37.1: This mostly works. You don't seem to be able to use
  238. expressions like Matrix(X*Y) in the middle of an expression and
  239. (Matrix)(X*Y) is unreliable. Gnu does not support IO manipulators. Gnu
  240. leaves rubbish on the heap. This is from output statements and not my
  241. package and may not be an error. The previous version of the package did
  242. not work under Gnu 1.37 or 1.39.
  243.  
  244.  
  245. ---------------------------------------------------------------------------
  246.  
  247.  
  248. Detailed Documentation
  249. ======================
  250.  
  251. Copyright (C) 1989,1990,1991: R B Davies and DSIR
  252.  
  253. Permission is granted to use but not to sell.
  254.  
  255.    --------------------------------------------------------------
  256.   | Please understand that this is a test version; there may     |
  257.   | still be bugs and errors. Use at your own risk. Neither I    |
  258.   | nor DSIR take any responsibility for any errors or omissions |
  259.   | in this package or for any misfortune that may befall you or |
  260.   | others as a result of its use.                               |
  261.    --------------------------------------------------------------
  262.  
  263. Please report bugs to me at
  264.  
  265.     robert@am.dsir.govt.nz
  266.  
  267. or
  268.  
  269.     Compuserve 72777,656
  270.  
  271. When reporting a bug please tell me which C++ compiler you are using (if
  272. known), and what version. Also give me details of your computer (if
  273. known). Tell me where you downloaded your version of my package from and
  274. its version number (eg newmat03 or newmat04). (There may be very minor
  275. differences between versions at different sites). Note any changes you
  276. have made to my code. If at all possible give me a piece of code
  277. illustrating the bug.
  278.  
  279. Please do report bugs to me.
  280.  
  281.  
  282. The matrix inverse routine and the sort routines are adapted from
  283. "Numerical Recipes in C" by Press, Flannery, Teukolsky, Vetterling,
  284. published by the Cambridge University Press.
  285.  
  286. Other code is adapted from routines in "Handbook for Automatic
  287. Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published
  288. by Springer Verlag. 
  289.  
  290.  
  291. Customising
  292. -----------
  293.  
  294. I use .hxx as the suffix of definition files and .cxx as the suffix of
  295. C++ source files. This does not cause any problems with the compilers I
  296. use except that Borland and Turbo need to be told to accept any suffix
  297. as meaning a C++ file rather than a C file.
  298.  
  299. Use the large model when you are using a PC. Do not "outline" inline
  300. functions.
  301.  
  302. Each file accessing the matrix package needs to have file newmat.hxx 
  303. #included  at the beginning. Files using matrix applications (Cholesky
  304. decomposition, Householder triangularisation) need newmatap.hxx instead
  305. (or as well). If you need the output functions you will also need
  306. newmatio.hxx. Glockenspiel also needs to have include.hxx #included before
  307. newmat.hxx.
  308.  
  309. The file  include.hxx  sets the options for the compiler. If you are using
  310. a compiler different from one I have worked with you may have to set up
  311. a new section in  include.hxx  appropriate for your compiler.
  312.  
  313. Borland, Turbo and Zortech are recognised automatically. If you using
  314. Glockenspiel on a PC, AT&T, or Gnu C++ activate the appropriate
  315. statement at the beginning of include.hxx.
  316.  
  317. Activate the appropriate statement to make the element type float or
  318. double.
  319.  
  320. The file (newmat9.cxx) containing the output routines can be used only
  321. with libraries that support the AT&T input/output routines including
  322. manipulators. It cannot be used with Zortech or Gnu.
  323.  
  324.  
  325. Constructors
  326. ------------
  327.  
  328. To construct an m x n matrix, A, (m and n are integers) use
  329.  
  330.     Matrix A(m,n);
  331.  
  332. The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and
  333. DiagonalMatrix types are symmetric. To construct an n x n matrix use,
  334. for example
  335.  
  336.     UpperTriangularMatrix U(n);
  337.  
  338. Likewise the RowVector and ColumnVector types take just one argument in
  339. their constructors:
  340.  
  341.     RowVector RV(n);
  342.  
  343. You can also construct vectors and matrices without specifying the
  344. dimension. For example
  345.  
  346.     Matrix A;
  347.  
  348. In this case the dimension must be set by an assignment statement or a
  349. re-dimension statement.
  350.  
  351. You can also use a constructor to set a matrix equal to another matrix
  352. or matrix expression.
  353.  
  354.     Matrix A = U;
  355.  
  356.     Matrix A = U * L;
  357.  
  358. Only conversions that don't lose information are supported - eg you
  359. cannot convert an upper triangular matrix into a diagonal matrix using =.
  360.  
  361.  
  362. Elements of matrices
  363. --------------------
  364.  
  365. Elements are accessed by expressions of the form A(i,j) where i and j
  366. run from 1 to the appropriate dimension. Access elements of vectors with
  367. just one argument. Diagonal matrices can accept one or two subscripts.
  368.  
  369. This is different from the earlier version of the package in which the
  370. subscripts ran from 0 to one less than the appropriate dimension. Use
  371. A.element(i,j) if you want this earlier convention.
  372.  
  373.  
  374. Matrix copy
  375. -----------
  376.  
  377. The operator = is used for copying matrices, converting matrices, or
  378. evaluating expressions. For example
  379.  
  380.     A = B;  A = L;  A = L * U;
  381.  
  382. Only conversions that don't lose information are supported. The
  383. dimensions of the matrix on the left hand side are adjusted to those of
  384. the matrix or expression on the right hand side. Elements on the right
  385. hand side which are not present on the left hand side are set to zero.
  386.  
  387. The operator << can be used in place of = where it is permissible for
  388. information to be lost.
  389.  
  390. For example
  391.  
  392.     SymmetricMatrix S; Matrix A;
  393.     ......
  394.     S << A.t() * A;
  395.  
  396. is acceptable whereas
  397.  
  398.     S = A.t() * A;                            // error
  399.  
  400. will cause a runtime error since the package does not (yet) recognise
  401. A.t()*A as symmetric.
  402.  
  403. Note that you can not use << with constructors. For example
  404.  
  405.     SymmetricMatrix S << A.t() * A;           // error
  406.  
  407. does not work.
  408.  
  409. A third copy routine is used in a similar role to =. Use
  410.  
  411.     A.Inject(D);
  412.  
  413. to copy the elements of D to the corresponding elements of A but leave
  414. the elements of A unchanged if there is no corresponding element of D
  415. (the = operator would set them to 0). This is useful, for example, for
  416. setting the diagonal elements of a matrix without disturbing the rest of
  417. the matrix. Unlike = and <<, Inject does not reset the dimensions of A, which
  418. must match those of D.
  419.  
  420. Both << and Inject can be used with submatrix expressions on the left
  421. hand side. See the section on submatrices.
  422.  
  423. To set the elements of a matrix to a scalar use operator =
  424.  
  425.     real r; Matrix A(m,n);
  426.     ......
  427.     Matrix A(m,n); A = r;
  428.  
  429. You can load the elements of a matrix from an array:
  430.  
  431.     Matrix A(3,2);
  432.     real a[] = { 11,12,21,22,31,33 };
  433.     A << a;
  434.  
  435. This construction cannot check that the numbers of elements match
  436. correctly. This version of << can be used with submatrices on the left
  437. hand side.
  438.  
  439.  
  440. Unary operators
  441. ---------------
  442.  
  443. The package supports unary operations
  444.  
  445.     change sign of elements            -A
  446.     transpose                          A.t()
  447.     inverse (of square matrix A)       A.i()
  448.  
  449.  
  450. Binary operations
  451. -----------------
  452.  
  453. The package supports binary operations
  454.  
  455.     matrix addition                    A+B
  456.     matrix subtraction                 A-B
  457.     matrix multiplication              A*B
  458.     equation solve (square matrix A)   A.i()*B
  459.  
  460. In the last case the inverse is not calculated.
  461.  
  462. Notes:
  463.  
  464. If you are doing repeated multiplication. For example A*B*C, use
  465. brackets to force the order to minimize the number of operations. If C
  466. is a column vector and A is not a vector, then it will usually reduce
  467. the number of operations to use A*(B*C) .
  468.  
  469. The package does not recognise B*A.i() as an equation solve. It is
  470. probably better to use (A.t().i()*B.t()).t() .
  471.  
  472.  
  473. Combination of a matrix and scalar
  474. ----------------------------------
  475.  
  476. The following expression multiplies the elements of a matrix A by a
  477. scalar f:  A * f; Likewise one can divide the elements of a matrix A by
  478. a scalar f:  A / f;
  479.  
  480. The expressions  A + f and A - f add or subtract a rectangular matrix of
  481. the same dimension as A with elements equal to f to or from the matrix
  482. A.
  483.  
  484. In each case the matrix must be the first term in the expression.
  485. Expressions such  f + A  or  f * A  are not recognised.
  486.  
  487.  
  488. Scalar functions of matrices
  489. ----------------------------
  490.             
  491.     int m = A.Nrows();                    // number of rows
  492.     int n = A.Ncols();                    // number of columns
  493.     real ssq = A.SumSquare();             // sum of squares of elements
  494.     real sav = A.SumAbsoluteValue();      // sum of absolute values
  495.     real mav = A.MaximumAbsoluteValue();  // maximum of absolute values
  496.     real norm = A.Norm1();                // maximum of sum of absolute
  497.                                              values of elements of a column
  498.     real norm = A.NormInfinity();         // maximum of sum of absolute
  499.                                              values of elements of a row
  500.     real t = A.Trace();                   // trace
  501.     LogandSign ld = A.LogDeterminant();   // log of determinant
  502.     BOOL z = A.IsZero();                  // test all elements zero
  503.     MatrixType mt = A.Type();             // type of matrix
  504.     real* s = Store();                    // pointer to array of elements
  505.     int l = Storage();                    // length of array of elements
  506.  
  507. A.LogDeterminant() returns a value of type LogandSign. If ld is of type 
  508. LogAndSign  use
  509.  
  510.     ld.Value()    to get the value of the determinant
  511.     ld.Sign()     to get the sign of the determinant (values 1, 0, -1)
  512.     ld.LogValue() to get the log of the absolute value.
  513.  
  514. A.IsZero() returns BOOL value TRUE if the matrix A has all elements
  515. equal to 0.0.
  516.  
  517. MatrixType mt = A.Type() returns the type of a matrix. Use (char*)mt to
  518. get a string  (UT, LT, Rect, Sym, Diag, RowV, ColV, Crout) showing the
  519. type.
  520.  
  521. SumSquare(A), SumAbsoluteValue(A), MaximumAbsoluteValue(A), Trace(A),
  522. LogDeterminant(A), Norm1(A), NormInfinity(A)  can be used in place of
  523. A.SumSquare(), A.SumAbsoluteValue(), A.MaximumAbsoluteValue(),
  524. A.Trace(), A.LogDeterminant(), A.Norm1(), A,NormInfinity().
  525.  
  526.  
  527. Submatrix operations
  528. --------------------
  529.  
  530. A.SubMatrix(fr,lr,fc,lc)
  531.  
  532. This selects a submatrix from A. the arguments  fr,lr,fc,lc  are the
  533. first row, last row, first column, last column of the submatrix with the
  534. numbering beginning at 1. This may be used in any matrix expression or
  535. on the left hand side of << or Inject. Inject does not check no
  536. information loss in this case. You can also use the construction
  537.  
  538.     real c; .... A.SubMatrix(fr,lr,fc,lc) << c;
  539.  
  540. to set a submatrix equal to a constant.
  541.  
  542. The follwing are variants of SubMatrix:
  543.  
  544.     A.SymSubMatrix(f,l)             //   This assumes fr=fc and lr=lc.
  545.     A.Rows(f,l)                     //   select rows
  546.     A.Row(f)                        //   select single row
  547.     A.Columns(f,l)                  //   select columns
  548.     A.Column(f)                     //   select single column
  549.  
  550. In each case f and l mean the first and last row or column to be
  551. selected (starting at 1).
  552.  
  553. If SubMatrix or its variant occurs on the right hand side of an = or <<
  554. or within an expression its type is as follows
  555.  
  556.     A.Submatrix(fr,lr,fc,lc):           If A is RowVector or
  557.                                         ColumnVector then same type
  558.                                         otherwise type Matrix
  559.     A.SymSubMatrix(f,l):                Same type as A
  560.     A.Rows(f,l):                        Type Matrix
  561.     A.Row(f):                           Type RowVector
  562.     A.Columns(f,l):                     Type Matrix
  563.     A.Column(f):                        Type ColumnVector
  564.  
  565.  
  566. If SubMatrix or its variant appears on the left hand side of  << , think
  567. of its type being Matrix. Thus L.Row(1) where L is LowerTriangularMatrix
  568. expects  L.Ncols()  elements even though it will use only one of them.
  569.  
  570.  
  571. Change dimensions
  572. -----------------
  573.  
  574. The following operations change the dimensions of a matrix. The values
  575. of the elements are lost.
  576.  
  577.     A.ReDimension(nrows,ncols);     // for type Matrix or nricMatrix
  578.     A.ReDimension(n);               // for all other types
  579.  
  580.  
  581. Change type
  582. -----------
  583.  
  584. The following functions interpret the elements of a matrix
  585. (stored row by row) to be a vector or matrix of a different type. Actual
  586. copying is usually avoided where these occur as part of a more
  587. complicated expression.
  588.  
  589.     A.CopyToRow()
  590.     A.CopyToColumn()
  591.     A.CopyToDiagonal()
  592.     A.CopyToMatrix(nrows,ncols)
  593.     A.c()
  594.     real(A)
  595.  
  596. The form .c() is used in matrix expressions when A is of a const
  597. type. The expression real(A) is used to convert a 1 x 1 matrix to a
  598. scalar.
  599.  
  600.  
  601. Multiple matrix solve
  602. ---------------------
  603.  
  604. If A is a square or symmetric matrix use
  605.  
  606.     CroutMatrix X = A;                // carries out LU decomposition
  607.     Matrix AP = X.i()*P; Matrix AQ = X.i()*Q;
  608.     LogAndSign ld = X.LogDeterminant();
  609.  
  610. rather than
  611.  
  612.     Matrix AP = A.i()*P; Matrix AQ = A.i()*Q;
  613.     LogAndSign ld = A.LogDeterminant();
  614.  
  615. since each operation will repeat the LU decompostion.
  616.  
  617.  
  618. Memory management
  619. -----------------
  620.  
  621. The package does not support delayed copy. Several strategies are
  622. required to prevent unnecessary matrix copies.
  623.  
  624. Where a matrix is called as a function argument use a constant
  625. reference. For example
  626.  
  627.     YourFunction(const Matrix& A)
  628.  
  629. rather than
  630.  
  631.     YourFunction(Matrix A)
  632.  
  633. Constant matrices cannot be used in matrix expressions so if you wish to
  634. use A in an expression within this function use A.c() rather than A.
  635.  
  636. Skip the rest of this section on your first reading.
  637.  
  638. A second place where it is desirable to avoid unnecessary copies is when
  639. a function is returning a matrix. Matrices can be returned from a
  640. function with the return command as you would expect. However these may
  641. incur one and possibly two copyings of the matrix. To avoid this use the
  642. following instructions.
  643.  
  644. Make your function of type  ReturnMatrix . Then precede the return
  645. statement with a Release statement (or a ReleaseAndDelete statement if
  646. the matrix was created with new). For example
  647.  
  648.  
  649.     ReturnMatrix MakeAMatrix()
  650.     {
  651.        Matrix A;
  652.        ......
  653.        A.Release(); return A;
  654.     }
  655.  
  656. or
  657.  
  658.     ReturnMatrix MakeAMatrix()
  659.     {
  660.        Matrix* m = new Matrix;
  661.        ......
  662.        m->ReleaseAndDelete(); return *m;
  663.     }
  664.  
  665. Note that .c() cannot be applied to a matrix following application of
  666. .Release() or ->ReleaseAndDelete() .
  667.  
  668.  --------------------------------------------------------------------- 
  669. | Do not forget to make the function of type ReturnMatrix; otherwise  |
  670. | incomprehensible run-time errors will occur with some compilers.    |
  671.  --------------------------------------------------------------------- 
  672.  
  673. You can also use .Release() or ->ReleaseAndDelete() to allow a matrix
  674. expression to recycle space. Suppose you call
  675.  
  676.     A.Release();
  677.  
  678. just before A is used just once in an expression. Then the memory used
  679. by A is either returned to the system or reused in the expression. In
  680. either case, A's memory is destroyed. This procedure can be used to
  681. imporve efficiency and reduce the use of memory.
  682.  
  683. Use ->ReleaseAndDelete for matrices created by new if you want to
  684. completely delete the matrix after it is accessed.
  685.  
  686.  
  687. Output
  688. ------
  689.  
  690. To print a matrix use an expression like
  691.  
  692.     Matrix A;
  693.     ......
  694.     cout << setw(10) << setprecision(5) << A;
  695.  
  696. This will work only with systems that support the AT&T input/output
  697. routines including manipulators.
  698.  
  699.  
  700. Accessing matrices of unspecified type
  701. --------------------------------------
  702.  
  703. Skip this section on your first reading.
  704.  
  705. Suppose you wish to write a function which accesses a matrix of unknown
  706. type including expressions (eg A*B). Then use a layout similar to the
  707. following:
  708.  
  709.    void YourFunction(BaseMatrix& X)
  710.    {
  711.       GeneralMatrix* gm = X.Evaluate();   // evaluate an expression
  712.                                           // if necessary
  713.       ........                            // operations on *gm
  714.       gm->tDelete();                      // delete *gm if a temporary
  715.    }
  716.  
  717. See, as an example, the definitions of operator<< in newmat9.cxx.
  718.  
  719.  
  720. Cholesky decomposition
  721. ----------------------
  722.  
  723. Suppose S is symmetric and positive definite. Then there exists a unique
  724. lower triangular matrix L such that L * L.t() = S. To calculate this use
  725.  
  726.     SymmetricMatrix S;
  727.     ......
  728.     LowerTriangularMatrix L = Cholesky(S);
  729.  
  730.  
  731. Householder triangularisation
  732. -----------------------------
  733.  
  734. Start with matrix
  735.  
  736.        / X    0 \      s
  737.        \ Y    0 /      t
  738.  
  739.          n    s
  740.  
  741. The Householder triangularisation post multiplies by an orthogonal
  742. matrix Q such that the matrix becomes
  743.  
  744.        / 0    L \      s
  745.        \ Z    M /      t
  746.  
  747.          n    s
  748.  
  749. where L is lower triangular. Note that X is the transpose of the matrix
  750. sometimes considered in this context.
  751.  
  752. This is good for solving least squares problems: choose b (matrix or row
  753. vector) to minimize the sum of the squares of the elements of
  754.  
  755.          Y - b*X
  756.  
  757. Then choose b = M * L.i();
  758.  
  759. Two routines are provided:
  760.  
  761.     HHDecompose(X, L);
  762.  
  763. replaces X by orthogonal columns and forms L.
  764.  
  765.     HHDecompose(X, Y, M);
  766.  
  767. uses X from the first routine, replaces Y by Z and forms M.
  768.  
  769.  
  770. Singular Value Decomposition
  771. ----------------------------
  772.  
  773. The singular value decomposition of an m x n matrix A ( where m >= n) is
  774. a decomposition
  775.  
  776.     A  = U * D * V.t()
  777.  
  778. where U is m x n with  U.t() * U  equalling the identity, D is an n x n
  779. diagonal matrix and V is an n x n orthogonal matrix.
  780.  
  781. Singular value decompositions are useful for understanding the structure
  782. of ill-conditioned matrices, solving least squares problems, and for
  783. finding the eigenvalues of A.t() * A.
  784.  
  785. To calculate the singular value decomposition of A (with m >= n) use one
  786. of
  787.  
  788.     SVD(A, D, U, V);                  // U (= A is OK)
  789.     SVD(A, D);
  790.     SVD(A, D, U);                     // U (= A is OK)
  791.     SVD(A, D, U, FALSE);              // U (can = A) for workspace only
  792.     SVD(A, D, U, V, FALSE);           // U (can = A) for workspace only
  793.  
  794. The values of A are not changed unless A is also inserted as the third
  795. argument.
  796.  
  797.  
  798. Eigenvalues
  799. -----------
  800.  
  801. An eigenvalue decomposition of a symmetric matrix A is a decomposition
  802.  
  803.     A  = V * D * V.t()
  804.  
  805. where V is an orthogonal matrix and D is a diagonal matrix.
  806.  
  807. Eigenvalue analyses are used in a wide variety of engineering,
  808. statistical and other mathematical analyses.
  809.  
  810. The package includes two algorithms: Jacobi and Householder. The first
  811. is extremely reliable but much slower than the second.
  812.  
  813. The code is adapted from routines in "Handbook for Automatic
  814. Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published
  815. by Springer Verlag. 
  816.  
  817.  
  818.     Jacobi(A,D,S,V);                  // A, S symmetric is for workspace,
  819.                                       //    S = A is OK
  820.     Jacobi(A,D);                      // A symmetric
  821.     Jacobi(A,D,S);                    // A, S symmetric is for workspace,
  822.                                       //    S = A is OK
  823.     Jacobi(A,D,V);                    // A symmetric
  824.  
  825.     EigenValues(A,D);                 // A symmetric
  826.     EigenValues(A,D,S);               // A, S symmetric is for back
  827.                                       //    transforming, S = A is OK
  828.     EigenValues(A,D,V);               // A symmetric
  829.  
  830.  
  831. Sorting
  832. -------
  833.  
  834. To sort the values in a matrix or vector, A, (in general this operation
  835. makes sense only for vectors and diagonal matrices) use
  836.  
  837.     SortAscending(A);
  838.  
  839. or
  840.  
  841.     SortDescending(A);
  842.  
  843.  
  844. Fast Fourier Transform
  845. ----------------------
  846.  
  847. FFT(CV1, CV2, CV3, CV4);       // CV3=CV1 and CV4=CV2 is OK
  848.  
  849. where CV1, CV2, CV3, CV4 are column vectors. CV1 and CV2 are the real
  850. and imaginary input vectors; CV3 and CV4 are the real and imaginary
  851. output vectors. The lengths of CV1 and CV2 must be equal and should be
  852. the product of numbers less than about 10 for fast execution.
  853.  
  854.  
  855. Interface to Numerical Recipes in C
  856. -----------------------------------
  857.  
  858. This package can be used with the vectors and matrices defined in
  859. "Numerical Recipes in C". You need to edit the routines in Numerical
  860. Recipes so that the elements are of the same type as used in this
  861. package. Eg replace float by double, vector by dvector and matrix by
  862. dmatrix, etc. You will also need to edit the function definitions to use
  863. the version acceptable to your compiler. Then enclose the code from
  864. Numerical Recipes in  extern "C" { ... }. You will also need to include
  865. the matrix and vector utility routines.
  866.  
  867. Then any vector in Numerical Recipes with subscripts starting from 1 in
  868. a function call can be replaced by a RowVector, ColumnVector or
  869. DiagonalMatrix in the present package. Similarly any matrix with
  870. subscripts starting from 1 can be replaced by an  nricMatrix  in the
  871. present package. The class nricMatrix is derived from Matrix and can be
  872. used in place of Matrix.
  873.  
  874. Numerical Recipes cannot change the dimensions of a matrix or vector. So
  875. matrices or vectors must be correctly dimensioned before a Numerical
  876. Recipes routine is called.
  877.  
  878. For example
  879.  
  880.    SymmetricMatrix B(44);
  881.    .....                             // load values into B
  882.    nricMatrix BX = B;                // copy values to an nricMatrix
  883.    DiagonalMatrix D(44);             // Matrices for output
  884.    nricMatrix V(44,44);              //    correctly dimensioned
  885.    int nrot;
  886.    jacobi(BX,44,D,V,&nrot);          // jacobi from NRIC
  887.    cout << D;                        // print eigenvalues
  888.  
  889. ---------------------------------------------------------------------------
  890.  
  891.  
  892. List of files
  893. =============
  894.  
  895. NEWMAT   TXT    documentation file
  896. NEWMAT   LIS    list of files
  897.  
  898. BOOLEAN  HXX    boolean class definition
  899. CONTROLW HXX    control word definition file
  900. INCLUDE  HXX    details of include files and options
  901. NEWMAT   HXX    main matrix clss definition file
  902. NEWMATAP HXX    applications definition file
  903. NEWMATIO HXX    input/output definition file
  904. NEWMATRC HXX    row/column functions definition files
  905. NEWMATRM HXX    rectangular matrix access definition files
  906. PRECISIO HXX    numerical precision constants
  907.  
  908. CHOLESKY CXX    Cholesky decomposition
  909. EVALUE   CXX    eigenvalues and eigenvector calculation
  910. FFT      CXX    fast Fourier transform
  911. HHOLDER  CXX    Householder triangularisation
  912. JACOBI   CXX    eigenvalues by the Jacobi method
  913. NEWMAT1  CXX    type manipulation routines
  914. NEWMAT2  CXX    row and column manipulation functions
  915. NEWMAT3  CXX    row and column access functions
  916. NEWMAT4  CXX    constructors, redimension, utilities
  917. NEWMAT5  CXX    transpose, evaluate, matrix functions
  918. NEWMAT6  CXX    operators, element access
  919. NEWMAT7  CXX    invert, solve, binary operations
  920. NEWMAT8  CXX    LU decomposition, scalar functions
  921. NEWMAT9  CXX    output routines
  922. NEWMATRM CXX    rectangular matrix access functions
  923. SORT     CXX    sorting functions
  924. SUBMAT   CXX    submatrix functions
  925. SVD      CXX    singular value decomposition
  926.  
  927. EXAMPLE  CXX    example of use of package
  928. EXAMPLE  TXT    output from example
  929. EXAMPLE  DEP    dependency file for example
  930.  
  931. ---------------------------------------------------------------------------
  932.  
  933.  
  934. Notes on the design of the package
  935. ==================================
  936.  
  937. Copyright (C) 1991: R B Davies and DSIR
  938.  
  939. Please treat this as an academic publication. You can use the ideas but
  940. please acknowledge in any publication.
  941.  
  942.  
  943. In this section, I describe some of the ideas behind this package, some
  944. of the decisions that I needed to make and give some details about the
  945. way it works. You don't need to read this section in order to use the
  946. package.
  947.  
  948. I don't think it is obvious what is the best way of going about
  949. structuring a matrix package. And I don't think you can figure this
  950. out with "thought" experiments. Different people have to try out
  951. different approaches. And someone else may have to figure out which is
  952. best. Or, more likely, the ultimate packages will lift some ideas from
  953. each of a variety of trial packages. So, I don't claim my package is an
  954. "ultimate" package, but simply a trial of a number of ideas.
  955.  
  956. But I do hope it will meet the immediate requirements of some people who
  957. need to carry out matrix calculations.
  958.  
  959.  
  960. What this is package for
  961. ------------------------
  962.  
  963. The package is used for the manipulation of matrices, including the
  964. standard operations such as multiplication as understood by numerical
  965. analysts, engineers and mathematicians. A matrix is a two dimensional
  966. array of numbers. However, very special operations such as matrix
  967. multiplication are defined specifically for matrices. This means that a
  968. "matrix" package tends to be different from a general "array" package.
  969.  
  970. I see a matrix package as providing the following
  971.  
  972. 1.   Evaluation of matrix expressions in a form familiar to
  973.      scientists and engineers. For example  A = B * (C + D.t());
  974. 2.   Access to the elements of a matrix;
  975. 3.   Access to submatrices;
  976. 4.   Common elementary matrix functions such as determinant and trace;
  977. 5.   A platform for developing advanced matrix functions such as eigen-
  978.      value analysis;
  979. 6.   Good efficiency and storage management;
  980. 7.   Graceful exit from errors (I don't provide this yet).
  981.  
  982. It may also provide
  983.  
  984. 8.   A variety of types of elements (eg real and complex);
  985. 9.   A variety of types of matrices (eg rectangular and symmetric).
  986.  
  987. In contrast an array package should provide
  988.  
  989. 1'.  Arrays can be copied with a single instruction; may have some
  990.      arithmetic operations, say + - and scalar + - * /, it may provide
  991.      matrix multiplication as a function;
  992. 2'.  High speed access to elements directly and with iterators;
  993. 3'.  Access to subarrays and rows (and columns?);
  994. 6'.  Good efficiency and storage management;
  995. 7'.  Graceful exit from errors;
  996. 8'.  A variety of types of elements (eg real and complex);
  997. 9'.  One, two and three dimensional arrays, at least, with starting
  998.      points of the indices set by user; may have arrays with ragged
  999.      rows.
  1000.  
  1001. It may be possible to amalgamate these two sets of requirements to some
  1002. extent. However my package is definitely oriented towards the first set.
  1003.  
  1004. Even within the bounds set by the first set of requirements there is a
  1005. substantial opportunity for variation between what different matrix
  1006. packages might provide.
  1007.  
  1008. It is not possible to build a matrix package that will meet everyone's
  1009. requirements. In many cases if you put in one facility, you impose
  1010. overheads on everyone using the package. This both is storage required
  1011. for the program and in efficiency. Likewise a package that is optimised
  1012. towards handling large matrices is likely to become less efficient for
  1013. very small matrices where the administration time for the matrix may
  1014. become significant compared with the time to carry out the operations.
  1015.  
  1016. It is better to provide a variety of packages (hopefully compatible) so
  1017. that most users can find one that meets their requirements. This package
  1018. is intended to be one of these packages; but not all of them.
  1019.  
  1020. Since my background is in statistical methods, this package is oriented
  1021. towards the kinds things you need for statistical analyses.
  1022.  
  1023.  
  1024. What size of matrices?
  1025. ----------------------
  1026.  
  1027. A matrix package may target small matrices (say 3 x 3), or medium sized
  1028. matrices, or very large matrices. A package targeting very small
  1029. matrices will seek to minimise administration. A package for medium
  1030. sized or very large matrices can spend more time on administration in
  1031. order to conserve space or optimise the evaluation of expressions. A
  1032. package for very large matrices will need to pay special attention to
  1033. storage and numerical properties.
  1034.  
  1035. This package is designed for medium sized matrices. This means it is
  1036. worth introducing some optimisations, but I don't have to worry about
  1037. the 64K limit that some compilers impose.
  1038.  
  1039.  
  1040. Allow matrix expressions?
  1041. -------------------------
  1042.  
  1043. I want to be able to write matrix expressions the way I would on paper.
  1044. So if I want to multiply two matrices and then add the transpose of a
  1045. third one I can write something like
  1046.  
  1047.    X = A * B + C.t();
  1048.  
  1049. I want this expression to be evaluated with close to the same efficiency
  1050. as a hand-coded version. This is not so much of a problem with
  1051. expressions including a multiply since the multiply will dominate the
  1052. time. However, it is not so easy to achieve with expressions with just +
  1053. and - .
  1054.  
  1055. A second requirement is that temporary matrices generated during the
  1056. evaluation of an expression are destroyed as quickly as possible.
  1057.  
  1058. A desirable feature is that a certain amount of "intelligence" be
  1059. displayed in the evaluation of an expression. For example, in the
  1060. expression
  1061.  
  1062.    X = A.i() * B;
  1063.  
  1064. where i() denotes inverse, it would be desirable if the inverse wasn't
  1065. explicitly calculated.
  1066.  
  1067.  
  1068. Which matrix types?
  1069. -------------------
  1070.  
  1071. As well as the usual rectangular matrices, matrices occuring repeatedly
  1072. in numerical calculations are upper and lower triangular matrices,
  1073. symmetric matrices and diagonal matrices. This is particularly the case
  1074. in calculations involving least squares and eigenvalue calculations. So
  1075. as a first stage these were the types I decided to include.
  1076.  
  1077. It is also necessary to have types row vector and column vector. In a
  1078. "matrix" package, in contrast to an "array" package, it is necessary to
  1079. have both these types since they behave differently in matrix
  1080. expressions. The vector types can be derived for the rectangular matrix
  1081. type, so having them does not greatly increase the complexity of the
  1082. package.
  1083.  
  1084. The problem with having several matrix types is the number of versions
  1085. of the binary operators one needs. If one has 5 distinct matrix types
  1086. then a simple package will need 25 versions of each of the binary
  1087. operators. In fact, we can evade this problem, but at the cost of some
  1088. complexity.
  1089.  
  1090.  
  1091. What element types?
  1092. -------------------
  1093.  
  1094. Ideally we would allow element types double, float, complex and int, at
  1095. least. It would be reasonably easy, using templates or equivalent, to
  1096. provide a package which could handle a variety of element types.
  1097. However, as soon as one starts implementing the binary operators between
  1098. matrices with different element types, again one gets an explosion in
  1099. the number of operations one needs to consider. Hence I decided to
  1100. implement only one element type. But the user can decide whether this is
  1101. float or double. The package assumes elements are of type real. The user
  1102. typedefs to float or double.
  1103.  
  1104. In retrospect, it would not be too hard to include matrices with complex
  1105. matrix type.
  1106.  
  1107. It might also be worth including symmetric and triangular matrices with
  1108. extra precision elements (double or long double) to be used for storage
  1109. only and with a minimum of operations defined. These would be used for
  1110. accumulating the results of sums of squares and product matrices or
  1111. multistage Householder triangularisations.
  1112.  
  1113.  
  1114. Naming convention
  1115. -----------------
  1116.  
  1117. How are classes and public member functions to be named? As a general
  1118. rule I have spelt identifiers out in full with individual words being
  1119. capitalised. For example "UpperTriangularMatrix". If you don't like this
  1120. you can #define or typedef shorter names. This convention means you can
  1121. select an abbreviation scheme that makes sense to you.
  1122.  
  1123. The convention causes problems for Glockenspiel C++ on a PC feeding into
  1124. Microsoft C. The names Glockenspiel generates exceed the the 32
  1125. characters recognised by Microsoft C and ambiguities result. So it is
  1126. necessary to #define shorter names.
  1127.  
  1128. Exceptions to the general rule are the functions for transpose and
  1129. inverse. To make matrix expressions more like the corresponding
  1130. mathematical formulae, I have used the single letter abbreviations, t()
  1131. and i() .
  1132.  
  1133.  
  1134. Row and Column index ranges
  1135. ---------------------------
  1136.  
  1137. In mathematical work matrix subscripts usually start at one. In C, array
  1138. subscripts start at zero. In Fortran, they start at one. Possibilities
  1139. for this package were to make them start at 0 or 1 or be arbitrary.
  1140. Alternatively one could specify an "index set" for indexing the rows and
  1141. columns of a matrix. One would be able to add or multiply matrices only
  1142. if the appropriate row and column index sets were identical.
  1143.  
  1144. In fact, I adopted the simpler convention of making the rows and columns
  1145. of a matrix be indexed by an integer starting at one, following the
  1146. traditional convention. In an earlier version of the package I had them
  1147. starting at zero, but even I was getting mixed up when trying to use
  1148. this earlier package. So I reverted to the more usual notation.
  1149.  
  1150.  
  1151. Structure of matrix objects
  1152. ---------------------------
  1153.  
  1154. Each matrix object contains the basic information such as the number of
  1155. rows and columns and a status variable plus a pointer to the data
  1156. array which is on the heap.
  1157.  
  1158.  
  1159. Data storage - one block or several
  1160. -----------------------------------
  1161.  
  1162. In this package the elements of the matrix are stored as a single array.
  1163. Alternatives would have been to store each row as a separate array or a
  1164. set of adjacent rows as a separate array. The present solution
  1165. simplifies the program but limits the size of matrices in systems that
  1166. have a 64k byte (or other) limit on the size of arrays. The large arrays
  1167. may also cause problems for memory management in smaller machines.
  1168.  
  1169.  
  1170. Data storage - by row or by column or other
  1171. -------------------------------------------
  1172.  
  1173. In Fortran two dimensional arrays are stored by column. In most other
  1174. systems they are stored by row. I have followed this later convention.
  1175. This makes it easier to interface with other packages written in C but
  1176. harder to interface with those written in Fortran.
  1177.  
  1178. An alternative would be to store the elements by mid-sized rectangular
  1179. blocks. This might impose less strain on memory management when one
  1180. needs to access both rows and columns.
  1181.  
  1182.  
  1183. Storage of symmetric matrices
  1184. -----------------------------
  1185.  
  1186. Symmetric matrices are stored as lower triangular matrices. The decision
  1187. was pretty arbitrary, but it does slightly simplify the Cholesky
  1188. decomposition program.
  1189.  
  1190.  
  1191. Element access - method and checking
  1192. ------------------------------------
  1193.  
  1194. We want to be able to use the notation A(i,j) to specify the (i,j)-th
  1195. element of a matrix. This is the way mathematicians expect to address
  1196. the elements of matrices. I didn't even consider using the totally alien
  1197. notation A[i][j]. There are two ways of working out the address of
  1198. A(i,j). One is using a "dope" vector which contains the first address of
  1199. each row. This is how C works when you use A[i][j]. Alternatively you
  1200. can calculate the address using the formula appropriate for the
  1201. structure of A. I use this second approach. It is probably slower, but
  1202. saves worrying about an extra bit of storage. The other question is
  1203. whether to check for i and j being in range. I do carry out this check
  1204. following years of experience with both systems that do and systems that
  1205. don't do this check.
  1206.  
  1207. I would hope that the routines I supply with this package will reduce
  1208. your need to access elements of matrices so speed of access is not a
  1209. high priority.
  1210.  
  1211.  
  1212. Use iterators?
  1213. --------------
  1214.  
  1215. Iterators are an alternative way of providing fast access to the
  1216. elements of an array or matrix when they are to be accessed
  1217. sequentially. They need to be customised for each type of matrix. I have
  1218. not implemented iterators in this package, although some iterator like
  1219. functions are used for some row and column functions.
  1220.  
  1221.  
  1222. Memory management - reference counting or status variable?
  1223. ----------------------------------------------------------
  1224.  
  1225. Consider the instruction
  1226.  
  1227.    X = A + B + C;
  1228.  
  1229. To evaluate this a simple program will add A to B putting the total in a
  1230. temporary T1. Then it will add T1 to C creating another temporary T2
  1231. which will be copied into X. T1 and T2 will sit around till the end of
  1232. the block. It would be faster if the program recognised that T1 was
  1233. temporary and stored the sum of T1 and C back into T1 instead of
  1234. creating T2 and then avoided the final copy by just assigning the
  1235. contents of T1 to X rather than copying. In this case there will be no
  1236. temporaries requiring deletion. (More precisely there will be a header
  1237. to be deleted but no contents).
  1238.  
  1239. For an instruction like
  1240.  
  1241.    X = (A * B) + (C * D);
  1242.  
  1243. we can't avoid one temporary being left over, so we would like this
  1244. temporary deleted as quickly as possible.
  1245.  
  1246. I provide the functionality for doing this by attaching a status
  1247. variable to each matrix. This indicates if the matrix is temporary so
  1248. that its memory is available for recycling or deleting. Any matrix
  1249. operation checks the status variables of the matrices it is working with
  1250. and recycles or deletes any temporary memory.
  1251.  
  1252. An alternative or additional approach would be to use delayed copying.
  1253. If a program requests a matrix to be copied, the copy is delayed until
  1254. an instruction is executed which modifies the memory of either the
  1255. original matrix or the copy. This saves the unnecessary copying in the
  1256. previous examples. However, it does not provide the additional
  1257. functionality of my approach.
  1258.  
  1259. It would be possible to have both delayed copy and tagging temporaries,
  1260. but this seemed an unnecessary complexity. In particular delayed copy
  1261. mechanisms seem to require two calls to the heap manager, using extra
  1262. time and making building a memory compacting mechanism more difficult.
  1263.  
  1264.  
  1265. Evaluation of expressions - use two stage method?
  1266. -------------------------------------------------
  1267.  
  1268. Consider the instruction
  1269.  
  1270.    X = B - X;
  1271.  
  1272. The simple program will subtract X from B, store the result in a
  1273. temporary T1 and copy T1 into X. It would be faster if the program
  1274. recognised that the result could be stored directly into X. This would
  1275. happen automatically if the program could look at the instruction first
  1276. and mark X as temporary.
  1277.  
  1278. C programmers would expect to avoid the same problem with
  1279.  
  1280.    X = X - B;
  1281.  
  1282. by using an operator -= (which I haven't provided, yet)
  1283.  
  1284.    X -= B;
  1285.  
  1286. However this is an unnatural notation for non C users and it is much
  1287. nicer to write
  1288.  
  1289.    X = X - B;
  1290.  
  1291. and know that the program will carry out the simplification.
  1292.  
  1293. Another example where this intelligent analysis of an instruction is
  1294. helpful is in
  1295.  
  1296.    X = A.i() * B;
  1297.  
  1298. where i() denotes inverse. Numerical analysts know it is inefficient to
  1299. evaluate this expression by carrying out the inverse operation and then
  1300. the multiply. Yet it is a convenient way of writing the instruction. It
  1301. would be helpful if the program recognised this expression and carried
  1302. out the more appropriate approach.
  1303.  
  1304. To carry out this "intelligent" analysis of an instruction  matrix
  1305. expressions are evaluated in two stages. In the the first stage a tree
  1306. representation of the expression is formed.
  1307.  
  1308. For example (A+B)*C is represented by a tree
  1309.  
  1310.                     *
  1311.                    / \
  1312.                   +   C
  1313.                  / \
  1314.                 A   B
  1315.  
  1316. Rather than adding A and B the + operator yields an object of a class
  1317. "AddedMatrix" which is just a pair of pointers to A and B. Then the *
  1318. operator yields a "MultipliedMatrix" which is a pair of pointers to the
  1319. "AddedMatrix" and C. The tree is examined for any simplifications and
  1320. then evaluated recursively.
  1321.  
  1322. Further possibilities not yet included are to recognise A.t()*A and
  1323. A.t()+A as symmetric or to improve the efficiency of evaluation of
  1324. expressions like A+B+C, A*B*C, A*B.t()  [t() denotes transpose].
  1325.  
  1326. One of the disadvantages of the two-stage approach is that the types of
  1327. matrix expressions are determined at run-time. So the compiler will not
  1328. detect errors of the type
  1329.  
  1330.    Matrix M; DiagonalMatrix D; ....; D = M;
  1331.  
  1332. We don't allow conversions using = when information would be lost. Such
  1333. errors will be detected when the statement is executed.
  1334.  
  1335.  
  1336. How to overcome an explosion in number of operations
  1337. ----------------------------------------------------
  1338.  
  1339. The package attempts to solve the problem of the large number of
  1340. versions of the binary operations required when one has a variety of
  1341. types. With n types of matrices the binary operations will each require
  1342. n-squared separate algorithms. Some reduction in the number may be
  1343. possible by carrying out conversions. However the situation rapidly
  1344. becomes impossible with more than 4 or 5 types.
  1345.  
  1346. Doug Lea told me that it was possible to avoid this problem. I don't
  1347. know what his solution is. Here's mine.
  1348.  
  1349. Each matrix type includes routines for extracting individual rows or
  1350. columns. I assume a row or column consists of a sequence of zeros, a
  1351. sequence of stored values and then another sequence of zeros. Only a
  1352. single algorithm is then required for each binary operation. The rows
  1353. can be located very quickly since most of the matrices are stored row by
  1354. row. Columns must be copied and so the access is somewhat slower. As far
  1355. as possible my algorithms access the matrices by row.
  1356.  
  1357. An alternative approach of using iterators will be slower since the
  1358. iterators will involve a virtual function access at each step.
  1359.  
  1360. In fact, I provide several algorithms for operations like + . If one is
  1361. adding two matrices of the same type then there is no need to access the
  1362. individual rows or columns and a faster general algorithm is
  1363. appropriate.
  1364.  
  1365. Generally the method works well. However symmetric matrices are not
  1366. always handled very efficiently (yet) since complete rows are not stored
  1367. explicitly.
  1368.  
  1369. The original version of the package did not use this access by row or
  1370. column method and provided the multitude of algorithms for the
  1371. combination of different matrix types. The code file length turned out
  1372. to be just a little longer than the present one when providing the same
  1373. facilities with 5 distinct types of matrices. It would have been very
  1374. difficult to increase the number of matrix types in the original
  1375. version. Apparently 4 to 5 types is about the break even point for
  1376. switching to the approach adopted in the present package.
  1377.  
  1378.  
  1379. Using const
  1380. -----------
  1381.  
  1382. The memory management scheme introduces a problem when a matrix is
  1383. declared const. Because an operator may want to recycle the memory of
  1384. its operands these operands cannot be declared const. It isn't
  1385. reasonable for a temporary matrix to be declared const. However, I don't
  1386. know how to tell this to the C++ compiler. One possibility is to provide
  1387. alternative versions of the operators for operands declared const. But
  1388. then one gets the explosion in the number of operators.
  1389.  
  1390. My solution is to include versions of the initialisers for matrices
  1391. declared const. Otherwise, you need to use A.c() in place of A if A is
  1392. declared const and you wish to use it in an expression.
  1393.  
  1394.  
  1395. A calculus of matrix types
  1396. --------------------------
  1397.  
  1398. The program needs to be able to work out the class of the result of a
  1399. matrix expression. This is to check that a conversion is legal or to
  1400. determine the class of a temporary. To assist with this, a class
  1401. MatrixType is defined. Operators +, -, *, >= are defined to calculate
  1402. the types of the results of expressions or to check that conversions are
  1403. legal.
  1404.  
  1405.  
  1406. Error handling
  1407. --------------
  1408.  
  1409. The package does not have graceful exit from errors. All errors are
  1410. treated as fatal. Originally I thought I would wait until exceptions
  1411. became available in C++. This now seems to have been delayed. In any
  1412. case I don't think exceptions will solve all the problems. Some clean up
  1413. of objects on the heap will often be required before one can exit via an
  1414. exception.
  1415.  
  1416. There are four categories of errors:
  1417.  
  1418.    Programming error - eg illegal conversion of a matrix, subscript out
  1419.    of bounds, matrix dimensions don't match;
  1420.  
  1421.    Illegal data error - eg Cholesky of a non-positive definite matrix;
  1422.  
  1423.    Out of space error - "new" returns a null pointer;
  1424.  
  1425.    Convergence failure - an iterative operation fails to converge.
  1426.  
  1427. For the first two of these, it might be sensible to terminate a program.
  1428. For the second two, one does want to return control to the user's
  1429. program in a convenient manner. I don't know a good way of doing this,
  1430. especially before exceptions are implemented.
  1431.  
  1432.  
  1433. Band and sparse matrices
  1434. ------------------------
  1435.  
  1436. The package does not yet support band or sparse types. At present the
  1437. package assumes that the structure of a matrix is determined by its
  1438. class and dimensions. This is not sufficient for band and sparse
  1439. matrices.
  1440.  
  1441. For band matrices one also needs to know the upper and lower band
  1442. widths. For sparse matrices there is going to be some kind of structure
  1443. vector. These are going to have to be calculated for the results of
  1444. expressions in much the same way that types are calculated. In addition,
  1445. a whole new set of row and column operations would have to be written
  1446. for sparse matrices. However the present ones will be fine for band
  1447. matrices.
  1448.  
  1449. Band and sparse matrices are important for people solving large
  1450. sets of differential equations. Sparse matrices are also important for
  1451. statistical and operational research applications.
  1452.  
  1453.  
  1454. -------------------------------------------------------------------------------
  1455.  
  1456.  
  1457.                    Matrix package problem report form
  1458.                    ----------------------------------
  1459.  
  1460. Version: ...............newmat03
  1461. Date of release: .......Nov 25th, 1991
  1462. Primary site: ..........Simtel20
  1463. Downloaded from: .......
  1464. Your email address: ....
  1465. Today's date: ..........
  1466. Your machine: ..........
  1467. Compiler & version: ....
  1468. Describe the problem - attach examples if possible:
  1469.  
  1470.  
  1471.  
  1472.  
  1473.  
  1474.  
  1475.  
  1476.  
  1477.  
  1478. Email to  robert@am.dsir.govt.nz  or  Compuserve 72777,656 
  1479.  
  1480. -------------------------------------------------------------------------------
  1481.