readme

basic linear algebra classes and applications (SVD,interpolation, multivariate optimization)

			Numerical Math Class Library
				version 4.3

***** Version history is at the very end

***** Comments, questions, problem reports, etc.
are all very welcome. Please mail me at
[e-mail] or [e-mail] or [e-mail]

***** Platforms
	I have personally compiled and tested LinAlg 4.3 on
	Linux 2.0.33 i586, gcc 2.7.2.1
	HP 9000/770, HP-UX 10.10, gcc 2.8.1 (with optimization off)
	UltraSparcII, Solaris 2.6, gcc 2.8.1 (with optimization off)
	PowerMac 8500/132, BeOS Release 4, Metrowerk's CodeWarrior C++
	Pentium, WinNT 4.0, Visual C++ 5.0
The previous (4.1) version also works on
	SunSparc20/Solaris 2.4, gcc 2.7.2 and SunWorkshopPro
	SGI, gcc 2.7.2 and SGI's native compiler
	PowerMac 7100/80, 8500/132,
	Metrowerk's CodeWarrior C++, v. 11-12
	DEC Alpha (Digital UNIX Version 5.60), gcc 2.7.2

Note: if g++ version 2.8.1 is used to compile LinAlg, optimization
must be turned off. Otherwise, the compiler crashes with internal
compiler errors. Other compilers as well as gcc v2.7.2 don't seem to have
this problem.

***** Verification files: vmatrix vvector vmatrix1 vmatrix2 vlastreams
			  vali vhjmin vfminbr vzeroin
			  vsvd vslesing

Don't forget to compile and run them, see comments in the Makefile for
details.  The verification code checks to see that all the functions
in this package have compiled and run well. The validation code can also
serve as an illustration of how package's classes and functions
may be employed.


***** Highlights and idioms

1. Data classes and view classes

	Objects of a class Matrix are the only ones that _own_ data, that
is, 2D grids of REALs. All other classes provide different views of/into
the grids: in 2D, in 1D (a view of all matrix elements linearly arranged
in a column major order), a view of a single matrix row, column, or the
diagonal, a view of a rectangular subgrid, etc.

	The views are designed to be as lightweight as possible: most of
the view functionality could be inlined. That is why views do not have any
virtual functions. View classes are supposed to be final, in Java parlance:
It makes little sense to inherit from them. Indeed, views are designed
to provide general and many special ways of accessing matrix elements, as
safely and efficiently as possible. Views turn out to be so flexible that
many former Matrix methods can now be implemented outside of the Matrix class.
These methods can do their job without a privileged access to Matrix'
internals while not sacrificing performance. Examples of such methods:
comparison of two matrices, multiplication of a matrix by a diagonal of
another matrix. 


2. Matrix streams

	Matrix streams provide a sequential view/access to a matrix or
its parts. Many of Linear Algebra algorithms actually require only sequential
access to a matrix or its rows/columns, which is simpler and faster than a
random access. That's why LinAlg stresses streams. There are two kinds of
sequential views: ElementWise streams are transient views that exist only
for the duration of an elementwise action. Hence these views don't require a
"locking" of a matrix. On the other hand, LAStream and the ilk are more
permanent views. An application traverses the stream at its own convenience,
bookmarking certain spots and possibly returning to them later. An LAStream
does assert a shared lock on the matrix, to prevent the grid of matrix
values from moving or disappearing.

	Again, LAStreamIn etc. are meant to denote a stream under user's
control, which the user can create and play with for some time, getting
elements, moving the current position, etc. That is, LAStreamIn and the 
Matrix object from which the stream was created (whose view the stream is)
can coexist. This arrangement is dangerous as one can potentially dispose
of the container before all of its views are closed. To guard against this,
a view asserts a "shared lock" on a matrix grid. Any operation that disposes 
of the grid or potentially moves it in memory (like resize_to()) checks this
lock. Actually, it tries to assert an exclusive lock. Exclusive locking
of a matrix object will always fail if there are any outstanding shared
locks.

	In contrast, ElementWiseConst and ElementWise are meant to be
transient, created only for the duration of a requested operation. That's
why it is syntactically impossible to dispose of the matrix that is
being viewed through a ElementWiseConst or ElementWise object. These
objects do not have any public constructors. This prevents a user from
building the view on his own (and forgetting to delete it). Therefore
creation of ElementWiseConst objects does not require asserting any locks,
which makes this class rather "light".

	Matrix streams may stride a matrix by an arbitrary amount. This
allows traversing of a matrix along the diagonal, by columns, by rows, etc.
Streams can be constructed of a Matrix itself, or from other matrix views
(MatrixColumn, MatrixRow, MatrixDiag). In the latter case, the streams are
confined only to specific portions of the matrix.

	Many methods and functions have been re-written to take advantage
of the streams. Notable examples:

  		// Vector norms are computed via the streams:
  inline double Vector::norm_1(void) const
	{ return of_every(*this).sum_abs(); }
  inline double Vector::norm_2_sqr(void) const
	{ return of_every(*this).sum_squares(); }
  inline double Vector::norm_inf(void) const
	{ return of_every(*this).max_abs(); }

