Andrew U. Frank
it appears possible to define matrix operations like LU-decomposition, SVD, inverse etc. in a recursive fashion (see transpose in the prelude).
is this possible? has anybody code in haskell which does this (not asking for high performance here, more instructive value).
I would like to see this, too. Myself I did some experiments in trying to implement the transformation to upper diagonal form by rational Givens transforms. The following code fragment does this recursively: -- rationalGivens :: Fractional a => a -> [a] -> [[a]] -> [[a]] -- marginal cases first rationalGivens qq [x] ((pp:[]):[]) = ( pp + qq * x * x : [] ) : [] rationalGivens _ [] [[]] = [[]] rationalGivens qq xy pmat | qq == 0 = pmat rationalGivens qq (x:y) (ppr:pmat) | x == 0 = ppr:rationalGivens qq y pmat rationalGivens qq (x:y) [[]] = ( qq * x * x : (1/x).*y ) : [] -- main case rationalGivens qq (x:y) ((pp:r):pmat) = let pp' = pp + qq * x * x qq' = pp * qq / pp' s = x * qq / pp' y' = y - x .* r r' = r + s .* y' in ( pp' : r' ) : rationalGivens qq' y' pmat -- rationalGivens 1 [2,1] [[1,0],[1]] == [[5.0,0.4],[1.2]] Arrays are double lists in this code, q a scale factor (starting with 1) xy a row vector to be added to the u.t. matrix pmat. The diagonal elements of pmat contain the square of the scale factor of the row, i.E. the matrix [[5.0,0.4],[1.2]] is to be interpreted as the product (sqrt(5) 0 ) ( 1 0.4 ) ( 0 sqrt(1.2)) ( 0 1 ) -- Dipl.-Math. Wilhelm Bernhard Kloke Institut fuer Arbeitsphysiologie an der Universitaet Dortmund Ardeystrasse 67, D-44139 Dortmund, Tel. 0231-1084-257