[clean-list] The Clas library

[Siegfried Gonzi]

Ladies and Gentlemen!


Attention please; Austria proudly presents:


The Clas-library from Thorsten Zoerner is a pretty first step when
you have to deal with linear algebra problems and one does not want
to use packages like Matlab or Yorick and so on.

Surely, they provide many, many built-in functions and compute different
defined problems very quickly. But I say it for the 1000 time: It is
impossible with this programs to write readable code. That said.


The Clean Clas-library  function "solveLSPartPiv" let you handle the
following
problem structure:

A x = y


where A is an array; x is the vector with solutions (which
"solveLSPartPiv"
will deliver you); and y is a know vector.

Clearly speaking: you know A and y and get x.

But often you want not only solve for one vector but for a set of
vectors. You
can then solve and call "solveLSPartPiv" for every vector again. But
this is
not very efficient, because you have to calculate the LU-factorization
again for
every same array, but one calculation for a set of vectors would be
enough.

Therefore:

A X = Y

can easily be solved. Here A is the array as before; X is the set of
vectors which
holds as colums the one special vector; and Y holds as columns the one
special known
vectors.


The following small program is a modified version of Zoerners'
"solveLSPartPiv"
function. Be aware: as opposed to the above definition Y expects the
know vectors
as rows (you may have to transpose it if needed!). And X deliver the
solution
vectors also in rows:

Xrow1: solution vector 1
Xrow2: solution vector 2
and so on


The reason is only due to efficiency exhaustion (because in Clean arrays
are stored
row after row).

==
module lin
import StdEnv,Clas3, MersenneTwister


// This program is freeware
// I do not take any responsibility which this program
// may cause to you
// Developed an a Macintosh with Clean 1.3.3
// Tested on a Macintosh and Sun workstation (Clean 1.3.3)
// IF YOU USE it on a Sun, please change all occurences of
// "#" in the type definition of both the Clas-library itself
// and my program to "!"
// The original Clas-library is copyright by Thorsten Zoerner

