------------------------->  GNU Sather - sourcefile  <-------------------------
-- Copyright (C) 199x by International Computer Science Institute            --
-- This file is part of the GNU Sather library. It is free software; you may --
-- redistribute  and/or modify it under the terms of the GNU Library General --
-- Public  License (LGPL)  as published  by the  Free  Software  Foundation; --
-- either version 2 of the license, or (at your option) any later version.   --
-- This  library  is distributed  in the  hope that it will  be  useful, but --
-- WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY --
-- or FITNESS FOR A PARTICULAR PURPOSE. See Doc/LGPL for more details.       --
-- The license text is also available from:  Free Software Foundation, Inc., --
-- 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA                     --
-------------->  Please email comments to <bug-sather@gnu.org>  <--------------

-- This class had calls to a c version that seems to no longer exist!
-- I commented that stuff out and added the test to test-oblig.sa
-- (gomes Nov21/95)
-- *************************************************************
-- Converted to use the new matrix and vector classes July 96
-- The test class for SVD seems to generate results that are slightly
-- off from the expected results in 0.2. 
-- Bo Hagstedt<bo.hagstedt@mailbox.swipnet.se> looked at the code
-- and found that it was good (and fixed several errors).
-- The IF statments protect accesses outside the matrices when
-- the loop index i points to the last row column as the sub loop will never 
-- start when the start index is greater han the stop indexand a sub loop uses
-- i+1 as its index.  
-- Thanks also to Alex Cozzi who looked at this code earlier.
-- gomes@icsi.berkeley.edu
-------------------------------------------------------------------
class NR_SVD is
   -- Singular value decomposition derived from Numerical Recipes.
   -- This algorithm is known to have a problem with an obscure case (but
   -- works almost all of the time). Better algorithms are being worked on.
   
   pythag(a,b:FLT):FLT is
      -- Square root of a^2+b^2 without destructive overflow or underflow.
      res:FLT;
      at:FLT:=a.abs;  bt:FLT:=b.abs; ct:FLT;
      if at>bt then ct:=bt/at; res:=at*(1.0+ct*ct).sqrt;
      elsif bt/=0.0 then ct:=at/bt; res:=bt*(1.0+ct*ct).sqrt;
      end; -- if
      return res;
   end; -- pythag
   
   svd(a:MAT, w:VEC, v:MAT)
      pre a.nr>=a.nc and w.dim=a.nc and v.nr=a.nc and v.nc=a.nc	is
	 -- Computes the singular value decomposition of a. When finished,
	 -- a w v^T equals the old a. Must have a.nr>=a.nc.
	 -- Based on Numerical Recipes in C, p. 68.
      m:INT:=a.nr; n:INT:=a.nc;
      -- Householder reduction to bidiagonal form (should be separate).
      rv1:VEC:=rv1.create(n); 
      g,scale,anorm,c,f,h,s,x,y,z:FLT; 
      i,its,j,jj,k,l,nm:INT;  flag:BOOL;
      -- Householder reduction to bidiagonal form
      i:=0; loop until!(i>=n);
	 l:=i+1;  rv1[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
	 if i<m then 
	    k:=i; loop until!(k>=m); scale:=scale+a[k,i].abs; k:=k+1 end; 
	    if scale/=0.0 then
	       k:=i; loop until!(k>=m); 
		  a[k,i]:=a[k,i]/scale; s:=s+a[k,i]*a[k,i]; k:=k+1 
	       end; -- loop
	       f:=a[i,i]; 
	       if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
	       h:=f*g-s; a[i,i]:=f-g;
	       j:=l; loop until!(j>=n);
		  s:=0.0; 
		  k:=i; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
		  f:=s/h;
		  k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i];k:=k+1 end;
		  j:=j+1
	       end; -- loop
	       k:=i; loop until!(k>=m); a[k,i]:=a[k,i]*scale; k:=k+1 end;
	    end; -- if
	 end; -- if
	 w[i]:=scale*g; g:=0.0; s:=0.0; scale:=0.0;
	 if i<m and i/=n-1 then
	    k:=l; loop until!(k>=n); scale:=scale+a[i,k].abs; k:=k+1 end;
	    if scale/=0.0 then
	       k:=l; loop until!(k>=n);
		  a[i,k]:=a[i,k]/scale; s:=s+a[i,k]*a[i,k]; k:=k+1
	       end; -- loop
	       f:=a[i,l];
	       if f<0.0 then g:=s.sqrt else g:=-s.sqrt end;
	       h:=f*g-s; a[i,l]:=f-g;
	       k:=l; loop until!(k>=n); 
		  rv1[k]:=a[i,k]/h;
		  k:=k+1 
	       end;
	       j:=l; loop until!(j>=m); 
		  s:=0.0;
		  k:=l; loop until!(k>=n); s:=s+a[j,k]*a[i,k]; k:=k+1 end;
		  k:=l; loop until!(k>=n); a[j,k]:=a[j,k]+s*rv1[k];k:=k+1 end;
		  j:=j+1
	       end; -- loop
	       k:=l; loop until!(k>=n); a[i,k]:=a[i,k]*scale; k:=k+1 end;
	    end; -- if
	 end; -- if
	 anorm:=anorm.max(w[i].abs+rv1[i].abs);
	 i:=i+1;
      end; -- loop
      -- Accumulation of right-hand transformations.
      i:=n-1; loop until!(i<0); 
	 if i<n-1 then
	    if g/=0.0 then
	       j:=l; loop until!(j>=n); 
		  -- double division to avoid possible underflow
		  v[j,i]:=(a[i,j]/a[i,l])/g; 
		  j:=j+1 
	       end;
	       j:=l; loop until!(j>=n);
		  s:=0.0;
		  k:=l; loop until!(k>=n); s:=s+a[i,k]*v[k,j]; k:=k+1 end;
		  k:=l; loop until!(k>=n); v[k,j]:=v[k,j]+s*v[k,i]; k:=k+1 end;
		  j:=j+1
	       end; -- loop
	    end; -- if
	    j:=l; loop until!(j>=n); v[i,j]:=0.0; v[j,i]:=0.0; j:=j+1 end;
	 end; -- if
	 v[i,i]:=1.0; g:=rv1[i]; l:=i;
	 i:=i-1;
      end; -- loop
      -- Accumulation of left-hand transformations.
      -- the new code is different : i=IMIN(m,n)
      --      i:=n-1; loop until!(i<0); 
      i:=n.min(m)-1; loop until!(i<0); 
	 l:=i+1; g:=w[i];
	 -- Redundant if removed by BH
	 j:=l; loop until!(j>=n); a[i,j]:=0.0; j:=j+1 end;
	 if g/=0.0 then
	    g:=1.0/g;
	    
	    -- if i/=n-1 then BH
	    j:=l; loop until!(j>=n);
	       s:=0.0;
	       k:=l; loop until!(k>=m); s:=s+a[k,i]*a[k,j]; k:=k+1 end;
	       f:=(s/a[i,i])*g;
	       k:=i; loop until!(k>=m); a[k,j]:=a[k,j]+f*a[k,i]; k:=k+1 end;
	       j:=j+1
	    end; -- loop
	    -- end; -- if
	    j:=i; loop until!(j>=m); a[j,i]:=a[j,i]*g; j:=j+1 end;
	 else
	    j:=i; loop until!(j>=m); a[j,i]:=0.0; j:=j+1 end;
	 end; -- if
	 a[i,i]:=a[i,i]+1.0;
	 i:=i-1;
      end; -- loop
      -- Diagonalization of the bidiagonal form.
      k:=n-1; loop until!(k<0);	-- loop over singular values
	 --its:=0; loop until!(its>=30); -- loop over allowed iterations
	 -- Should read as follows in order to get the "its" error message. /BH
	 its:=1; loop until!(its>30); -- loop over allowed iterations
	    flag:=true; 
	    l:=k; loop until!(l<0); -- test for splitting
	       nm:=l-1;		-- rv[1] is always zero
	       if rv1[l].abs.flt+anorm=anorm then flag:=false; break! end;
	       if w[nm].abs.flt+anorm=anorm then break! end;
	       l:=l-1
	    end; -- loop
	    if flag then
	       c:=0.0; s:=1.0; 
	       i:=l; loop until!(i>k);
		  f:=s*rv1[i];
		  -- Corrected the typo, rv1[1] => rv1[i]. /BH
		  rv1[i]:=c*rv1[i];
		  if f.abs.flt+anorm/=anorm then
		     g:=w[i]; h:=pythag(f,g); w[i]:=h; h:=1.0/h;
		     c:=g*h; s:=-f*h;
		     j:=0; loop until!(j>=m); 
			y:=a[j,nm]; z:=a[j,i]; 
			a[j,nm]:=y*c+z*s; a[j,i]:=z*c-y*s;
			j:=j+1
		     end; -- loop
		  end; -- if
		  i:=i+1
	       end; -- loop
	    end; -- if
	    z:=w[k];
	    if l=k then		-- convergence
	       if z<0.0 then	-- singular value is made non-negative
		  w[k]:=-z;
		  j:=0; loop until!(j>=n); v[j,k]:=-v[j,k]; j:=j+1 end;
	       end; -- if
	       break!;
	    end; -- if
	    if its=30 then #ERR+"SVD: 30 iterations, no convergence\n"; end;
	    -- shift from bottom 2 by 2 minor
	    x:=w[l]; nm:=k-1; y:=w[nm]; g:=rv1[nm]; h:=rv1[k];
	    f:=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g:=pythag(f,1.0);
	    if f>=0.0 then f:=((x-z)*(x+z)+h*((y/(f+g.abs))-h))/x
	    else f:=((x-z)*(x+z)+h*((y/(f-g.abs))-h))/x end;
	    -- Next QR transformation
	    c:=1.0; s:=1.0; 
	    j:=l; loop until!(j>=nm+1);
	       i:=j+1; g:=rv1[i]; y:=w[i]; h:=s*g; g:=c*g; z:=pythag(f,h);
	       rv1[j]:=z; c:=f/z; s:=h/z; f:=x*c+g*s; g:=g*c-x*s;
	       h:=y*s; y:=y*c; 
	       jj:=0; loop until!(jj>=n);
		  x:=v[jj,j]; z:=v[jj,i]; v[jj,j]:=x*c+z*s; v[jj,i]:=z*c-x*s;
		  jj:=jj+1
	       end; -- loop
	       z:=pythag(f,h); w[j]:=z;
	       if z/=0.0 then z:=1.0/z; c:=f*z; s:=h*z end;
	       f:=(c*g)+(s*y); x:=(c*y)-(s*g);
	       jj:=0; loop until!(jj>=m); 
		  y:=a[jj,j]; z:=a[jj,i]; a[jj,j]:=y*c+z*s; a[jj,i]:=z*c-y*s;
		  jj:=jj+1;
	       end; -- loop
	       j:=j+1
	    end; -- loop
	    rv1[l]:=0.0; rv1[k]:=f; w[k]:=x;
	    its:=its+1;
	 end; -- loop
	 k:=k-1
      end; -- loop
   end; -- svd

   
