Artifact ef5997c8815407d6b3f6be669e439cabe5c0c143a05ba6090bcc894605f2a2d4:
- Executable file
r37/packages/sparse/spsvd.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 17590) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/sparse/spsvd.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 17590) [annotate] [blame] [check-ins using]
%**********************************************************************% % % % Computation of the Singular Value Decomposition of sparse matrices % % containing numeric entries. Uses specific rounded number routines to % % speed things up. % % % % Author: Stephen Scowcroft. Date: June 1995 % % (based on code by Matt Rebbeck) % % % % The algorithm was taken from "Linear Algebra" - J.H.Wilkinson % % & C. Reinsch % % % %**********************************************************************% module spsvd; symbolic procedure spsvd(A); % % Computation of the singular values and complete orthogonal % decomposition of a real rectangular matrix A. % % A = tp(U) diag(q) V, U tp(U) = V tp(V) = I, % % and q contains the singular values along the diagonal. % (tp => transpose). % % begin scalar ee,U,V,g,x,eps,tolerance,q,s,f,h,y,test_f_splitting, cancellation,test_f_convergence,convergence,c,z,denom,q_mat, I_rounded_turned_on,trans_done,val,val2,cols,cols2,tmpu ,tmpv; integer i,j,k,l,l1,m,n; trans_done := I_rounded_turned_on := nil; if not !*rounded then << on rounded; I_rounded_turned_on := t; >>; if not matrixp(A) then rederr "Error in spsvd: non matrix input."; % The value of eps can be decreased to increase accuracy. % As usual, doing this will slow things down (and vice versa). % It should not be made smaller than the value of rd!-tolerance!*. eps := get_num_part(my_reval({'times,1.5,{'expt,10,-8}})); tolerance := get_num_part(my_reval({'expt,10,-31})); % Algorithm requires m >= n. If this is not the case then transpose % the input and swap U and V in the output (as A = tp(U) diag(q) V % but tp(A) = tp(V) diag(q) U ). if sprow_dim(A) < spcol_dim(A) then << A := algebraic tp(A); trans_done := t; >>; m := sprow_dim(A); n := spcol_dim(A); U := copy_vect(A,nil); V := mkempspmat(n,list('spm,n,n)); ee := mkvect(n); q := mkvect(n); % Householder's reduction to bidiagonal form: g := x := 0; for i:=1:n do <<tmpu:=copy_vect(smtp (u,nil),nil); putv(ee,i,g); s := 0; l := i+1; cols:=findrow(tmpu,i); if cols then <<for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j>=i and j<=m then s := specrd!:plus(s,specrd!:expt(val,2)); >>; >>; if get_num_part(s) < tolerance then g := 0 else << f := findelem2(U,i,i); if get_num_part(f)<0 then g := specrd!:sqrt(s) else g := my_minus(specrd!:sqrt(s)); h := specrd!:plus(specrd!:times(f,g),my_minus(s)); letmtr3(list(U,i,i),specrd!:plus(f,my_minus(g)),u,nil); tmpu:=copy_vect(smtp (u,nil),nil); cols:=findrow(tmpu,i); for j:=l:n do <<cols2:=findrow(tmpu,j); s := 0; for each xx in cdr cols do <<val:=cdr xx; k:=car xx; if k>=i and k<=m then <<val2:=atsoc(k,cols2); if val2 then s := specrd!:plus(s,specrd!:times(val,cdr val2)); >>; >>; f := specrd!:quotient(s,h); for each xx in cdr cols do <<val:=cdr xx; k:=car xx; if k>=i and k<=m then <<val2:=atsoc(k,cols2); if val2=nil then val2:=(0 . 0); if not (f='(!:rd!: . 0.0)) then letmtr3(list(U,k,j),specrd!:plus(cdr val2, specrd!:times(f,val)),u,nil); >>; >>; >>; >>; putv(q,i,g); s := 0; cols:=findrow(u,i); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j>=l and j<=n then s := specrd!:plus(s,specrd!:expt(val,2)); >>; if get_num_part(s) < tolerance then g := 0 else <<f := findelem2(U,i,i+1); if get_num_part(f)<0 then g := specrd!:sqrt(s) else g := my_minus(specrd!:sqrt(s)); h := specrd!:plus(specrd!:times(f,g),my_minus(s)); letmtr3(list(U,i,i+1),specrd!:plus(f,my_minus(g)),u,nil); cols:=findrow(u,i); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j>=l and j<=n then putv(ee,j,specrd!:quotient(val,h)); >>; for j:=l:m do <<cols2:=findrow(u,j); s := 0; for each xx in cdr cols do <<val:=cdr xx; k:=car xx; if k>=l and k<=n then <<val2:=atsoc(k,cols2); if val2 then s := specrd!:plus(s,specrd!:times(val,cdr val2)); >>; >>; for each xx in cdr cols2 do <<k:=car xx; val2:=cdr xx; if k>=l and k<=n then <<val:=getv(ee,k); if val=nil then val:=0; letmtr3(list(U,j,k),specrd!:plus(val2, specrd!:times(s,val)),u,nil); >>; >>; >>; >>; y := specrd!:plus(abs(get_num_part(getv(q,i))), abs(get_num_part(getv(ee,i)))); if get_num_part(y) > get_num_part(x) then x := y; >>; % Accumulation of right hand transformations: for i:=n step -1 until 1 do << cols:=findrow(u,i); if get_num_part(g) neq 0 then <<val:=findelem2(u,i,i+1); h := specrd!:times(val,g); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j>=l and j<=n then letmtr3(list(V,j,i),specrd!:quotient(val,h),v,nil); >>; cols:=findrow(u,i); tmpv:=copy_vect(smtp(v,nil),nil); for j:=l:n do <<cols2:=findrow(tmpv,j); s:=0; for each xx in cdr cols do <<k:=car xx; val:=cdr xx; if k>=l and k<=n then <<val2:=atsoc(k,cols2); if val2 then s := specrd!:plus(s,specrd!:times(val,cdr val2)); >>; >>; cols:=findrow(tmpv,i); for each xx in cdr cols do <<k:=car xx; val:=cdr xx; if k>=l and k<=n then <<val2:=atsoc(k,cols2); if val2=nil then val2:=(0 . 0); letmtr3(list(V,k,j),specrd!:plus(cdr val2, specrd!:times(s,val)),v,nil); >>; >>; >>; >>; for j:=l:n do << letmtr3(list(V,i,j),0,v,nil); letmtr3(list(V,j,i),0,v,nil); >>; letmtr3(list(V,i,i),1,v,nil); g := getv(ee,i); l := i; >>; % Accumulation of left hand transformations: for i:=n step -1 until 1 do <<tmpu:=copy_vect(smtp (u,nil),nil); tmpv:=copy_vect(smtp (v,nil),nil); l := i+1; g := getv(q,i); cols:=findrow(u,i); for each xx in cdr cols do <<j:=car xx; if j>=l and j<=n then letmtr3(list(U,i,j),0,u,nil); >>; if get_num_part(g) neq 0 then <<h := specrd!:times(findelem2(U,i,i),g); tmpu:=copy_vect(smtp (u,nil),nil); cols:=findrow(tmpu,i); for j:=l:n do <<cols2:=findrow(tmpu,j); s := 0; for each xx in cdr cols do <<val:=cdr xx; k:=car xx; if k>=l and k<=m then <<val2:=atsoc(k,cols2); if val2 then s := specrd!:plus(s,specrd!:times(val,cdr val2)); >>; >>; f := specrd!:quotient(s,h); for each xx in cdr cols do <<val:=cdr xx; k:=car xx; if k>=i and k<=m then <<val2:=atsoc(k,cols2); if val2=nil then val2:=(0 . 0); if not (f='(minus (!:rd!: . 0.0))) then letmtr3(list(U,k,j),specrd!:plus(cdr val2, specrd!:times(f,val)),u,nil); >>; >>; >>; tmpu:=copy_vect(smtp (u,nil),nil); cols:=findrow(tmpu,i); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j>=i and j<=m then letmtr3(list(U,j,i),specrd!:quotient(val,g),u,nil); >>; >> else for each xx in cdr cols do << j:=car xx; if j>=i and j<=m then letmtr3(list(U,j,i),0,u,nil); >>; letmtr3(list(U,i,i),specrd!:plus(findelem2(U,i,i),1),u,nil); >>; % Diagonalisation of the bidiagonal form: eps := get_num_part(specrd!:times(eps,x)); test_f_splitting := t; k := n; while k>=1 do << convergence := nil; if test_f_splitting then << l := k; test_f_convergence := cancellation := nil; while l>=1 and not (test_f_convergence or cancellation) do << if abs(get_num_part(getv(ee,l))) <= eps then test_f_convergence := t else if abs(get_num_part(getv(q,l-1))) <= eps then cancellation := t else l := l-1; >>; >>; tmpu:=copy_vect(smtp (u,nil),nil); % Cancellation of e[l] if l>1: if not test_f_convergence then << c := 0; s := 1; l1 := l-1; i := l; while i<=k and not test_f_convergence do <<cols:=findrow(tmpu,i); f := specrd!:times(s,getv(ee,i)); putv(ee,i,specrd!:times(c,getv(ee,i))); if abs(get_num_part(f)) <= eps then test_f_convergence := t else << g := getv(q,i); h := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(g,g))); putv(q,i,h); c := specrd!:quotient(g,h); s := specrd!:quotient(my_minus(f),h); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j<=m then << y := findelem2(U,j,l1); letmtr3(list(U,j,l1),specrd!:plus(specrd!:times(y,c), specrd!:times(val,s)),u,nil); >>; >>; i := i+1; >>; >>; >>; z := getv(q,k); if l = k then convergence := t; if not convergence then << % Shift from bottom 2x2 minor: x := getv(q,l); y := getv(q,k-1); g := getv(ee,k-1); h := getv(ee,k); f := specrd!:quotient(specrd!:plus(specrd!:times( specrd!:plus(y,my_minus(z)),specrd!:plus(y,z)), specrd!:times(specrd!:plus(g,my_minus(h)), specrd!:plus(g,h))),specrd!:times( specrd!:times(2,h),y)); g := specrd!:sqrt(specrd!:plus(specrd!:times(f,f),1)); % Needed to change < here to <=. if get_num_part(f)<=0 then denom := specrd!:plus(f,my_minus(g)) else denom := specrd!:plus(f,g); f := specrd!:quotient(specrd!:plus(specrd!:times( specrd!:plus(x,my_minus(z)),specrd!:plus(x,z)), specrd!:times(h,specrd!:quotient(y, specrd!:plus(denom,my_minus(h))))),x); % Next QR transformation: c := s := 1; for i:=l+1:k do <<g := getv(ee,i); y := getv(q,i); h := specrd!:times(s,g); g := specrd!:times(c,g); z := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(h,h))); putv(ee,i-1,z); c := specrd!:quotient(f,z); s := specrd!:quotient(h,z); f := specrd!:plus(specrd!:times(x,c),specrd!:times(g,s)); g := specrd!:plus(specrd!:times(my_minus(x),s), specrd!:times(g,c)); h := specrd!:times(y,s); y := specrd!:times(y,c); for j:=1:m do << z := findelem2(V,j,i); x :=findelem2(v,j,i-1); letmtr3(list(V,j,i-1),specrd!:plus(specrd!:times(x,c), specrd!:times(z,s)),v,nil); letmtr3(list(V,j,i),specrd!:difference(specrd!:times (z,c), specrd!:times(x,s)),v,nil); >>; z := specrd!:sqrt(specrd!:plus(specrd!:times(f,f), specrd!:times(h,h))); putv(q,i-1,z); c := specrd!:quotient(f,z); s := specrd!:quotient(h,z); f := specrd!:plus(specrd!:times(c,g),specrd!:times(s,y)); x := specrd!:plus(specrd!:times(my_minus(s),g), specrd!:times(c,y)); for j:=1:m do << y := findelem2(U,j,i-1); z := findelem2(u,j,i); letmtr3(list(U,j,i-1),specrd!:plus(specrd!:times(y,c), specrd!:times(z,s)),u,nil); letmtr3(list(U,j,i),specrd!:difference(specrd!:times(z,c), specrd!:times(y,s)),u,nil); >>; >>; putv(ee,l,0); putv(ee,k,f); putv(q,k,x); >> else % convergence: <<tmpv:=copy_vect(smtp (v,nil),nil); if get_num_part(z)<0 then << % q[k] is made non-negative: putv(q,k,my_minus(z)); cols:=findrow(tmpv,k); for each xx in cdr cols do <<j:=car xx; val:=cdr xx; if j<=n then letmtr3(list(V,j,k),my_minus(val),v,nil); >>; >>; k := k-1; >>; >>; q_mat := spq_to_diag_matrix(q); if I_rounded_turned_on then off rounded; v:=spden_to_sp(v); % to print it out in Sparse manner u:=spden_to_sp(u); if trans_done then return {'list,algebraic tp V,q_mat,algebraic tp U} else return {'list,algebraic tp U,q_mat,algebraic tp V}; end; flag('(spsvd),'opfn); % To make it available from algebraic (user) mode. symbolic procedure spq_to_diag_matrix(q); % % Converts q (a vector) to a diagonal matrix with the elements of % q on the diagonal. % begin scalar q_mat; integer i,sq_dim_q,val; sq_dim_q := upbv(q); q_mat := mkempspmat(sq_dim_q,list('spm,sq_dim_q,sq_dim_q)); for i:=1:sq_dim_q do << val:=getv(q,i); if val='(!:rd!: . 0.0) then nil else letmtr3(list(q_mat,i,i),val,q_mat,nil); >>; return q_mat; end; % The lists are then re-written into desired sparse list format ready % for printing. symbolic procedure spden_to_sp(list); begin scalar tl,nmat,val,cols,j; tl:=caddr list; nmat:=mkempspmat(cadr tl,tl); for i:=1:cadr tl do << cols:=findrow(list,i); for each xx in cdr cols do <<j:=car xx; val:=reval cdr xx; if val='(!:rd!: . 0.0) then nil else letmtr3(list(nmat,i,j),val,nmat,nil); >>; >>; return nmat; end; symbolic procedure sprd_copy_mat(A); % % Creates a copy of the input matrix and returns it as well as % reval-ing each elt to get them in !:rd!: form; % begin scalar C; integer row_dim,column_dim; row_dim := sprow_dim(A); column_dim := spcol_dim(A); C := list('sparsemat,list('u,row_dim,column_dim)); for i:=1:row_dim do << for j:=1:column_dim do << letmtr3(list(C,i,j),my_reval(findelem2(A,i,j)),c,nil); >>; >>; return C; end; symbolic procedure sppseudo_inverse(in_mat); % % Also known as the Moore-Penrose Inverse. % % Given the singular value decomposition A := tp(U) diag(q) V % the pseudo inverse A^(-1) is defined as % % A^(-1) = tp(V) (diag(q))^(-1) U. % % NB: this can be quite handy as we can take the inverse of non % square matrices (A * pseudo_inverse(A) = identity). % begin scalar psu_inv,svd_list,a,b,c; svd_list := cdr spsvd(in_mat); a:=car svd_list; c:=caddr svd_list; b:=cadr svd_list; c:=algebraic tp c; b:=algebraic (1/b); psu_inv := algebraic (c * b * a); return psu_inv; end; flag('(sppseudo_inverse),'opfn); endmodule; end; %*********************************************************************** %======================================================================= % % End of Code. % %======================================================================= %***********************************************************************