solveLSPartPivM ::  {#*{#Real}} {#*{#Real}} -> {#{#Real}}
solveLSPartPivM   a x               
	#! (lu, p) = luPartPiv a
	= solveLSPartPiv_i 0 lu p x
where
	solveLSPartPiv_i :: !Int *{#{#Real}} *{#Int} *{#{#Real}} -> *{#{#Real}}
	solveLSPartPiv_i  n lu p x=:{ [n] = x_row }
		| (n == size x) = x
		#! pv = perVec p x_row
		#! y = forwardSubstUnit lu  pv 
		#! yb = backwardSubst lu y 
		= solveLSPartPiv_i (n+1) lu p {x & [n] = yb} 
 
// Generation of a random number array	
// The MersenneTwister package for generating
// random numbers is available from the Clean homepage
makeMat:: !Int !Int !Int {#*{#Real}} -> {#*{#Real}}
makeMat n i j mat
	| n == (size mat) = mat
	#! el = taken j (genRandReal (4397+n) ) 
	#! ell = { (elm+1.0) \\ elm <- el }
	 = makeMat (n+1) i j { mat & [n] = ell }
where
	taken:: !Int [!Real] -> [!Real]
	taken n [] = []
	taken n [x:r] 
		| n == 0 = []
 		= [x:taken (n-1) r]

makeMatt:: !Int !Int -> {#*{#Real}}
makeMatt i j = { { 0.0 \\ j <- [0..j-1] } \\ j <- [0..i-1] }



TakeOne:: {#{#Real}} -> !Real
TakeOne a = a.[0,0]

Start:: !Real
Start = TakeOne( solveLSPartPivM  a b ) 		
where
	a :: {#*{#Real}}
	//a = {{0.0,23.45,2.34},{3.45,0.0,0.0},{9.23,3.45,12.34}}
	a = makeMat 0 300 300 (makeMatt 300 300)
	b :: {#*{#Real}}
	//b = {{2.3,4.3,5.6},{2.56,2.45,6.78},{9.78,54.6,3.4}}
	b = makeMat 0 300 300 (makeMatt 300 300)
	//b = {{1.0,2.0,34.0}}
==


You can test with an array of random numbers. This is done with:

makeMatt:: !Int !Int !Int {#*{#Real}} -> {#*{#Real}}

you pass the counter index (starting at 0), i (rows), j (columns) and 
the scaffolding-array and you get an array with random numbers. 
In order to use the random number generator you have to download 
the so called "MerseneTwister"-package. You get it from the Clean
homepage.


Timings:


Macintosh:

The above method works and delivers for test cases the same solution as
compared to Yorick. But there remains one fact which bothers me:

Yorick is about 2 times faster!

It is quite interesting that Yorick and Clean takes for the
LU-factorization
with a great array (300x300) nearly the same time (about 10sec). But
then
Yorick takes only about 5sec for the forward- and backward substitution
as
opposed to Clean with 20sec. This all on a 100MHz Macintosh (603e
processor).

This is remarkable, because it seems that with Clean the limiting factor
is
memory-access and Yorick does not show that problem. This is even more
of interest,
because as I can remember Yorick also features a limiting-memory problem
when it comes to the FFT and matrix-matrix-multiplication (then Yorick
and
Clean show the same timings).

I think that Yorick uses some tuned or quite efficient forward- and
backward
substitution algorithm. I even cannot check it, because the Yorick
function
"LUsolve" is based on an external lin-algebra-compiled routine (maybe in
C?).

I tried also the version without pivoting and even then Clean shows the
same 
timing. This was a reason for me to think on memory as the limiting
factor.
Why Yorick cope the problem that cute is out of my vision.

But nobody should now draw conclusions from this observation before one
know
the exact Yorick-algorithm:


Sun OS:

Here the situation is as follows: The random 300x300 array takes on the
(old) Sun
16sec and Yorick takes 10sec. Also noteable is the fact that Clean can
handle
strict arrays on Unix only. Maybe unboxed arrays will give another 20%
and then Yorick
and Clean would show about the same execution time.


And now to another strangeness:

Try this (subfunction in solveLSPartPivM):

solveLSPartPiv_i :: !Int *{#{#Real}} *{#Int} *{#{#Real}} -> *{#{#Real}}

with:

solveLSPartPiv_i :: !Int *{#{#Real}} *{#Int} {#*{#Real}} -> {#*{#Real}}
	

It only emerged after many hours of trial and error that:

{#*{#Real}}

does not work and one has to supersede it instead with:

*{#{#Real}}

otherwise one will always get the error message that it cannot be
coarced:

I have stopped smoking since 5 years now. A doctor will say this is good
and
will save me 7 years more on life. But the doctor does not know Clean
and
also does not know that I sometimes program in Clean. One can then sum
up
and show that 7 years more on life will equal out, because Clean array
programming
will give you also 7 years but in the negative sense.


Anybody installed Matlab on his computer? You could put-in in a good
duty when you
post the execution time of the solution of the above 300x300 random
matrix 
compared to Clean on your machine.

As I have said at other places: the Macintosh series 5300 (which I also
posses) 
have been a crap (a rather weak transition from the 68k to the PPC
architecture).
This means on other machines the situation can be very different (as has
been
verified by the SunOS). Newer Macintoshes will perform better in this
respect.


NOTE: If you want to use the above program on a Sun, you must change all
occurences of "#" to "!" and you must also change all occurences in the
Clas libraries to "!" (please, only in the type definitions!). 
Because there is a bug with unboxed arrays on the Sun.
But it works with strict arrays.

NOTE ALSO:

You can also use the above program to find the inverse of a matrix which
has got zeros in the diagonal. The algorithm for finding the inverse 
in the original Clas-library does not include pivoting (pivoting is only
used with LU-factorization).

But you can circumvent this with (see my above program): 

solveLSPartPivM ::  {#*{#Real}} {#*{#Real}} -> {#{#Real}}

e.g:

solveLSPartPivM A Y

where A is the array from which you want the inverse and Y is a 
diagonal-1-array (all elements are zero except the diagonals which are 
equal 1.0).


Thank you fan for your attention and autographs will follow tomorrow,
S. Gonzi