readme

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

{
  const double a = 0, const double b = 1.0;
  F1 f1("from the Forsythe book");	// Instantiating the functor
  const double found_location_of_min = fminbr(a,b,f1);
  printf("\nIt took %d iterations\n",f1.get_counter());
}

Note how similar this looks to the previous example. The way fminbr()
is called hasn't changed at all! This is to be expected: from the
operational point of view, a function class is almost identical to a
function pointer. In a sense, a function class is a syntactic wrapper
around a function pointer, which is tucked away into a virtual table
of UnivariateFunctor class. However, there are some differences: For
one, main() does not care any more about global data f1() may use, and
how they should be initialized. It's the job of a functor's
constructor. Also, since the function pointer itself as well as
function's private environment ('counter', for example) are hidden,
there is no way one can mess them up, leave uninitialized, or use
inappropriately.

The validation code vfminbr.cc and vzeroin.cc shows further refinement
of this idea. Please mail me if you have any question or comment.



***** Grand plans
	Computing of a random orthogonal matrix
	Saving/loading of a matrix
	Finding matrix Extrema (and extrema of abs(matrix))
	Compute X'AX for a square matrix
	Compute x'Ax (a square form)
	Asymmetry index
	Add an ArithmeticProgression class ("virtual" constant vector)
	Recompile and beautify FFT and SVD parts, taking advantage of
	streams as much as possible. For example, FFT should take and
	transform a StrideStream rather than a vector. In that case, one can
	apply FFT to the whole matrix, or a separate matrix row or
	a column.
	Make FFT itself a stream class (a Complex stream). This will
	obviate the need for FFT:real(), FFT::abs(), FFT::abs_half() etc.
	methods.
	In FFT, add functions to compute 2D sin/cos transforms.
	Make an FFT FAQ (questions I answered through e-mail,e.g., 2D FFT,
	DFT vs. Fourier integral, etc). Describe in detail how to perform an
	inverse transform.
	Make a special class for SymmetricMatrix
	When gcc starts supporting member function templates, make
	iterator classes for iterating over a vector, two vectors, Matrix
	slices, etc.
	Code to verify a matrix (pseudo)inverse, that is,
	test Moore-Penrose conditions XAX = X, AXA = A; AX and XA are
	symmetric (that is, AX = (AX)', etc.) where X is a (pseudo)inverse
	of A.
	Another SVD application: compute a covariance matrix for a
	given design matrix X, i.e. find the inverse of X'X for a
	rectangular matrix X.
	SVD optimization: when solving Ax=b via SVD, one doesn't need
	to accumulate matrix U: one can apply the transformations directly
	to b. That is, make SVD take matrices U,V as input (which must be
	appropriately initialized), so SVD would accumulate transformations
	in them. U,V don't have to be initialized as unit matrices at all
	(and may be omitted).
	Add ispline(): building a spline and integrating it
	Add slehol: solve a set of SLE with a symmetric matrix
	Lazy matrix operators A+B, multtranspose(A,B), etc.
	Overload operator "," for inline filling in of vectors, as in
		Vector v = VectorInline, 0, 1.0, 3;
	which creates a vector of three elements and assigns initial
	values to them.
	Introduce matrix streams of the second order, that is, a
	stream that traverses a matrix column-by-column, and can create
	a column stream to access the current column. The same for the
	rows.
	Note an article "Scientific Computing: C++ vs Fortran", Nov 1997
	issue of DDJ (described in my DDJ notes)
	Add dump(), validate() methods to every class!
	Note suggestions by Marco Tenuti (in particular, computing
	the rank of a matrix)
	Maybe I have to make a (protected) method
	"bool Matrix::q_within(const REAL * const p) const"
	and use it in LAStreamIn, etc. wherever I need to make sure the pointer
	is within matrix's grid. Maybe this should be Matrix::ConstReference's
	method.
	Consider adding MAC (multiply-accumulate) "instruction" similar
	to the one in DSP. Note that streaming is a typical memory
	access model of DSPs: that's why DSP almost never incorporate a
	a data cache (See "DSP Processors Hit the Mainstream,"
	Jennifer Eyre and Jeff Bier, Computer, Vol. 31, No. 8, August 1998,
	pp. 51-59). Also refer to an article "Smarter Memory: Improving
	Bandwidth for Streamed References", Computer, July 1998, pp.54-63. A
	memory architecture designed in the article achieves low
	overall latencies because it was told by a compiler that
	a stream operation is to follow. Thus streaming is a deep
	software/hardware concept.
	Add streams over Matrix triangles (although they can easily be
	emulated with regular or block streams -- and appropriate 'seek'-ing).
	Construct ElementWiseConst from LAStreamIn or LAStreamOut;
	Make all ElementWise iterators (to_every(), of_every() constructors)
	accept an IRange argument.
	Note mentioning of LinAlg in 
		http://z.ca.sandia.gov/~if-forum/list/msg00000.html
	Replace all 'double' with REAL_EXT
	Think about combinators of a form of_every(of_every1,of_every2);
	In particular, "ElementWiseAction& apply(ElementWiseAction1& functor,
	ConstElementWise&);" where functor takes one argument and
	returns one value. Generalize this 'apply' to take two
	ConstStreams and a functor that takes two arguments. This is
	similar to Scheme's map. Thus it should be possible to
	write "to_every(V).apply(adder(x),of_every(V1));" to
	compute V += V1*x;
	How to calculate: M = N*H*3.5; without modifying the matrices N and H?
	(Boris Holscher's question).
	A member function to return the first element of the matrix, so
	one can easily convert a 1x1 matrix to a scalar
	(Boris Holscher's suggestion)
	Note a similarity between streams and restricted pointers (a
	proposed addition to ANSI C by a Numerical C Extensions Group).
	Restricted pointers are a means of asserting an aliasing constraints
	(if an array is not being modified via a restricted pointer, a
	compiler may assume the array is not being modified at all).


***** Revision history

Version 4.3 (Christmas 1998)

	Improved portability: LinAlg now compiles with gcc 2.8.1, Visual
	C++ and Metrowerk's C++. I heard that SGI's native compiler
	takes LinAlg 4.3 as well.

	The package is complete now: SVD and FFT have been finally
	ported from the 3.2 version.

	Genuine Lambda-abstractions: see vzeroin.cc and vfminbr.cc

	MinMax class (in LinALg.h) is extended with methods to permit
	accumulation of	min and max; see svd.cc for a usage example.

	Introducing IRange and IRangeStride classes in LinALg.h

	Added an explanation of differences between Lazy Matrices and
	expression template styles. See this README file, above.

	Introducing streams over an arbitrary rectangular block of a
	matrix

	All streams have now an ability to 'seek(offset, seek_dir)' in 
	a manner similar to that of ordinary iostreams.
	In particular, one can 'ignore(how_many)' elements in a stream.
	See the explanation above for more details.

	One can _subrange_ a stream: create a new stream that spans over
	a part of another. The new stream may not _extend_ the old one.
	When specifying range boundaries, -INF and INF are permitted, and
	interpreted intelligently. One may create an input stream as a
	subrange of an output stream. See notes above for more details.

	Changed ElementAction and ConstElementAction: they are now
	pure interfaces; the index of the current element is passed
	as an argument to the interface's 'operation'

	Note, SVD::rotate is actually a general-purpose function.

	builtin.h provides now its own implementation of div() on Linux.
	The implementation of div() in linux's libc.a does not agree with
	gcc's optimizations (thanks to Erik Kruus for pointing this out).

	Added pow(double,long) and pow(long,long) to myenv.cc
	They used to be in libg++; yet libstdc++ has dropped them: they
	are not standard yet often convenient

	LAStreams.h is separated from LinAlg.h to make LinAlg.h smaller
	and more manageable

	Trivial changes to the FFT package: taking advantage of
	LAStreams on few occasions, making sure everything works and
	verifies.


Version 4.1

	Completely redesigned the entire hierarchy, based on container/view/
	stream principles. See above for more details.

	Matrix class per se no longer provides a direct access to matrix
	elements. Use MatrixDA, MatrixRow, MatrixColumn, or MatrixDiag
	views.

	MatrixRow, MatrixCol, MatrixDiag, LazyMatrix all inherit from
	DimSpec, that is, all answer queries about their dimensions
	(row/col lower and upper bounds).

	Added lazy matrix operator * (const Matrix&B, const Matrix&C), so
	that it is possible to write now Matrix D = A * B;

	operator *= (Matrix& m, const ConstMatrixDiag& diag)
	is not a member of the Matrix class any more: it does not
	need any privileged access to a matrix object (nor ConstMatrixDiag
	class). The multiplication itself is implemented with streams.
	(and makes a good example of their usage). It is almost just as
	efficient as the priviledged implementation: all loops
	are inlined, no extra procedures/functions are called 
	while traversing the matrices and multiplying their elements.

	By the same token, computing of Vector norms no longer requires
	any direct access to vector's elements. Entire computation is
	a trivial application of an ElementWiseConst iterator.

	No need for separate operations on MatrixColumn, etc:
	it's enough that MatrixColumn can be tacitly converted to
	ElementWise stream, meaning that the whole suite of elementwise
	operations immediately applies.

	Glorified Matrix constructors are gone. They are subsumed by Lazy
	matrices (see the discussion above).

	The matrix inversion algorithm is tinkered with to avoid a direct
	access to matrix' elements (that is, avoid using a column matrix
	index). That's why a method invert() can be applied to a
	Matrix itself, rather than MatrixDA.

	Made constructors (iterators?) for a Vector from MatrixRow/Col/Diag

	ali.cc and hjmin.cc are re-written to take advantage of LAStreams:
	most of the time they use only sequential access to vector
	elements.

	New style casts (const_cast, static_cast) are used throughout.

	Use of mixin inheritance throughout (DimSpec, AREALStream, etc.)

	Tests (validation code) remained mostly unchanged. The only
	modifications were insertions of to_every()/of_every() in a few
	places. This shows how much of the LinAlg interface remained the same.

	ElementPrimAction is replaced with ElementWiseAction/
	ElementWiseConstAction, which are to be applied to a
	(resp, writable and const) ElementWise stream.

	The code and comments in zeroin.cc and fminbr.cc are greatly
	embellished: The code is now more in the C++ true spirit.
	The function under investigation is now a functor
	(a "clause") rather than a mere function pointer;
	DBL_EPSILON is used throughout, instead of my own
	EPSILON (which is deprecated)
	Declarations in math_num.h are adjusted accordingly.

	Validation code vzeroin.cc and vfminbr.cc is almost completely
	re-written to take advantage and show off functors
	(function "clauses") as arguments to zeroin() and fminbr().
	Furthermore, test cases themselves are declared as anonymous,
	"lambda"-like functions, similar to Scheme/Dylan even in
	appearance.

	The validation code is also expanded with new test cases,
	from Matlab:Demos:zerodemo.m

	use M_PI instead of PI (in the fft package)
	FFT package embellished for the new gcc
	determinant.cc is rewritten using only _serial_ access to Matrix
	elements. The algorithm is faster and clearer now.
	Makefile is beautified: c++l is gotten rid of
	Corrected spellings for Householder and Hooke-Jeeves
	corrected an insidious bug in SVD::rip_through() method of QR part
	of svd.cc:
	only abs values of quantities (quantity f) makes sense to compare
	with eps; many thanks to W. Meier <[e-mail]>
	vsvd.cc now includes the test case W. Meier used to show the bug
	fabs() is replaced with abs() throughout the package (which does
	not actually change functionality, but makes the code more generic)
	
Version 3.2
	hjmin(), Hooke-Jeeves optimization, embellished and tested
	ali.cc beautified, using nested functions, etc. (nice style)
	Added SVD, singular value decomposition, and a some code
	to use it (Solving Ax=b, where A doesn't have to be rectangular,
	and b doesn't have to be a vector)
	Minor embellishments
	using bool datatype wherever appropriate
	short replaced by 'int' as a datatype for indices, etc.: all
	modern CPUs handle int more efficiently than short
	(and memory isn't that of a problem any more)
	Forcing promise (LazyMatrix) on assignment
	Testing Matrix(Matrix::Inverted,m) glorified constructor
	added retrieving an element from row/col/diag
	added Matrix C.mult(A,B);	// Product A*B
	*= operation on Matrix slices
	Making a vector from a LazyMatrix
	made sure that if memory is allocated via new, it's disposed
	of via delete; if it was allocated via malloc/calloc, it's
	disposed of via free(). It's important for CodeWarrior, where
	new and malloc() are completely different beasts.
	

Version 3.1
	Deleted dummy destructors (they're better left inferred: it results
	in more optimal code, especially in a class with virtual functions)
	Minor tweaking and embellishments to please gcc 2.6.3
	#included <float.h> and <minmax.h> where missing

Version 3.0 (beta): only Linear Algebra package was changed
	got rid of a g++ extension when returning objects (in a_columns, etc)
	removed "inline MatrixColumn Matrix::a_column(const int coln) const"
	  and likes: less problems with non-g++ compilers, better portability
	Matrix(operation::Transpose) constructors and likes, 
	  (plus like-another constructor) Zero:, Unit, constructors
	invert() matrix
	true copy constructor (*this=another; at the very end)
	fill in the matrix (with env) by doing ElementAction

	cleaned Makefile (pruned dead wood, don't growl creating libla
	  for the first time), a lots of comments as to how to make things
	used M_PI instead of PI
	#ifndef around GNU #pragma's
	REAL is introduced via 'typedef' rather than via '#define': more
	  portable and flexible
	added verify_element_value(), verify_matrix_identity() to the main
	  package (used to be local functions). They are useful in
	  writing the validation code
	added inplace and regular matrix multiplication: computing A*B and A'*B
	all matrix multiplication code is tested now
	implemented A*x = y (inplace (w/resizing))
	Matrix::allocate() made virtual
	improved allocation of vectors (more optimal now)
	added standard function to make an orthogonal Haar matrix
	  (useful for testing/validating purposes)
	Resizing a vector keeps old info now (expansion adds zeros (but still
	  keeps old elements)
	universal matrix creator from a special class: Lazy Matrices

Version 2.0, Mar 1994: Only LinAlg package was changed significantly
	Linear Algebra library was made more portable and compiles now
	  under gcc 2.5.8 without a single warning.
	Added comparisons between a matrix and a scalar (for-every/all -type
	  comparisons)
	More matrix slice operators implemented
	  (operations on a col/row/diag of a matrix, assignments them
	  to/from a vector)
	Other modules weren't changed (at least, significantly), but work
	  fine with the updated libla library

Version 1.1, Mar 1992 (initial revision)