end; -- class NR_SVD

-------------------------------------------------------------------

class SVD_TEST is
   include TEST;

   main is
      class_name("SVD");
      a:MAT:=MAT::create(3,4);
      a.inplace_uniform_random;
      svd_matrix_test(1, a);
      a:=MAT::create(3,4);
      a.inplace_uniform_random;
      svd_matrix_test(2, a);
      a:=MAT::create(4,3); a.inplace_uniform_random;
      svd_matrix_test(3, a);
      a:=MAT::create(5, 5); a.inplace_uniform_random;
      svd_matrix_test(4, a);
      a:=MAT::create(||0.0,1.0,0.0|,|0.0,1.0,1.0|,|0.0,0.0,0.0||);
      svd_matrix_test(5, a);	-- the bad example for Numerical Recipes

      svd_back_sub_test(1,
	    MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
	    VEC::create(|1.0,0.0,0.0|));
      svd_back_sub_test(2,
	    MAT::create(||1.0,0.0,0.0|,|0.0,2.0,0.0|,|0.0,0.0,3.0||),
	    VEC::create(|1.0,1.0,1.0|));
      svd_back_sub_test(3,
	    MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
	    VEC::create(|1.0,1.0,1.0|));
      svd_back_sub_test(4,
	    MAT::create(||1.0,2.0,3.0|,|2.0,2.0,3.0|,|3.0,2.0,3.0||),
	    VEC::create(|1.0,2.0,3.0|));
      svd_back_sub_test(5,
	    MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
	    ,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
	    VEC::create(|1.0,1.0,1.0,1.0,1.0|));
      svd_back_sub_test(6,
	    MAT::create(||1.0,2.0,3.0,4.0,5.0|,|2.0,3.0,4.0,5.0,1.0|
	    ,|3.0,4.0,5.0,1.0,2.0|,|4.0,5.0,1.0,2.0,3.0|,|5.0,1.0,2.0,3.0,4.0||),
	    VEC::create(|1.0,2.0,3.0,4.0,5.0|));
      svd_back_sub_test(7,
	    MAT::create(||1.4, 2.1, 2.1, 7.4, 9.6|
	    , |1.6, 1.5, 1.1, 0.7, 5.0|, |3.8, 8.0, 9.6, 5.4, 8.8|
	    , |4.6, 8.2, 8.4, 0.4, 8.0|, |2.6, 2.9, 0.1, 9.6, 7.7||),
	    VEC::create(|1.1, 1.6, 4.7, 9.1, 0.1|));

      lfm:MAT:=MAT::create(||1.0,2.0|,|3.0,4.0|,|5.0,6.0||);
      fl_lfi:FLIST{VEC}:=FLIST{VEC}::create;
      fl_lfi:=fl_lfi.push(VEC::create(|1.0,0.0|));
      fl_lfi:=fl_lfi.push(VEC::create(|1.0,1.0|));
      fl_lfi:=fl_lfi.push(VEC::create(|0.0,1.0|));
      fl_lfi:=fl_lfi.push(VEC::create(|2.0,1.0|));
      lfi:ARRAY{VEC} := fl_lfi.array;
      
      fl_lfo:FLIST{VEC}:=FLIST{VEC}::create; 
      loop
	 fl_lfo:=fl_lfo.push(lfm.times_vec(lfi.elt!));
      end;
      lfo: ARRAY{VEC} := fl_lfo.array;
      loop
	 #OUT +"lfo: "+lfo.elt!.str+"\n";
      end;
      