The methods ElementWiseConst::sum_abs(), sum_squares(), and max_abs()
can also be applied to two collections. In this case, they work through
element-by-element differences between the collections.

	An example of using LAStreams:

		// Aitken-Lagrange interpolation (see ali.cc)
  double ALInterp::interpolate()
  {
    LAStreamIn args(arg);	// arg and val are vectors
    LAStreamOut vals(val);
    register double valp = vals.peek();	// The best approximation found so far
    register double diffp = DBL_MAX;	// abs(valp - prev. to valp)

			// Compute the j-th row of the Aitken scheme and
			// place it in the 'val' array
    for(register int j=2; j<=val.q_upb(); j++)
    {
      register double argj = (args.get(), args.peek());
      register REAL&  valj = (vals.get(), vals.peek());
      args.rewind(); vals.rewind();	// rewind the vector streams
      for(register int i=1; i<=j-1; i++)
      {
        double argi = args.get();
        valj = ( vals.get()*(q-argj) - valj*(q-argi) ) / (argi - argj);
      }
      ....
    }
    return valp;
  }
Note, all the streams are light-weight: they use no heap storage and leave no
garbage. All the operations in the snippets above are inlined.

	A stride of a stride stream doesn't have to be an even multiple
of the number of elements in the corresponding collection, as the 
following not-so-trivial example shows. It places a 'value' at the
_anti_-diagonal of a matrix m:	

  m.clear();
  LAStrideStreamOut m_str(m,m.q_nrows()-1);
  m_str.get();	// Skip the first element
  for(register int i=0; i<m.q_nrows(); i++)
    m_str.get() = value;
Again, the matrix is accessed only sequentially. All operations in this snippet
are inlined.

	If an operation involves two ElementWiseStride streams, they don't
have to have the same stride. Note that an ElementWise iterator can always
be converted into an ElementWiseStride stream -- an operation the compiler
does tacitly. This allows an assignment of a vector to a matrix row or
column or a diagonal, or the other way around:
			// Checking the difference between a matrix row and
			// a given vector vr
  assert( of_every(ConstMatrixRow(m,i)).max_abs(of_every(vr)) == 2.0 );

			// multiplying the i-th matrix column by a vector,
			// elementwise
  to_every(MatrixColumn(m,i)) *= of_every(vc);


	Streams can also be constructed implicitly: For example, a function
  bool FPoint::is_step_relatively_small(LAStreamIn hs, const double tau);
in hjmin.cc is declared to take a LAStreamIn as an argument. It is called
however as 
  if( pbase.is_step_relatively_small(h,tau) ) ...
where h is a vector. The vector is converted to a "LAStreamIn hs" when it is
passed as a parameter. Note one can just as easily pass an entire matrix
or a matrix column to is_step_relatively_small(), and the function will take
it. Again, no heap storage is allocated in the process, and all the functions
are inlined. The FFT package uses this technique extensively.

	All streams have now an ability to 'seek(offset, seek_dir)'. In
particular, one can 'ignore(how_many)' elements in a stream. how_many and
the offset can be positive as well as zero or negative quantities. 'seek_dir'
specifies that the offset will count from the beginning of a stream,
end of the stream, or from the stream's current position. This is very
similar to seeking in an ordinary iostream. Trying to seek past the upper
end of a stream does _not_ generate an error: it merely sets the EOF
condition on the stream. This condition can be checked later by an eof()
method. This lets one check to see if a given offset is within a stream,
without incurring a stream exception. This shows an easy way to implement
'has(offset);' and 'peek(offset);' functions suggested by Erik Kruus.
Note, trying to get() or peek() from a stream at the EOF results in an
immediate error. See vlastreams.cc for more details and examples of many
a seek.

	As yet another example of clarity and efficiency of streams,
consider the following method of a SVD package. In LinAlg 3.2, it is
used to read:	
				// Carry out U * S with a rotation matrix
				// S (which combines i-th and j-th columns
				// of U, i>j)
inline void SVD::rotate(Matrix& U, const int i, const int j,
		   const double cos_ph, const double sin_ph)
{
  MatrixColumn Ui(U,i), Uj(U,j);
  for(register int l=1; l<=U.q_row_upb(); l++)
  {
    REAL& Uil = Ui(l); REAL& Ujl = Uj(l);
    const REAL Ujl_was = Ujl;
    Ujl =  cos_ph*Ujl_was + sin_ph*Uil;
    Uil = -sin_ph*Ujl_was + cos_ph*Uil;
  }
}

