Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

LinAlg.h

Go to the documentation of this file.
00001 // This may look like C code, but it is really -*- C++ -*-
00002 /*
00003  ************************************************************************
00004  *
00005  *                        Linear Algebra Package
00006  *
00007  * The present package implements all the basic algorithms dealing
00008  * with vectors, matrices, matrix columns, etc.
00009  * Matrix is a basic object in the package; vectors, symmetric matrices,
00010  * etc. are considered matrices of a special type.
00011  *
00012  * Matrix elements are arranged in memory in a COLUMN-wise
00013  * fashion (in FORTRAN's spirit). In fact, this makes it very easy to
00014  * pass the matrices to FORTRAN procedures implementing more
00015  * elaborate algorithms.
00016  *
00017  * Matrix and vector indices always start with 1, spanning up to 
00018  * a given limit, unless specified otherwise by a user. 
00019  *
00020  * The present package provides all facilities to completely AVOID returning
00021  * matrices. If one really needs to return a matrix, return
00022  * a LazyMatrix object instead. The conversion is completely transparent
00023  * to the end user, e.g. "Matrix m = haar_matrix(5);" and _is_ efficient.
00024  *
00025  * Container and views: object of a class Matrix is the only one that
00026  * owns the data, that is, a 2D grid of values. All other classes provide
00027  * different views of this grid (as a long 1D array of a matrix stretched
00028  * column-by-column, a view of a single matrix row, column, or a diagonal,
00029  * a view of a certain matrix block, etc).
00030  *
00031  * The views are designed to be as lightweight as possible: most of
00032  * the view functionality could be inlined. That is why views do not
00033  * have any virtual functions. The view classes are supposed to be
00034  * final, in JAVA parlance. It makes little sense to inherit from them.
00035  * It is unnecessary as views are designed to be flexible and efficient,
00036  * they provide a wide variety of access to the matrix with safety and
00037  * as efficient as possible. So flexible that many former Matrix methods
00038  * can be implemented outside of the Matrix class as they no longer need
00039  * a priviledged access to Matrix internals, without any loss of performance.
00040  *
00041  * Importance of Matrix streams: a sequential view/access to a matrix or
00042  * its parts. Many of Linear Algebra algorithms actually require only
00043  * sequential access to a matrix or its rows/columns, which is simpler and
00044  * faster than a random access. That's why LinAlg stresses streams. There are
00045  * two kinds of sequential views: ElementWise are transient views that
00046  * typically exist only for the duration of a elementwise action.
00047  * That's why it doesn't require a "locking" of a matrix. On the other hand,
00048  * LAStream and the ilk are more permanent views of the matrix, with
00049  * roughly the same functionality. An application itself can traverse the
00050  * stream at its own convenience, bookmarking certain spots and possibly
00051  * returning to them later. LAStream does assert a shared lock on the matrix,
00052  * to prevent the element matrix from moving/disappearing.
00053  *
00054  * $Id: LinAlg.h,v 1.2 2004/10/02 22:23:32 maxx Exp $
00055  *
00056  ************************************************************************
00057  */
00058 
00059 #ifndef __GNUC__
00060 #pragma once
00061 #endif
00062 #ifndef _LinAlg_h
00063 #define _LinAlg_h 1
00064 
00065 #if defined(__GNUC__)
00066 #pragma interface
00067 #define __unused(X) X __attribute__((unused))
00068 //SB+:
00069 #include <iostream.h>     // When ostream is a template, can't make it opaque
00070 using namespace std;
00071 //SB-:class ostream;                    // An Opaque class
00072 #else
00073 #include <iostream.h>     // When ostream is a template, can't make it opaque
00074 using namespace std;
00075 #define __unused(X) X
00076 #endif
00077 
00078 #include "myenv.h"
00079 #include "std.h"
00080 #if defined(__GNUC__)
00081 #include <values.h>
00082 #ifndef MAXINT
00083 #include <limits.h>
00084 #define MAXINT INT_MAX
00085 #endif
00086 #else
00087 #include <limits.h>
00088 #define MAXINT INT_MAX
00089 #endif
00090 #include <math.h>
00091 #include "builtin.h"
00092 #include "minmax.h"
00093 #if defined(index)
00094 #undef index
00095 #endif
00096 
00097 /*
00098  *------------------------------------------------------------------------
00099  *                      Auxiliary classes
00100  */
00101 
00102 typedef float REAL;                     // Scalar field of the Linear Vector
00103                                         // space
00104 namespace linalg 
00105 {
00106     using namespace linalg;
00107     
00108 
00109                         // Dimensions specifier of a 2D object
00110 class DimSpec
00111 {
00112 protected:
00113   int nrows;                            // No. of rows
00114   int ncols;                            // No. of columns
00115   int row_lwb;                          // Lower bound of the row index
00116   int col_lwb;                          // Lower bound of the col index
00117 
00118 public:
00119 
00120   DimSpec(const int _nrows, const int _ncols)   // Indices start with 1
00121         : nrows(_nrows), ncols(_ncols), row_lwb(1), col_lwb(1)
00122           { assert(nrows>0 && ncols>0); }
00123   DimSpec(const int _row_lwb, const int _row_upb,       // Or with low:upper
00124           const int _col_lwb, const int _col_upb)       // boundary specif.
00125         : nrows(_row_upb-_row_lwb+1), 
00126           ncols(_col_upb-_col_lwb+1),
00127           row_lwb(_row_lwb), col_lwb(_col_lwb)
00128           { assert(nrows>0 && ncols>0); }
00129   
00130                                 // Status inquires
00131   int q_row_lwb(void) const                     { return row_lwb; }
00132   int q_row_upb(void) const                     { return nrows+row_lwb-1; }
00133   int q_nrows(void) const                       { return nrows; }
00134   int q_col_lwb(void) const                     { return col_lwb; }
00135   int q_col_upb(void) const                     { return ncols+col_lwb-1; }
00136   int q_ncols(void) const                       { return ncols; }
00137   
00138   bool operator == (const DimSpec& another_dim_spec) const
00139         { return nrows == another_dim_spec.nrows &&
00140                  ncols == another_dim_spec.ncols &&
00141                  row_lwb == another_dim_spec.row_lwb &&
00142                  col_lwb == another_dim_spec.col_lwb; }
00143  
00144   friend ostream& operator << (ostream& os, const DimSpec& dimspec);
00145 };
00146 
00147                                 // A 2D index, a pair representing
00148                                 // a row,col location in some 2D structure
00149 struct rowcol
00150 {
00151   int row, col;
00152   rowcol(const int _row, const int _col) : row(_row), col(_col) {}
00153   bool operator == (const rowcol another) const
00154         { return row == another.row && col == another.col; }
00155   friend ostream& operator << (ostream& os, const rowcol& rc);
00156 };
00157 
00158                 // A Range of indices: lwb:upb
00159                 // if lwb > upb, the range is "empty"
00160                 // If lwb is not specified, it defaults to -IRange::INF
00161                 // (meaning the range spans from the earliest possible
00162                 // element in a collection the range is being applied to)
00163                 // If upb is not specified, it defaults to IRange::INF
00164                 // meaning the range spans through the end of the collection,
00165                 // whatever that may be.
00166 struct IRange
00167 {
00168   const int lwb, upb;
00169   enum { INF=MAXINT };
00170   IRange(const int _lwb, const int _upb) : lwb(_lwb), upb(_upb) {}
00171   static IRange from(const int lwb) { return IRange(lwb,INF); }
00172   static IRange through(const int upb) { return IRange(-INF,upb); }
00173   bool q_empty(void) const      { return lwb > upb; }
00174   bool q_proper(void) const
00175     // ms: cast for gcc 3.0.4 (IRIX)
00176         { return fabs((double)lwb) != INF && fabs((double)upb) != INF && lwb <= upb; }
00177   friend ostream& operator << (ostream& os, const IRange range);
00178   int no_less_than(const int strict_lwb) const
00179         { if( lwb == -IRange::INF ) return strict_lwb;
00180           assert(lwb >= strict_lwb); return lwb; }
00181   int no_more_than(const int strict_upb) const
00182         { if( upb == IRange::INF ) return strict_upb;
00183           assert(upb <= strict_upb); return upb; }
00184 };
00185 
00186                 // The same as above, but with a given stride
00187 struct IRangeStride
00188 {
00189   const int lwb, upb, stride;
00190   IRangeStride(const int _lwb, const int _upb, const int _stride)
00191         : lwb(_lwb), upb(_upb), stride(_stride) {}
00192   IRangeStride(const IRange range)
00193         : lwb(range.lwb), upb(range.upb), stride(1) {}
00194   bool q_empty(void) const      { return lwb > upb; }
00195   bool q_proper(void) const
00196     // ms: cast for gcc 3.0.4 (IRIX)
00197         { return fabs((double)lwb) != IRange::INF && fabs((double)upb) != IRange::INF
00198                 && lwb <= upb && stride > 0; }
00199 };
00200 
00201 
00202                 // A cut from DimSpec in two dimensions
00203 class DimSpecSubranged : public DimSpec
00204 {
00205 protected:
00206    int min_offset;      // An offset to a[imin,jmin] from the original DimSpec
00207    int max_offset;      // An offset to a[imax,jmax] from the original DimSpec
00208    const int original_nrows;    // The number of rows in the original DimSpec
00209 
00210 public:
00211  DimSpecSubranged(const DimSpec& orig, const IRange row_range,
00212                   const IRange col_range)
00213         : DimSpec(row_range.no_less_than(orig.q_row_lwb()),
00214                   row_range.no_more_than(orig.q_row_upb()),
00215                   col_range.no_less_than(orig.q_col_lwb()),
00216                   col_range.no_more_than(orig.q_col_upb())),
00217           original_nrows(orig.q_nrows())
00218           { min_offset = (col_lwb - orig.q_col_lwb())*original_nrows +
00219                          (row_lwb - orig.q_row_lwb());
00220             max_offset = min_offset + (ncols-1)*original_nrows + nrows - 1;
00221           }
00222                 // we are sure row_lwb is within [orig.row_lwb,orig.row_upb]
00223                 // and so is the new row_upb (and the same for col)
00224 
00225   int q_min_offset(void) const          { return min_offset; }
00226   int q_max_offset(void) const          { return max_offset; }
00227   int q_row_diff(void) const            { return original_nrows - nrows; }
00228 };
00229                                 // An interface describing an operation to apply
00230                                 // to a matrix element.
00231                                 // This interface is used by Matrix::apply()
00232                                 // and similar for-each-type iterators. They
00233                                 // call a concrete instance of this interface
00234                                 // on each element of a matrix etc. collection
00235                                 // under consideration.
00236                                 // The operation receives the (i,j) index of
00237                                 // the current element matrix and its _value_:
00238                                 // thus this interface may not modify the
00239                                 // element itself.
00240 class ConstElementAction
00241 {
00242 public:
00243   virtual void operation(const REAL element, const int i, const int j) = 0;
00244 };
00245 
00246                                 // The same as above, only an operation
00247                                 // is allowed to modify a matrix element
00248                                 // it is being applied to.
00249 class ElementAction
00250 {
00251 public:
00252   virtual void operation(REAL& element, const int i, const int j) = 0;
00253 
00254 //private:                      // Those aren't implemented; but making them
00255                                 // private forbids the assignement
00256 //  ElementAction(const ElementAction&);
00257 //  void operator= (const ElementAction&);
00258 };
00259 
00260 
00261                                 // An interface that performs an operation
00262                                 // on two elements of Matrices being jointly
00263                                 // traversed. The operation receives the values
00264                                 // of the two current elements in both matrices
00265                                 // and the element's index (which must be
00266                                 // the same for each matrix involved).
00267                                 // As interface recieves elements' _values_
00268                                 // it may not modify the elements.
00269                                 // This interface is being used to compute
00270                                 // some property of two matrices being jointly
00271                                 // traversed.
00272 class ConstElementAction2
00273 {
00274 public:
00275   virtual void operation(const REAL element1, const REAL element2,
00276                          const int i, const int j) = 0;
00277 };
00278 
00279                                 // Lazy matrix constructor
00280                                 // That is, a promise of a matrix rather than
00281                                 // a matrix itself. The promise is forced
00282                                 // by a Matrix constructor or assignment
00283                                 // operator, see below.
00284                                 // It's highly recommended a function never
00285                                 // return Matrix itself. Use LazyMatrix
00286                                 // instead
00287 class LazyMatrix : public DimSpec
00288 {
00289   friend class Matrix;
00290   virtual void fill_in(Matrix& m) const = 0;
00291 
00292                                 // Not implemented: cloning is prohibited
00293   //  LazyMatrix(const LazyMatrix&);
00294   LazyMatrix& operator = (const LazyMatrix&);
00295 public:
00296   LazyMatrix(const int nrows=1, const int ncols=1)      // Indices start with 1
00297     : DimSpec(nrows,ncols) {}
00298   LazyMatrix(const int row_lwb, const int row_upb,// Or with low:upper
00299              const int col_lwb, const int col_upb)// boundary specif.
00300     : DimSpec(row_lwb,row_upb,col_lwb,col_upb) {}
00301   LazyMatrix(const DimSpec& dims) : DimSpec(dims) {}
00302 };
00303 
00304                                 // A pair specifying extreme values
00305                                 // of some domain/result
00306 class MinMax
00307 {
00308   REAL min_val, max_val;                // Invariant: max_val >= min_val
00309 public:
00310   MinMax(REAL _min_val, REAL _max_val)
00311     : min_val(_min_val), max_val(_max_val) {}
00312   explicit MinMax(REAL val)
00313     : min_val(val), max_val(val) {}
00314   MinMax& operator << (REAL val)
00315         { if( val > max_val ) max_val = val;
00316           else if(val < min_val) min_val = val; return *this; }
00317   REAL min(void) const          { return min_val; }
00318   REAL max(void) const          { return max_val; }
00319   double ratio(void) const      { return max_val/min_val; }
00320   friend ostream& operator << (ostream& os, const MinMax& mx);
00321 };
00322 
00323                 // The following class watches over exclusive (write)
00324                 // and shared (read) access to a data structure the object
00325                 // is a member of.
00326                 // The ref_count field tells how the object is being
00327                 // currently accessed:
00328                 //      = 0: the object is not engaged
00329                 //      =-1: exclusive (write access)
00330                 //      >0 : the object is being currently read (viewed)
00331                 //      by one or several readers
00332                 // The watchdog watches over obvious blunders as trying
00333                 // to get a read access while the object is being exclusively
00334                 // held, or trying to grab the write access or destroy the
00335                 // object while it is engaged in some other activity.
00336                 // Note, if an object contains a pointer to some data
00337                 // structure, using this pointer requires a read (shared)
00338                 // access, even if the data structure would be modified.
00339                 // Indeed, as long as all 'readers' access the data structure
00340                 // at the same place in memory, they certainly are manipulating
00341                 // the same (and the real) thing. Modifying a data structure
00342                 // _pointer_ on the other hand requires an exclusive access
00343                 // to the object (pointer holder). Thus read/write access
00344                 // is meant a shared/exclusive access to
00345                 // object's data and object's layout (but not the contents
00346                 // of the data itself)
00347 class RWWatchDog
00348 {
00349   int ref_count;                        // see explanation above
00350 
00351   RWWatchDog(const RWWatchDog&);        // Deliberately unimplemented:
00352   void operator = (const RWWatchDog&);  // no cloning allowed!
00353 
00354   volatile void access_violation(const char reason []);
00355 
00356 public:
00357   bool q_engaged(void) const            { return ref_count != 0; }
00358   bool q_shared(void) const             { return ref_count > 0; }
00359   bool q_exclusive(void) const          { return ref_count == -1; }
00360 
00361   void get_exclusive(void)
00362         { if( q_engaged() ) access_violation("get_exclusive"); ref_count = -1; }
00363   void release_exclusive(void)
00364         { if(!q_exclusive()) access_violation("release_exclusive");
00365           ref_count = 0; }
00366 
00367   void get_shared(void)
00368         { if( q_exclusive() ) access_violation("get_shared"); ref_count++; }
00369   void release_shared(void)
00370         { if( !q_shared() ) access_violation("release_shared"); ref_count--; }
00371 
00372   RWWatchDog(void) : ref_count(0) {}
00373   ~RWWatchDog(void)
00374         { if( ref_count!=0 ) access_violation("destroying"); 
00375           ref_count = -1; }
00376 
00377                                 // Dump the current status
00378   friend ostream& operator << (ostream& os, const RWWatchDog& wd);
00379 };
00380 
00381 #if 0
00382 class MinMaxLoc : public MinMax
00383 {
00384 };
00385 #endif
00386 
00387 
00388 /*
00389  *------------------------------------------------------------------------
00390  *      Special kinds of matrices that have to be forward-declared
00391  */
00392  
00393 class Matrix;
00394 class MatrixColumn;
00395 class ConstMatrixColumn;
00396 class MatrixRow;
00397 class MatrixDiag;
00398                                 // Create an orthogonal (2^n)*(no_cols) Haar
00399                                 // (sub)matrix, whose columns are Haar
00400                                 // functions
00401                                 // If no_cols is 0, create the complete matrix
00402                                 // with 2^n columns
00403 class haar_matrix : public LazyMatrix
00404 {
00405 public:
00406   void fill_in(Matrix& m) const;
00407   haar_matrix(const int n, const int no_cols=0);
00408 };
00409 
00410                                 // Lazy transposition
00411 class LazyTransposedMatrix : public LazyMatrix
00412 {
00413   const Matrix& proto;
00414   void fill_in(Matrix& m) const;
00415   inline LazyTransposedMatrix(const Matrix & _proto);
00416 public:
00417   friend inline const LazyTransposedMatrix transposed(const Matrix & proto);
00418 };
00419 
00420 /*
00421  *------------------------------------------------------------------------
00422  *      Data classes: Matrix and (its particular case) Vector
00423  *
00424  * Note that view classes (below) may "borrow" matrix data area
00425  * for a while; to make sure they "return" the data, we use a limited form
00426  * of reference counting. That is, the destruction/reallocation of a data
00427  * class are only possible when the reference counter is 1
00428  * (that is, there are no oustanding views). Otherwise it is a grave error.
00429  *
00430  */
00431 
00432 class Matrix : public DimSpec
00433 {
00434 public:
00435   friend class ElementWiseConst;
00436   friend class ElementWiseStrideConst;
00437   friend class LAStreamIn;
00438   friend class LAStreamOut;
00439   friend class LAStrideStreamIn;
00440   friend class LAStrideStreamOut;
00441   friend class LABlockStreamIn;
00442   friend class LABlockStreamOut;
00443   friend class MatrixDABase;
00444   class ConstReference;                 // Forward decl of nested classes
00445   class Reference;
00446   friend class ConstReference;
00447   
00448   friend class Vector;
00449   friend class ConstMatrixColumn;
00450   friend class ConstMatrixRow;
00451   friend class ConstMatrixDiag;
00452 
00453 private:                        // Private part
00454   int valid_code;                       // Validation code
00455   enum { MATRIX_val_code = 5757 };      // Matrix validation code value
00456   RWWatchDog ref_counter;
00457 
00458 protected:                      // May be used in derived classes
00459   const char * name;                            // Name for the matrix
00460   int nelems;                           // Total no of elements, nrows*ncols
00461   REAL * elements;                      // Elements themselves
00462 
00463   void allocate(void);
00464 
00465   friend void haar_matrix::fill_in(Matrix& m) const;
00466   friend void LazyTransposedMatrix::fill_in(Matrix& m) const;
00467 
00468 public:                 // Public interface
00469 
00470                                 // Constructors and destructors
00471                                         // Make a blank matrix
00472   Matrix(const int nrows, const int ncols);     // Indices start with 1
00473   Matrix(const int row_lwb, const int row_upb,  // Or with low:upper
00474          const int col_lwb, const int col_upb); // boundary specif.
00475   Matrix(const DimSpec& dimension_specs);
00476   inline Matrix(const Matrix&  another); // A real copy constructor, expensive
00477 
00478 
00479                                         // Construct a matrix applying a spec
00480                                         // operation to two prototypes
00481                                         // Example: Matrix A(10,12), B(12,5);
00482                                         // Matrix C(A,Matrix::Mult,B);
00483 //  enum MATRIX_CREATORS_2op { Mult,            // A*B
00484 //                           TransposeMult,     // A'*B
00485 //                           InvMult,           // A^(-1)*B
00486 //                           AtBA };            // A'*B*A
00487 //  Matrix(const Matrix& A, const MATRIX_CREATORS_2op op, const Matrix& B);
00488 //  Matrix(const Vector& x, const Vector& y);   // x'*y (diad) matrix
00489 
00490   Matrix(const LazyMatrix& lazy_constructor);//Make a matrix using given recipe
00491   explicit Matrix(const char * file_name);      // Read the matrix from the file
00492                                         // (not yet implemented!)
00493 
00494   ~Matrix(void);                        // Destructor
00495 
00496   void set_name(const char * name);     // Set a new matrix name
00497 
00498                                         // Erase the old matrix and create a
00499                                         // new one according to new boundaries
00500   void resize_to(const int nrows, const int ncols);     // Indexation @ 1
00501   void resize_to(const int row_lwb, const int row_upb,  // Or with low:upper
00502                  const int col_lwb, const int col_upb); // boundary specif.
00503   void resize_to(const DimSpec& dimension_specs);       // Like other matrix m
00504 
00505 
00506   void is_valid(void) const
00507   { if( valid_code != MATRIX_val_code ) invalid_matrix(); }
00508 
00509                                 // Status inquires
00510   int q_no_elems(void) const            { return nelems; }
00511   const char * q_name(void) const       { return name; }
00512 
00513   class ConstReference
00514   {
00515     friend class Reference;
00516     Matrix& m;
00517     ConstReference(const ConstReference&);      // Not implemented and forbidden
00518     void operator = (const ConstReference&);    // reference isn't transferable
00519   public:
00520     ConstReference(const Matrix& _m)
00521         : m(const_cast<Matrix&>(_m))
00522         { m.is_valid(); m.ref_counter.get_shared(); }
00523    ~ConstReference(void) { m.is_valid(); m.ref_counter.release_shared(); }
00524     /*explicit*/ operator const Matrix& (void) const { return m; }
00525     const Matrix& ref(void) const { return m; }
00526   };
00527 
00528   class Reference : public ConstReference
00529   {
00530     Reference(const Reference&);        // Not implemented and forbidden
00531     void operator = (const Reference&); // reference isn't transferable
00532   public:
00533     Reference(Matrix& _m) : ConstReference(_m) {}
00534     /*explicit*/ operator Matrix& (void) const { return m; }
00535     Matrix& ref(void) const { return m; }
00536   };
00537 
00538                                         // Apply a user-defined action
00539   Matrix& apply(ElementAction& action);         // to each matrix element
00540   const Matrix& apply(ConstElementAction& action) const; // to each matrix element
00541 
00542                                         // Although the following method does
00543                                         // not need a priviledged access to a
00544                                         // Matrix, it is required to be a member
00545                                         // of the class
00546   inline Matrix& operator = (const REAL val);
00547                                         // Invert the matrix returning the
00548                                         // determinant if desired
00549                                         // determinant = 0 if the matrix is
00550                                         // singular
00551                                         // If determ_ptr=dummy_determinant_ref
00552                                         // and the matrix *is*
00553                                         // singular, throw up
00554   static double dummy_determinant_ref;
00555   Matrix& invert(double &determ_ref = dummy_determinant_ref);
00556 
00557                                 // Element-wise operations on two matrices
00558   inline Matrix& operator = (const Matrix& source);     // Assignment
00559   inline Matrix& operator = (const ConstReference& ref)
00560                 { return *this = ref.ref(); }
00561   Matrix& operator = (const LazyMatrix& source);// Force the promise of a matrix
00562   Matrix& clear(void);                  // Clear the matrix (fill out with 0)
00563  
00564 
00565                                         // Comparisons
00566   bool operator == (const Matrix& im2) const;
00567   friend inline void are_compatible(const Matrix& im1, const Matrix& im2);
00568 
00569 
00570                                 // True matrix operations
00571                                 // (on a matrix as a whole)  
00572   Matrix& operator *= (const Matrix& source);   // Inplace multiplication
00573                                                 // possible only for square src
00574   inline Matrix& operator *= (const ConstReference& ref)
00575                 { return *this *= ref.ref(); }
00576 
00577   
00578   void mult(const Matrix& A, const Matrix& B);  // Compute A*B
00579  
00580                                 // Compute matrix norms
00581   double row_norm(void) const;          // MAX{ SUM{ |M(i,j)|, over j}, over i}
00582   double norm_inf(void) const           // Alias, shows the norm is induced
00583          { return row_norm(); }         //      by the vector inf-norm
00584   double col_norm(void) const;          // MAX{ SUM{ |M(i,j)|, over i}, over j}
00585   double norm_1(void) const             // Alias, shows the norm is induced
00586          { return col_norm(); }         //      by the vector 1-norm
00587   double e2_norm(void) const;           // SUM{ m(i,j)^2 }, Note it's square
00588                                         // of the Frobenius rather than 2. norm
00589 
00590   friend double e2_norm(const Matrix& m1, const Matrix& m2);
00591                                         // e2_norm(m1-m2)
00592 
00593   double determinant(void) const;       // Matrix must be a square one
00594 
00595   double asymmetry_index(void) const;   // How far is the matrix from being
00596                                         // symmetrical (0 means complete symm)
00597                                         // (not yet implemented)
00598 
00599                                 // Make matrices of a special kind
00600   Matrix& unit_matrix(void);            // Matrix needn't be a square
00601   Matrix& hilbert_matrix(void);         // Hilb[i,j] = 1/(i+j-1); i,j=1..max
00602 
00603                                 // I/O: write, read, print 
00604                                         // Write to a file
00605                                         // "| command name" is OK as a file
00606                                         // name
00607   void write(const char * file_name,const char * title = "") const;
00608   void info(void) const;                // Print the info about the Matrix
00609   void print(const char title []) const;// Print the Matrix as a table
00610   volatile void invalid_matrix(void) const;   // The matrix in question is not valid
00611                                         // ==> crash the program
00612 };
00613 
00614 void compare(const Matrix& im1, const Matrix& im2, const char title[]);
00615 
00616 
00617 /*
00618  *------------------------------------------------------------------------
00619  *                 Vector as a n*1 matrix (that is, a col-matrix)
00620  */
00621 class Vector : public Matrix
00622 {
00623 public:
00624   Vector(const int n);          // Create a blank vector of a given size
00625                                 // (indexation starts at 1)
00626                                 // Create a general lwb:upb vector blank vector
00627   Vector(const int lwb, const int upb);
00628 //  Vector(const Vector& another);      // a copy constructor (to be inferred)
00629 
00630                                         // Make a vector and assign init vals
00631   Vector(const int lwb, const int upb,  // lower and upper bounds
00632          double iv1, ...                // DOUBLE numbers. The last arg of
00633          );                             // the list must be string "END"
00634                                         // E.g: Vector f(1,2,0.0,1.5,"END");
00635 
00636   Vector(const LazyMatrix& lazy_constructor);//Make a vector using given recipe
00637 
00638                                 // Resize the vector (keeping as many old
00639                                 // elements as possible), expand by zeros
00640   void resize_to(const int n);                  // Indexation @ 1
00641   void resize_to(const int lwb, const int upb);         // lwb:upb specif
00642   void resize_to(const Vector& v);                      // like other vector
00643 
00644                                 // Unlike Matrix, Vector permits a direct
00645                                 // access to its elements without further ado
00646                                 // Still, streams/of_every() are more efficient
00647                                 // for uniform access
00648   inline REAL& operator () (const int index);
00649   inline REAL  operator () (const int index) const
00650         { return (const_cast<Vector&>(*this))(index); }
00651 
00652                                 // Listed below are specific vector operations
00653                                 // (unlike n*1 matrices)
00654 
00655                                         // Status inquires
00656   int q_lwb(void) const         { return q_row_lwb(); }
00657   int q_upb(void) const         { return q_row_upb(); }
00658 
00659                                         // Compute the scalar product
00660   friend double operator * (const Vector& v1, const Vector& v2);
00661 
00662                                         // "Inplace" multiplication
00663                                         // target = A*target
00664                                         // A needn't be a square one (the
00665                                         // target will be resized to fit)
00666   Vector& operator *= (const Matrix& A);
00667 
00668                                         // Vector norms
00669   inline double norm_1(void) const;             // SUM{ |v[i]| }
00670   inline double norm_2_sqr(void) const;         // SUM{ v[i]^2 }
00671   inline double norm_inf(void) const;           // MAX{ |v[i]| }
00672 
00673   Vector& operator = (const Vector& v)
00674         { Matrix::operator =(v); return *this; }
00675   Vector& operator = (const LazyMatrix& source);// Force the promise of a vector
00676 
00677                                         // Although the following methods do
00678                                         // not need a priviledged access to
00679                                         // Vector, they are required to be
00680                                         // members of the class
00681   inline Vector& operator = (const ElementWiseConst& stream);
00682   inline Vector& operator = (const ElementWiseStrideConst& stream);
00683   inline Vector& operator = (const REAL val)
00684         { Matrix::operator =(val); return *this; }
00685 };
00686 
00687                                 // Service functions (useful in the
00688                                 // verification code). They print some detail
00689                                 // info if the validation condition fails
00690 void verify_element_value(const Matrix& m,const REAL val);
00691 void verify_matrix_identity(const Matrix& m1, const Matrix& m2);
00692 inline void verify_element_value(const Matrix::ConstReference& m,const REAL val)
00693         { verify_element_value(m.ref(),val); }
00694 inline void verify_matrix_identity(const Matrix::ConstReference& m1,
00695                                    const Matrix& m2)
00696         { verify_matrix_identity(m1.ref(),m2); }
00697 inline void verify_matrix_identity(const Matrix& m1,
00698                                    const Matrix::ConstReference& m2)
00699         { verify_matrix_identity(m1,m2.ref()); }
00700 inline void verify_matrix_identity(const Matrix::ConstReference& m1,
00701                                    const Matrix::ConstReference& m2)
00702         { verify_matrix_identity(m1.ref(),m2.ref()); }
00703 
00704                                 // Multiply by the diagonal
00705                                 // of another matrix
00706 Matrix& operator *= (Matrix& m, const ConstMatrixDiag& diag);
00707 
00708                                 // Compute the product of two matrices
00709 class LazyMatrixProduct : public linalg::LazyMatrix
00710 {
00711 private:
00712   const Matrix& A;
00713   const Matrix& B;
00714   void fill_in(Matrix& m) const { m.mult(A,B); };
00715     
00716   LazyMatrixProduct(const Matrix & __A, 
00717                     const Matrix & __B)
00718         : LazyMatrix(__A.q_row_lwb(),__A.q_row_upb(),
00719                      __B.q_col_lwb(),__B.q_col_upb()),
00720           A(__A), B(__B) {};
00721     
00722   
00723 public:
00724   friend const LazyMatrixProduct operator * (const Matrix& A, const Matrix& B)
00725         { return LazyMatrixProduct(A,B); };
00726     
00727 };
00728 
00729                                 // Make a zero matrix
00730 class LazyZeroMatrix : public LazyMatrix
00731 {
00732   void fill_in(Matrix& m) const { __unused(Matrix& m1) = m; }
00733   LazyZeroMatrix(const Matrix& proto)
00734     : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
00735 public:
00736   friend const LazyZeroMatrix zero(const Matrix& proto)
00737     { proto.is_valid(); return proto; }
00738 };
00739 
00740                                 // Make a unit matrix
00741 class LazyUnitMatrix : public LazyMatrix
00742 {
00743   void fill_in(Matrix& m) const { m.unit_matrix(); }
00744   LazyUnitMatrix(const Matrix& proto)
00745     : LazyMatrix(static_cast<const DimSpec&>(proto)) {}
00746 public:
00747   friend const LazyUnitMatrix unit(const Matrix& proto)
00748     { proto.is_valid(); return proto; }
00749 };
00750 
00751                                 // Make an inverse matrix
00752 class LazyInverseMatrix : public LazyMatrix
00753 {
00754   const Matrix& proto;
00755   void fill_in(Matrix& m) const { m = proto; m.invert(); }
00756   LazyInverseMatrix(const Matrix& _proto)
00757     : LazyMatrix(static_cast<const DimSpec&>(_proto)),
00758       proto(_proto) {}
00759 public:
00760   friend const LazyInverseMatrix inverse(const Matrix& proto)
00761     { proto.is_valid(); return proto; }
00762 };
00763 
00764 
00765                                 // Make a transposed matrix
00766 inline LazyTransposedMatrix::LazyTransposedMatrix(const Matrix & _proto)
00767         : LazyMatrix(_proto.q_col_lwb(),_proto.q_col_upb(),
00768                      _proto.q_row_lwb(),_proto.q_row_upb()),
00769           proto(_proto) {}
00770 inline const LazyTransposedMatrix transposed(const Matrix & proto)
00771         { proto.is_valid(); return proto; }
00772 
00773 /*
00774  *------------------------------------------------------------------------
00775  *      Matrix Views: provide access to elements of a Matrix
00776  *      (or groups of elements: rows, columns, diagonals, etc)
00777  * Although views can be created and exist for some time, one should
00778  * not attempt to pass around views (especially return them as a result
00779  * of the function) when the data object itself (a matrix) could be
00780  * disposed of. Thanks to the reference counting, this situation would
00781  * be detected, and the program would crash.
00782  */
00783  
00784 /*
00785  *------------------------------------------------------------------------
00786  */
00787 
00788 /*
00789  * Functions/procedures (generally grouped together under a ElementWise
00790  * class umbrella) that perform a variety of operation on each element
00791  * of some structure in turn.
00792  *
00793  * ElementWiseConst and ElementWise classes are largely a syntactic
00794  * sugar; you can't explicitly construct an object of this class, on stack
00795  * or on heap. The only way the class can be used is within a pattern
00796  *      to_every(matrix) = 1;
00797  *      if( of_every(matrix) > 1 ) ...
00798  * etc. Besides, this is the only way it makes sense.
00799  *
00800  * Note that the ElementWise classes are friends of Matrix, so they are
00801  * privy of Matrix internals and can use them to make processing faster.
00802  */
00803 
00804 
00805                                 // A functor describing an operation to apply 
00806                                 // to every element of a collection in
00807                                 // turn
00808                                 // This operation receives only the elements's
00809                                 // value, and thus it may not modify the
00810                                 // element under consideration.
00811 class ElementWiseConstAction
00812 {
00813                                 // Those aren't implemented; but making them
00814                                 // private forbids the assignement
00815   ElementWiseConstAction(const ElementWiseConstAction&);
00816   ElementWiseConstAction& operator= (const ElementWiseConstAction&);
00817 public:
00818   ElementWiseConstAction(void) {}
00819   virtual void operator () (const REAL element) = 0;
00820 };
00821 
00822                                 // The same as above but this functor is given
00823                                 // a reference to an element, which it can
00824                                 // modify
00825 class ElementWiseAction
00826 {
00827                                 // Those aren't implemented; but making them
00828                                 // private forbids the assignement
00829   ElementWiseAction(const ElementWiseAction&);
00830   ElementWiseAction& operator= (const ElementWiseAction&);
00831 public:
00832   ElementWiseAction(void) {}
00833   virtual void operator () (REAL& element) = 0;
00834 };
00835 
00836 class ElementWiseConst
00837 {
00838 public:
00839   friend class ElementWise;
00840   friend class ElementWiseStrideConst;
00841   friend class ElementWiseStride;
00842 
00843   REAL * const start_ptr;               // Pointer to the beginning of
00844                                         // the group of elements
00845   REAL * const end_ptr;                 // Points after the end of the group
00846   
00847   void operator=(const ElementWiseConst&);// is not implemented, ergo, is
00848                                         // not allowed
00849   ElementWiseConst(const ElementWiseConst&);// Cloning is not allowed, either
00850   
00851                                 // A private constructor, to make
00852                                 // sure the object can't be constructed
00853                                 // and left, while the matrix would disappear
00854                                 
00855                                 // Use of_every() friend to do
00856                                 // the actual construction
00857   ElementWiseConst(const Matrix& m)
00858         : start_ptr(const_cast<REAL*>(m.elements)),
00859           end_ptr(const_cast<REAL*>(m.elements+m.nelems))
00860         { m.is_valid(); }
00861 
00862   inline ElementWiseConst(const ConstMatrixColumn& mc);
00863   
00864                                 // Just a helper class to disambiguate 
00865                                 // type conversion from ElementWiseConst
00866                                 // to ElementWiseStrideConst and others
00867   struct ElementWiseConstRef
00868   {
00869     const ElementWiseConst& ref;
00870     ElementWiseConstRef(const ElementWiseConstRef&);
00871     void operator = (const ElementWiseConstRef&);
00872     ElementWiseConstRef(const ElementWiseConst& _ref) : ref(_ref) {}
00873   };
00874 
00875 public:
00876 
00877                                         // No public constructors...
00878 
00879                                         // ElementWiseConst(Matrix& m) would be
00880                                         // implicitly called
00881   friend inline ElementWiseConst of_every(const Matrix& m) { return m; }
00882   friend inline ElementWiseConst of_every(const Matrix::ConstReference& m)
00883         { return m.ref(); }
00884   
00885   friend inline ElementWiseConst of_every(const ConstMatrixColumn& mc);
00886                                 
00887                                 // Every ElementWiseConst is a
00888                                 // ElementWiseStrideConst with a stride of 1
00889   inline operator ElementWiseStrideConst (void);
00890   
00891                                 // Comparisons
00892                                 // Find out if the predicate
00893                                 // "element op val" is true for ALL
00894                                 // elements of the group
00895   bool    operator ==  (const REAL val) const;  // ? all elems == val
00896   bool    operator !=  (const REAL val) const;  // ? all elems != val
00897   bool    operator <   (const REAL val) const;  // ? all elems <  val
00898   bool    operator <=  (const REAL val) const;  // ? all elems <= val
00899   bool    operator >   (const REAL val) const;  // ? all elems >  val
00900   bool    operator >=  (const REAL val) const;  // ? all elems >= val
00901   
00902                                 // Find out if the predicate
00903                                 // "element op another.elem" is true for ALL
00904                                 // elements of the two groups
00905   bool    operator ==  (const ElementWiseConst& another) const;
00906   bool    operator !=  (const ElementWiseConst& another) const;
00907   bool    operator <   (const ElementWiseConst& another) const;
00908   bool    operator <=  (const ElementWiseConst& another) const;
00909   bool    operator >   (const ElementWiseConst& another) const;
00910   bool    operator >=  (const ElementWiseConst& another) const;
00911                                 
00912   void sure_compatible_with(const ElementWiseConst& another) const
00913         { assure(end_ptr-start_ptr == another.end_ptr-another.start_ptr,
00914                 "Two groups of elements to operate on have different sizes"); }
00915   
00916                                 // Data reduction: reduce a collection
00917                                 // to one value: a "norm" of the collection
00918   double sum(void) const;
00919   double sum_squares(void) const;
00920   double sum_abs(void) const;
00921   double max_abs(void) const;
00922   //double foldr(const double seed, Reductor& reductor);
00923   
00924                                 // Reduce a difference between two collections
00925                                 // to a single value: a "norm" of the
00926                                 // difference
00927   double sum_squares(const ElementWiseConst& another) const;
00928   double sum_abs(const ElementWiseConst& another) const;
00929   double max_abs(const ElementWiseConst& another) const;
00930   //double foldr(const double seed, Reductor2& reductor);
00931   
00932                                 // Just let the user do what he wants
00933                                 // with each element in turn
00934                                 // (without modifying them, of course)
00935   ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
00936 };
00937 
00938 
00939 class ElementWise : public ElementWiseConst
00940 {
00941   void operator=(const ElementWise&);   // is not implemented, ergo, is
00942                                         // not allowed
00943   ElementWise(const ElementWise&);      // Cloning is not allowed, either
00944   
00945                                 // A private constructor, to make
00946                                 // sure the object can't be constructed
00947                                 // and left, while the matrix would disappear
00948                                 
00949                                 // Use a to_every() friend to do
00950                                 // the actual construction
00951   ElementWise(Matrix& m) : ElementWiseConst(m) {}
00952   inline ElementWise(const MatrixColumn& mc);
00953 public:
00954 
00955                                         // No public constructors...
00956 
00957                                         // ElementWise(Matrix& m) would be
00958                                         // implicitly called
00959   friend inline ElementWise to_every(Matrix& m) { return m; }
00960   friend inline ElementWise to_every(const Matrix::Reference& m)
00961         { return m.ref(); }
00962 //  friend inline ElementWiseConst of_every(const ConstMatrixColumn& mc);
00963   friend inline ElementWise to_every(const MatrixColumn& mc); // see below
00964 
00965                                 // Every ElementWise is a
00966                                 // ElementWiseStride with a stride of 1
00967   inline operator ElementWiseStride (void);
00968                                 // The same as the conversion operator above,
00969                                 // use when a compiler (gcc 2.8.1) fails to
00970                                 // apply it
00971   inline ElementWiseStride with_stride (void);
00972 
00973                                 // group-scalar arithmetics
00974                                 // Modify every element of the group
00975                                 // according to the operation
00976   void operator =   (const REAL val);   // Assignment to all the elems
00977   void operator +=  (const double val); // Add to elements
00978   void operator -=  (const double val); // Take from elements
00979   void operator *=  (const double val); // Multiply elements by a val
00980 
00981                                 // Other element-wise matrix operations
00982   void abs(void);                       // Take an absolute value of a matrix
00983   void sqr(void);                       // Square each element
00984   void sqrt(void);                      // Take the square root
00985 
00986   void operator = (const ElementWiseConst& another);
00987   void operator += (const ElementWiseConst& another);
00988   void operator -= (const ElementWiseConst& another);
00989   void operator *= (const ElementWiseConst& another);
00990   void operator /= (const ElementWiseConst& another);
00991 
00992   ElementWiseAction& apply(ElementWiseAction& functor);
00993 };
00994 
00995 inline Vector& Vector::operator = (const ElementWiseConst& stream)
00996         { to_every(*this) = stream; return *this; }
00997 
00998 
00999 inline double Vector::norm_1(void) const
01000         { return of_every(*this).sum_abs(); }
01001 inline double Vector::norm_2_sqr(void) const
01002         { return of_every(*this).sum_squares(); }
01003 inline double Vector::norm_inf(void) const
01004         { return of_every(*this).max_abs(); }
01005 
01006 //---------------------------------------------------------------------------
01007 //                      The same, only with a stride
01008 
01009 class ElementWiseStrideConst
01010 {
01011   friend class ElementWiseStride;
01012   friend class ElementWiseConst;
01013 
01014   REAL * const start_ptr;               // Pointer to the beginning of
01015                                         // the group of elements
01016   REAL * const end_ptr;                 // Points after the end of the group
01017   const int stride;
01018   
01019   void operator=(const ElementWiseStrideConst&);// is not implemented, ergo, is
01020                                         // not allowed
01021   ElementWiseStrideConst(const ElementWiseStrideConst&);// Cloning is not allowed, either
01022   
01023                                 // A private constructor, to make
01024                                 // sure the object can't be constructed
01025                                 // and left, while the matrix would disappear
01026                                 
01027                                 // Use of_every() friend to do
01028                                 // the actual construction
01029   ElementWiseStrideConst(const ElementWiseConst::ElementWiseConstRef& ewc)
01030         : start_ptr(ewc.ref.start_ptr),
01031           end_ptr(ewc.ref.end_ptr),
01032           stride(1)
01033           {}
01034           
01035   inline ElementWiseStrideConst(const ConstMatrixRow& mr);
01036   inline ElementWiseStrideConst(const ConstMatrixDiag& md);
01037 
01038 public:
01039 
01040                                         // No public constructors...
01041  
01042  friend inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr);
01043  friend inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md);
01044 
01045                                 // Comparisons
01046                                 // Find out if the predicate
01047                                 // "element op val" is true for ALL
01048                                 // elements of the group
01049   bool    operator ==  (const REAL val) const;  // ? all elems == val
01050   bool    operator !=  (const REAL val) const;  // ? all elems != val
01051   bool    operator <   (const REAL val) const;  // ? all elems <  val
01052   bool    operator <=  (const REAL val) const;  // ? all elems <= val
01053   bool    operator >   (const REAL val) const;  // ? all elems >  val
01054   bool    operator >=  (const REAL val) const;  // ? all elems >= val
01055   
01056                                 // Find out if the predicate
01057                                 // "element op another.elem" is true for ALL
01058                                 // elements of the two groups
01059   bool    operator ==  (const ElementWiseStrideConst& another) const;
01060   bool    operator !=  (const ElementWiseStrideConst& another) const;
01061   bool    operator <   (const ElementWiseStrideConst& another) const;
01062   bool    operator <=  (const ElementWiseStrideConst& another) const;
01063   bool    operator >   (const ElementWiseStrideConst& another) const;
01064   bool    operator >=  (const ElementWiseStrideConst& another) const;
01065                                 
01066 //  void sure_compatible_with(const ElementWiseStrideConst& another) const
01067 //      { assure(end_ptr-start_ptr == another.end_ptr-another.start_ptr,
01068 //              "Two groups of elements to operate on have different sizes"); }
01069   
01070                                 // Data reduction: reduce a collection
01071                                 // to one value: a "norm" of the collection
01072   double sum(void) const;
01073   double sum_squares(void) const;
01074   double sum_abs(void) const;
01075   double max_abs(void) const;
01076   //double foldr(const double seed, Reductor& reductor);
01077   
01078                                 // Reduce a difference between two collections
01079                                 // to a single value: a "norm" of the
01080                                 // difference
01081   double sum_squares(const ElementWiseStrideConst& another) const;
01082   double sum_abs(const ElementWiseStrideConst& another) const;
01083   double max_abs(const ElementWiseStrideConst& another) const;
01084   //double foldr(const double seed, Reductor2& reductor);
01085   
01086                                 // Just let the user do what he wants
01087                                 // with each element in turn
01088                                 // (without modifying them, of course)
01089   ElementWiseConstAction& apply(ElementWiseConstAction& functor) const;
01090 };
01091   
01092 class ElementWiseStride : public ElementWiseStrideConst
01093 {
01094   friend class ElementWise;
01095   void operator=(const ElementWiseStride&);     // is not implemented, ergo, is
01096                                         // not allowed
01097   ElementWiseStride(const ElementWiseStride&);  // Cloning is not allowed, either
01098   
01099                                 // A private constructor, to make
01100                                 // sure the object can't be constructed
01101                                 // and left, while the matrix would disappear
01102                                 
01103                                 // Use a to_every() friend to do
01104                                 // the actual construction
01105   inline ElementWiseStride(const MatrixRow& mr);
01106   inline ElementWiseStride(const MatrixDiag& md);
01107   ElementWiseStride(const ElementWiseConst::ElementWiseConstRef& ewc)
01108         : ElementWiseStrideConst(ewc) {}
01109 public:
01110 
01111                                         // No public constructors...
01112 
01113                                 // The return statement of the friends below
01114                                 // would implicitly call a constructor
01115   friend inline ElementWiseStride to_every(const MatrixRow& mr);
01116   friend inline ElementWiseStride to_every(const MatrixDiag& md);
01117  
01118                                 // group-scalar arithmetics
01119                                 // Modify every element of the group
01120                                 // according to the operation
01121   void operator =   (const REAL val);   // Assignment to all the elems
01122   void operator +=  (const double val); // Add to elements
01123   void operator -=  (const double val); // Take from elements
01124   void operator *=  (const double val); // Multiply elements by a val
01125 
01126                                 // Other element-wise matrix operations
01127   void abs(void);                       // Take an absolute value of a coll
01128   void sqr(void);                       // Square each element
01129   void sqrt(void);                      // Take the square root
01130 
01131   void operator = (const ElementWiseStrideConst& another);
01132   void operator += (const ElementWiseStrideConst& another);
01133   void operator -= (const ElementWiseStrideConst& another);
01134   void operator *= (const ElementWiseStrideConst& another);
01135   void operator /= (const ElementWiseStrideConst& another);
01136 
01137   ElementWiseAction& apply(ElementWiseAction& functor);
01138 };
01139 
01140 inline ElementWiseConst::operator ElementWiseStrideConst (void)
01141         { return ElementWiseConstRef(*this); }
01142 
01143 inline ElementWise::operator ElementWiseStride (void)
01144         { return ElementWiseConstRef(*this); }
01145 
01146 inline ElementWiseStride ElementWise::with_stride (void)
01147         { return ElementWiseConstRef(*this); }
01148 
01149 inline Vector& Vector::operator = (const ElementWiseStrideConst& stream)
01150         {  to_every(*this).with_stride() = stream; return *this; }
01151 
01152 
01153 /*
01154  *------------------------------------------------------------------------
01155  *      MatrixDA: a view of a matrix that provides a direct access
01156  *      to its elements. A column index is built to speed up access
01157  */
01158 
01159                                 // A mixin class that builds a row index
01160                                 // for a matrix grid
01161 class MatrixDABase : public DimSpec
01162 {
01163   REAL * const * const index;           // index[i] = &matrix(0,i) (col index)
01164   MatrixDABase(const MatrixDABase&);    // Not implemented and forbidden:
01165   void operator = (const MatrixDABase&);// no cloning/assignment allowed
01166   
01167   static REAL * const * build_index(const Matrix& m);
01168 
01169 protected:
01170   MatrixDABase(const Matrix& m)
01171         : DimSpec(m),
01172           index( build_index(m) )
01173           {}
01174   inline const REAL& operator () (const int rown, const int coln) const;
01175 public:
01176   ~MatrixDABase(void);
01177 };
01178 
01179 class ConstMatrixDA : public Matrix::ConstReference, public MatrixDABase
01180 {
01181   ConstMatrixDA(const ConstMatrixDA&);  // Not implemented and forbidden:
01182   void operator = (const ConstMatrixDA&);// no cloning/assignment allowed
01183 public:
01184   ConstMatrixDA(const Matrix& m)
01185         : Matrix::ConstReference(m), MatrixDABase(m)
01186           {}
01187                                 // Individual element manipulations
01188   const REAL operator () (const int rown, const int coln) const
01189                         { return MatrixDABase::operator() (rown,coln); }
01190   const REAL operator () (const rowcol rc) const
01191                 { return operator () (rc.row,rc.col); }
01192                                 // I wish the cast to const Matrix& worked,
01193                                 // as given by the following operator. Alas,
01194                                 // an explicit member function seems to
01195                                 // be needed    
01196 //  const Matrix& get_container(void) const
01197 //      { return Matrix::ConstReference::operator const Matrix& (); }
01198 };
01199  
01200 class MatrixDA : public Matrix::Reference, public MatrixDABase
01201 {
01202   MatrixDA(const MatrixDA&);            // Not implemented and forbidden:
01203   void operator = (const MatrixDA&);    // no cloning/assignment allowed
01204 
01205 public:
01206   MatrixDA(Matrix& m)
01207         : Matrix::Reference(m), MatrixDABase(m)
01208           {}
01209                                 // Individual element manipulations
01210   REAL& operator () (const int rown, const int coln) const
01211         { return const_cast<REAL&>(MatrixDABase::operator()(rown,coln)); }
01212   REAL& operator () (const rowcol rc) const
01213                 { return operator () (rc.row,rc.col); }
01214         
01215                                 // I wish the cast to Matrix& worked by itself,
01216                                 // as given by the following operator. Alas,
01217                                 // an explicit member function seems to
01218                                 // be needed
01219 #if defined (__GNUC__)
01220   Matrix& get_container(void) { return Matrix::Reference::operator Matrix& (); }
01221 #else
01222   Matrix& get_container(void) { return *this; }
01223 #endif
01224 };
01225  
01226 /*
01227  *------------------------------------------------------------------------
01228  *              MatrixRow, MatrixCol, MatrixDiag
01229  * They are special kind of views, of selected parts of the matrix
01230  * they are implemented as special kind of streams that walk only selected
01231  * parts of the matrix
01232  */
01233 
01234 class ConstMatrixColumn : protected Matrix::ConstReference,
01235                           public DimSpec
01236 {
01237   friend class ElementWiseConst;
01238   friend class LAStreamIn;
01239   friend class LAStreamOut;
01240   const REAL * const col_ptr;           // Pointer to the column under
01241                                         // consideration
01242   ConstMatrixColumn(const ConstMatrixColumn&);  // Not implemented and forbidden:
01243   void operator = (const ConstMatrixColumn&);// no cloning/assignment allowed
01244 
01245 protected:
01246                                 // Access the i-th element of the column
01247   inline const REAL& get_ref(const int index) const;
01248 
01249 public:                                 // Take a col of the matrix
01250   ConstMatrixColumn(const Matrix& m, const int col);
01251 
01252   const REAL operator () (const int i) const
01253                         { return get_ref(i); }
01254 };
01255  
01256 inline ElementWiseConst::ElementWiseConst(const ConstMatrixColumn& mc)
01257         : start_ptr(const_cast<REAL*>(mc.col_ptr)),
01258           end_ptr(const_cast<REAL*>(mc.col_ptr)+mc.nrows)
01259         { }
01260 
01261                                 // Access the i-th element of the column
01262 inline const REAL& ConstMatrixColumn::get_ref(const int index) const
01263 { register const int offset = index - row_lwb;
01264   if( offset >= nrows || offset < 0 )
01265       _error("MatrixColumn index %d is out of row boundaries [%d,%d]",
01266               index,q_row_lwb(),q_row_upb());
01267   return col_ptr[offset];
01268 }
01269 
01270 
01271 inline ElementWiseConst of_every(const ConstMatrixColumn& mc)
01272                 { return mc; }
01273 
01274 class MatrixColumn : public ConstMatrixColumn
01275 {
01276   friend class ElementWise;
01277   friend class LAStreamOut;
01278   void operator = (const MatrixColumn&);// no cloning/assignment allowed
01279 public:                                 // Take a col of the matrix
01280   MatrixColumn(Matrix& m, const int col) :
01281         ConstMatrixColumn(m,col) {}
01282 
01283   REAL& operator () (const int i)
01284                         { return const_cast<REAL&>(get_ref(i)); }
01285 };
01286 
01287 inline ElementWise::ElementWise(const MatrixColumn& mc) : ElementWiseConst(mc) {}
01288 inline ElementWise to_every(const MatrixColumn& mc) { return mc; }
01289 
01290                         // A view of a single row of a matrix
01291 class ConstMatrixRow : protected Matrix::ConstReference,
01292                        public DimSpec
01293 {
01294   friend class ElementWiseStrideConst;
01295   friend class LAStrideStreamIn;
01296   friend class LAStrideStreamOut;
01297   const REAL * const row_ptr;           // Pointer to the row under
01298                                         // consideration
01299   const int stride;                     // if ptr=@a[row,i], then
01300                                         //    ptr+stride = @a[row,i+1]
01301                                         // Since elements of a[] are stored
01302                                         // col after col, stride = nrows
01303   const REAL * const end_ptr;           // Points after the end of the matrix
01304   
01305   ConstMatrixRow(const ConstMatrixRow&);        // Not implemented and forbidden:
01306   void operator = (const ConstMatrixRow&);      // no cloning/assignment allowed
01307 
01308 protected:
01309                                 // Access the i-th element of the row
01310   inline const REAL& get_ref(const int index) const;
01311 
01312 public:                                 // Take a row of the matrix
01313   ConstMatrixRow(const Matrix& m, const int row);
01314 
01315   const REAL operator () (const int i) const
01316                         { return get_ref(i); }
01317 };
01318  
01319 inline
01320 ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixRow& mr)
01321         : start_ptr(const_cast<REAL*>(mr.row_ptr)),
01322           end_ptr(const_cast<REAL*>(mr.end_ptr)),
01323           stride(mr.stride)
01324         { }
01325 
01326                                 // The return statement of_every below
01327                                 // would implicitly call a constructor
01328 inline ElementWiseStrideConst of_every(const ConstMatrixRow& mr)
01329         { return mr; }
01330 
01331 
01332                                 // Access the i-th element of the row
01333 inline const REAL& ConstMatrixRow::get_ref(const int index) const
01334 { register const int offset = index - col_lwb;
01335   if( offset >= ncols || offset < 0 )
01336       _error("MatrixRow index %d is out of column boundaries [%d,%d]",
01337               index,q_col_lwb(),q_col_upb());
01338   return row_ptr[stride*offset];
01339 }
01340 
01341 
01342 class MatrixRow : public ConstMatrixRow
01343 {
01344   friend class ElementWiseStride;
01345   friend class LAStrideStreamOut;
01346   MatrixRow(const MatrixRow&);          // Not implemented and forbidden:
01347   void operator = (const MatrixRow&);   // no cloning/assignment allowed
01348 public:                                 // Take a col of the matrix
01349   MatrixRow(Matrix& m, const int row) :
01350         ConstMatrixRow(m,row) {}
01351 
01352   REAL& operator () (const int i)
01353                         { return const_cast<REAL&>(get_ref(i)); }
01354 };
01355 
01356 inline ElementWiseStride::ElementWiseStride(const MatrixRow& mr)
01357         : ElementWiseStrideConst(mr) {}
01358 
01359 inline ElementWiseStride to_every(const MatrixRow& mr) { return mr; }
01360 
01361                                 // A view of the matrix' main diagonal
01362                                 // ConstMatrixDiag is described as a Matrix
01363                                 // 1:n x 1:n, where n=min(nrows,ncols)
01364 class ConstMatrixDiag : protected Matrix::ConstReference,
01365                         public DimSpec
01366 {
01367   friend class ElementWiseStrideConst;
01368   friend class LAStrideStreamIn;
01369   friend class LAStrideStreamOut;
01370   const REAL * const start_ptr;         // Pointer to the upper left corner
01371   const int stride;                     // if ptr=@a[i,i], then
01372                                         //    ptr+stride = @a[i+1,i+1]
01373                                         // Since elements of a[] are stored
01374                                         // col after col, stride = nrows+1
01375   const REAL * const end_ptr;           // Points after the end of the matrix
01376   
01377   ConstMatrixDiag(const ConstMatrixDiag&);      // Not implemented and forbidden:
01378   void operator = (const ConstMatrixDiag&);     // no cloning/assignment allowed
01379 
01380 protected:
01381                                 // Access the i-th element of the diag
01382                                 // Note that numbering always starts with 1,
01383                                 // regardless of col_lwb and row_lwb of the
01384                                 // original matrix
01385   inline const REAL& get_ref(const int index) const;
01386 
01387 public:
01388   explicit ConstMatrixDiag(const Matrix& m);
01389 
01390   const REAL operator () (const int i) const
01391                         { return get_ref(i); }
01392 };
01393  
01394 inline
01395 ElementWiseStrideConst::ElementWiseStrideConst(const ConstMatrixDiag& md)
01396         : start_ptr(const_cast<REAL*>(md.start_ptr)),
01397           end_ptr(const_cast<REAL*>(md.end_ptr)),
01398           stride(md.stride)
01399         { }
01400 
01401                                 // The return statement of of_every below
01402                                 // would implicitly call a constructor
01403 inline ElementWiseStrideConst of_every(const ConstMatrixDiag& md)
01404         { return md; }
01405 
01406 
01407                                 // Access the i-th element of the diagonal
01408 inline const REAL& ConstMatrixDiag::get_ref(const int index) const
01409 { if( index > ncols || index < 1 )
01410       _error("MatrixDiag index %d is out of diag boundaries [%d,%d]",
01411               index,1,ncols);
01412   return start_ptr[stride*(index-1)];
01413 }
01414 
01415 
01416 class MatrixDiag : public ConstMatrixDiag
01417 {
01418   friend class ElementWiseStride;
01419   friend class LAStrideStreamOut;
01420   MatrixDiag(const MatrixDiag&);        // Not implemented and forbidden:
01421   void operator = (const MatrixDiag&);  // no cloning/assignment allowed
01422 public:                                 // Take a col of the matrix
01423   explicit MatrixDiag(Matrix& m) : ConstMatrixDiag(m) {}
01424 
01425   REAL& operator () (const int i)
01426                         { return const_cast<REAL&>(get_ref(i)); }
01427 };
01428 
01429 inline ElementWiseStride::ElementWiseStride(const MatrixDiag& md)
01430         : ElementWiseStrideConst(md) {}
01431 
01432 inline ElementWiseStride to_every(const MatrixDiag& md) { return md; }
01433 
01434 
01435 /*
01436  *------------------------------------------------------------------------
01437  *                      Inline Matrix operations
01438  */
01439 
01440 inline Matrix::Matrix(const int no_rows, const int no_cols)
01441         : DimSpec(no_rows,no_cols)
01442 {
01443   allocate();
01444 }
01445 
01446 inline Matrix::Matrix(const int row_lwb, const int row_upb,
01447                       const int col_lwb, const int col_upb)
01448         : DimSpec(row_lwb,row_upb,col_lwb,col_upb)
01449 {
01450   allocate();
01451 }
01452 
01453 inline Matrix::Matrix(const DimSpec& dimension_specs)
01454         : DimSpec(dimension_specs)
01455 {
01456   allocate();
01457 }
01458 
01459 inline Matrix::Matrix(const LazyMatrix& lazy_constructor)
01460         : DimSpec(lazy_constructor)
01461 {
01462   allocate();
01463   lazy_constructor.fill_in(*this);
01464 }
01465 
01466 
01467                                 // Force a promise of a matrix
01468                                 // That is, apply a lazy_constructor
01469                                 // to the current matrix
01470 inline Matrix& Matrix::operator = (const LazyMatrix& lazy_constructor)
01471 {
01472   is_valid();
01473   if( !(static_cast<const DimSpec&>(lazy_constructor) == *this) )
01474     info(),
01475     _error("The matrix above is incompatible with the assigned "
01476            "Lazy matrix");
01477       
01478   lazy_constructor.fill_in(*this);
01479   return *this;
01480 }
01481                                         // Copy constructor, expensive: use
01482                                         // sparingly
01483 inline Matrix::Matrix(const Matrix& another)
01484         : DimSpec(another)
01485 {
01486   another.is_valid();
01487   allocate();
01488   *this = another;
01489 }
01490 
01491                                 // Resize the matrix to new dimensions
01492 inline void Matrix::resize_to(const int row_lwb, const int row_upb,
01493                  const int col_lwb, const int col_upb)
01494 {
01495   resize_to(DimSpec(row_lwb,row_upb,col_lwb,col_upb));
01496 }
01497 
01498 inline const REAL& MatrixDABase::operator () (const int rown, const int coln) const
01499 {
01500   register int arown = rown-row_lwb;            // Effective indices
01501   register int acoln = coln-col_lwb;
01502 
01503   if( arown >= nrows || arown < 0 )
01504     _error("Row index %d is out of Matrix boundaries [%d,%d]",
01505            rown,row_lwb,nrows+row_lwb-1);
01506   if( acoln >= ncols || acoln < 0 )
01507     _error("Col index %d is out of Matrix boundaries [%d,%d]",
01508            coln,col_lwb,ncols+col_lwb-1);
01509   
01510   return (index[acoln])[arown];
01511 }
01512 
01513 inline Matrix& Matrix::clear(void)      // Clean the Matrix
01514 {
01515   is_valid();
01516   memset(elements,0,nelems*sizeof(REAL));
01517   return *this;
01518 }
01519 
01520 inline void are_compatible(const Matrix& m1, const Matrix& m2)
01521 {
01522   m1.is_valid();
01523   m2.is_valid();
01524   
01525   if( !(static_cast<const DimSpec&>(m1) == m2) )
01526     m1.info(), m2.info(), _error("The matrices above are incompatible");
01527 }
01528 
01529                                 // Assignment
01530 inline Matrix& Matrix::operator = (const Matrix& source)
01531 {
01532   are_compatible(*this,source);
01533   memcpy(elements,source.elements,nelems*sizeof(REAL));
01534   return *this;
01535 }
01536 
01537                                 // Apply a user-defined action to each
01538                                 // element of a collection
01539 inline ElementWiseConstAction&
01540 ElementWiseConst::apply(ElementWiseConstAction& functor) const
01541 {
01542   register const REAL * ep = start_ptr;
01543   while( ep < end_ptr )
01544     functor(*ep++);
01545   return functor;
01546 }
01547 
01548 inline ElementWiseAction&
01549 ElementWise::apply(ElementWiseAction& functor)
01550 {
01551   register REAL * ep = start_ptr;
01552   while( ep < end_ptr )
01553     functor(*ep++);
01554   return functor;
01555 }
01556 
01557 inline ElementWiseConstAction&
01558 ElementWiseStrideConst::apply(ElementWiseConstAction& functor) const
01559 {
01560   for(register const REAL * ep=start_ptr; ep < end_ptr; ep += stride )
01561     functor(*ep);
01562   return functor;
01563 }
01564 
01565 inline ElementWiseAction&
01566 ElementWiseStride::apply(ElementWiseAction& functor)
01567 {
01568   for(register REAL * ep=start_ptr; ep < end_ptr; ep += stride )
01569     functor(*ep);
01570   return functor;
01571 }
01572 
01573                         // Create a blank vector of a given size
01574                         // (indexation starts at 1)
01575 inline Vector::Vector(const int n)
01576         : Matrix(n,1) {}
01577   
01578                         // Create a general lwb:upb vector blank vector
01579 inline Vector::Vector(const int lwb, const int upb) : Matrix(lwb,upb,1,1) {}
01580 
01581                         // Force the promise of a vector
01582 inline Vector& Vector::operator = (const LazyMatrix& source)
01583         { Matrix::operator =(source); return *this; }
01584 
01585                                 // Resize the vector for a specified number
01586                                 // of elements, trying to keep intact as many
01587                                 // elements of the old vector as possible.
01588                                 // If the vector is expanded, the new elements
01589                                 // will be zeroes
01590 inline void Vector::resize_to(const int n)      { Vector::resize_to(1,n); }
01591 
01592 inline void Vector::resize_to(const Vector& v)
01593 {
01594   Vector::resize_to(v.q_lwb(),v.q_upb());
01595 }
01596 
01597                                 // Get access to a vector element
01598 inline REAL& Vector::operator () (const int ind)
01599 {
01600   is_valid();
01601   register const int aind = ind - row_lwb;
01602   if( aind >= nelems || aind < 0 )
01603     _error("Requested element %d is out of Vector boundaries [%d,%d]",
01604            ind,row_lwb,nrows+row_lwb-1);
01605   
01606   return elements[aind];
01607 }
01608 
01609                                 // Make a vector following a recipe
01610 inline Vector::Vector(const LazyMatrix& lazy_constructor)
01611   : Matrix(lazy_constructor)
01612 {
01613   assure(ncols == 1,
01614          "cannot make a vector from a promise of a full matrix");
01615 }
01616 
01617                                 // The following is just a syntactic sugar
01618                                 // to avoid typing of_every and to_every
01619 inline Matrix& Matrix::operator = (const REAL val)
01620         { to_every(*this) = val; return *this; }
01621         
01622 inline Matrix& operator -= (Matrix& m, const double val)
01623         { to_every(m) -= val; return m; }
01624 inline Matrix& operator += (Matrix& m, const double val)
01625         { to_every(m) += val; return m; }
01626 inline Matrix& operator *= (Matrix& m, const double val)
01627         { to_every(m) *= val; return m; }
01628         
01629 inline Matrix::Reference& operator -= (Matrix::Reference& m, const double val)
01630         { to_every(m) -= val; return m; }
01631 inline Matrix::Reference& operator += (Matrix::Reference& m, const double val)
01632         { to_every(m) += val; return m; }
01633 inline Matrix::Reference& operator *= (Matrix::Reference& m, const double val)
01634         { to_every(m) *= val; return m; }
01635 
01636 
01637 inline bool operator ==  (const Matrix& m, const REAL val)
01638         { return of_every(m) == val; }
01639 inline bool operator !=  (const Matrix& m, const REAL val)
01640         { return of_every(m) != val; }
01641 inline bool operator ==  (const Matrix::ConstReference& m, const REAL val)
01642         { return of_every(m) == val; }
01643 inline bool operator !=  (const Matrix::ConstReference& m, const REAL val)
01644         { return of_every(m) != val; }
01645 }
01646 
01647 
01648 #endif
01649 

Generated on Wed Dec 15 21:20:29 2004 for vuVolume by  doxygen 1.3.9.1