--d      test("to_linear_fit_of",
--d	    lfm.copy.inplace_zero.to_linear_fit_of(lfi,lfo).str, lfm.str);
      fl_lfoa:FLIST{VEC}:=FLIST{VEC}::create;
      loop
	 fl_lfoa:=fl_lfoa.push(lfo.elt!+(VEC::create(|1.0,2.0,3.0|)));
      end;
      lfoa:ARRAY{VEC} := fl_lfoa.array;
      
      lfma:MAT:=MAT::create
	    (||1.0,2.0,1.0|,|3.0,4.0,2.0|,|5.0,6.0,3.0||);
--d      test("to_affine_fit_of",
--d	    lfma.copy.inplace_zero.inplace_affine_fit_of(lfi,lfoa).str, lfma.str);
--d      wt:FLIST{FLT}:=FLIST{FLT}::create;
--d      wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0); wt:=wt.push(1.0);
--d      lfmcopy ::= lfm.copy;
--d      lfmcopy.inplace_zero;
--d      lfmcopy.inplace_weighted_linear_fit_of(lfi,lfo,wt);
--d      test("to_weighted_linear_fit_of",lfmcopy.str,lfm.str);
--d      test("to_weighted_affine_fit_of",
--d	    lfma.copy.inplace_zero.inplace_weighted_affine_fit_of(lfi,lfoa,wt).str
--d	    , lfma.str);
      finish;
   end;

   svd_matrix_test(n:INT, a:MAT) is
      -- Test the singular value decomposition of `a' labelling the test
      -- by `n'.
      u:MAT:=u.create(a.nr.max(a.nc),a.nc);
      v:MAT:=v.create(a.nc, a.nc); w:VEC:=w.create(a.nc);
      a.svd_in(u,w,v);
      u.inplace_times_diagonal(w);
      v.inplace_trans;
      tmp:MAT:=u.copy;
      tmp.inplace_arg_times_arg(u,v);
      acpy::= a.copy;
      acpy.inplace_zero;
      acpy.inplace_portion_of_arg(tmp);
      cmp:MAT:=acpy;
      unchecked_test("svd_in "+n, cmp.str, a.str);
   end; -- svd_test_matrix

   svd_back_sub_test(n:INT, a:MAT, b:VEC) is
      -- Test `svd_back_sub' on the problem `a.x=b', label with `n'.
      size:INT:=a.nr;
      u:MAT:=u.create(size,size); v:MAT:=v.create(size,size);
      w:VEC:=w.create(size); a.svd_in(u,w,v);
      wmax:FLT:=w.max_value; wmin:FLT:=wmax*(0.000001);
      k:INT:=0; loop until!(k=size);
	 if w[k]<=wmin then w[k]:=0.0 end; k:=k+1 
      end; -- loop
      x:VEC:=VEC::create(size);
      a.svd_back_sub(u,w,v,b,x);
      c:VEC:=a.times_vec(x);
      unchecked_test("svd_back_sub "+n, c.str, b.str);
   end; -- svd_back_sub_test



end; -- class TEST_SVD

-------------------------------------------------------------------