The new version of the package, v4.3, rewrites this code in the
following way, avoiding random access to Matrix elements (and incident
range checks for MatrixColumns' indices).

inline void SVD::rotate(Matrix& U, const int i, const int j,
		   const double cos_ph, const double sin_ph)
{
  LAStreamOut Ui(MatrixColumn (U,i));
  LAStreamOut Uj(MatrixColumn (U,j));
  while( !Ui.eof() )
  {
    REAL& Uil = Ui.get(); REAL& Ujl = Uj.get();
    const REAL Ujl_was = Ujl;
    Ujl =  cos_ph*Ujl_was + sin_ph*Uil;
    Uil = -sin_ph*Ujl_was + cos_ph*Uil;
  }
}



LinAlg's implementation and validation code provides many more examples
of flexibility, clarity, and efficiency of streams.


4. Subranging a stream

	LinAlg 4.3 implements subranging of a stream: creation of a new
stream that spans over a part of another. The latter can be either
a stride 1 or an arbitrary stride stream. The new stream may not _extend_
the old one: subranging is therefore a safe operation. A child stream
is always confined within its father's boundaries. A substream, once
created, lives independently of its parent: either of them can go out of
scope without affecting the other.

	To create a substream, one has to specify its range as a 
[lwb,upb] pair: an all-inclusive range. The 'lwb' may also be
given as -IRange::INF, and 'upb' as IRange::INF. These special values
are interpreted intelligently. Given an output stream, one may create
either an immutable, input substream, or allow modifications to a part
of the original stream. Obviously, given a read-only stream, only
read-only substreams may be created: the compiler will see to that at
_compile_ time. File sample_ult.cc provides a good example of using
substreams to efficiently reflect the upper triangle of a square matrix
onto the lower one (yielding a symmetric matrix).  See vlastreams.cc and
SVD for more extensive examples of subranging.

	As a good illustration to the power and efficiency of streams,
consider the following snippet from SVD::left_householder() (file svd.cc).
In LinAlg 3.2, the snippet was programmed as:

   for(j=i+1; j<=N; j++)		// Transform i+1:N columns of A
   {
     MatrixColumn Aj(A,j);
     REAL factor = 0;
     for(k=i; k<=M; k++)
       factor += UPi(k) * Aj(k);	// Compute UPi' * A[,j]
     factor /= beta;
     for(k=i; k<=M; k++)
       Aj(k) -= UPi(k) * factor;
   }

Version 4.3 uses substreams (which span only a part of matrix'
columns):

   IRange range = IRange::from(i - A.q_col_lwb());
   LAStreamOut UPi_str(MatrixColumn(A,i),range);
   ...
   for(register int j=i+1; j<=N; j++)	// Transform i+1:N columns of A
   {
     LAStreamOut Aj_str(MatrixColumn(A,j),range);
     REAL factor = 0;
     while( !UPi_str.eof() )
       factor += UPi_str.get() * Aj_str.get();	// Compute UPi' * A[,j]
     factor /= beta;
     for(UPi_str.rewind(), Aj_str.rewind(); !UPi_str.eof(); )
       Aj_str.get() -= UPi_str.get() * factor;
     UPi_str.rewind();
   }
 

5. Streams over an arbitrary rectangular block of a matrix

	LinAlg 4.3 permits LAstreams that span an arbitrary  rectangular
block of a matrix (including the whole matrix, a single matrix element,
a matrix row or a column, or a part thereof). You simply specify a matrix, 
and the range of rows and columns to include in the stream.

vlastreams.cc verification code shows very many examples. The following
is a annotated snippet from vlastreams' output (vlastreams.dat)

	---> Test LABlockStreams for 2:12x0:20
        checking accessing of the whole matrix as a block stream...
        checking modifying the whole matrix as a block stream...
					# The first column of m
        checking LABlockStreams clipped as [-INF,INF] x [-INF,0]
					# The second column of m
        checking LABlockStreams clipped as [2,INF] x [1,1]
					# A part of the last column of m
        checking LABlockStreams clipped as [3,12] x [20,INF]
					# The first row of m
        checking LABlockStreams clipped as [-INF,2] x [0,INF]
					# The second row of m
        checking LABlockStreams clipped as [3,3] x [-INF,20]
					# A part of the last row of m
        checking LABlockStreams clipped as [12,INF] x [-INF,19]
					# The first matrix element
        checking LABlockStreams clipped as [2,2] x [0,0]
					# The last matrix element
        checking LABlockStreams clipped as [12,INF] x [20,INF]
					# A 2x3 block at the upper-right corner
        checking LABlockStreams clipped as [3,4] x [18,INF]
					# A 3x2 block at the lower-left corner
        checking LABlockStreams clipped as [9,11] x [2,3]
        checking modifying LABlockStreams clipped as [4,4] x [3,3]

With block streams, it is trivial to assign a block of one matrix
to a block of another matrix. See vlastreams.cc for examples.


6. Direct access to matrix elements

	In LinAlg 3.2 and earlier, a Matrix contained a column index.
The index was used to speed up direct access to matrix elements 
(like in m(i,j)). The index was allocated on heap and initialized upon
Matrix construction. In LinAlg 4+, a Matrix object no longer allocates
any index. The Matrix class therefore does not allow direct access to
matrix elements. This is done on purpose, to slim down matrices and make