(FILECREATED "30-Jan-86 13:10:23" {QV}<PEDERSEN>LISP>VMATRIXOPS.;2 8483
changes to: (VARS VMATRIXOPSCOMS MATRIXOPSCOMS)
(FNS MAKESAVINGS MAKESTACK COL2NORM COLDOTPROD COLROT COLSCALE COLSWAP
SVDNASH.SAVE SVDNASH SVDNASH.FAST SVDTEST RSOLV)
previous date: "28-Jan-86 23:36:57" {QV}<PEDERSEN>LISP>MATRIXOPS.;1)
(PRETTYCOMPRINT VMATRIXOPSCOMS)
(RPAQQ VMATRIXOPSCOMS ((FNS MAKESAVINGS MAKESTACK RSOLV SVDNASH SVDTEST)
(VARS STACK)))
(DEFINEQ
(MAKESAVINGS
[LAMBDA (LST) (* jop: "29-Jan-86 12:34")
(MAKE-ARRAY (QUOTE (50 5))
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)
(QUOTE :INITIAL-CONTENTS)
LST])
(MAKESTACK
[LAMBDA (LST) (* jop: "28-Jan-86 22:56")
(MAKE-ARRAY (QUOTE (21 3))
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)
(QUOTE :INITIAL-CONTENTS)
LST])
(RSOLV
[LAMBDA (RMATRIX BVECTOR CVECTOR) (* jop: "29-Jan-86 17:45")
(* * Calcluate the solution vector BVECTOR for the system of linear equations R*B=C, where RMATRIX is upper
triangular with non-zero diagonal elements. The vector BVECTOR must be created by the CALLER!)
(LET* ((N (ARRAY-DIMENSION RMATRIX 0))
(M (ARRAY-DIMENSION RMATRIX 1))
(INDEXLIMIT (SUB1 N)))
(* * Solution by backsubstitution.)
(BLAS.ARRAYBLT CVECTOR 0 1 BVECTOR 0 1 N)
(LASET (FQUOTIENT (LAREF BVECTOR INDEXLIMIT)
(LAREF RMATRIX INDEXLIMIT INDEXLIMIT))
B INDEXLIMIT)
[if (IGREATERP N 1)
then (bind J FTEMP for JJ from 1 to INDEXLIMIT
do (SETQ J (IDIFFERENCE N JJ))
(SETQ FTEMP (FMINUS (LAREF B J)))
(BLAS.AXPY FTEMP RMATRIX J M BVECTOR 0 1 J)
(LASET (FQUOTIENT (LAREF BVECTOR (SUB1 J))
(LAREF RMATRIX (SUB1 J)
(SUB1 J)))
BVECTOR
(SUB1 J]
BVECTOR])
(SVDNASH
[LAMBDA (UMATRIX SVECTOR VMATRIX) (* jop: "29-Jan-86 14:39")
(* * Singular-value decomposition by means of orthogonalization by plane rotations. Taken from Nash and Shlien:
"Partial svd algorithms." On entry UMATRIX contains the m by n matrix to be decomposed, SVECTOR must be a vector of
length n, and VMATRIX must be a square n by n matrix. On return UMATRIX has been overwritten by the left singular
vectors, SVECTOR contains the singular values, and VMATRIX contains the right singular vectors.)
(LET* ((M (ARRAY-DIMENSION UMATRIX 0))
(N (ARRAY-DIMENSION UMATRIX 1))
(EPS 1.0E-6)
(SLIMIT (IMAX (IQUOTIENT N 4)
6))
(NT N)
(TOL EPS))
(* * Initialize VMATRIX to identity matrix.)
(BLAS.ARRAYFILL 0.0 VMATRIX)
(for I from 0 to (SUB1 N) do (LASET 1.0 VMATRIX I I))
(* * The main loop: repeatedly sweep over all pairs of columns in U, rotating as needed, until no rotations in a
complete sweep are effective. Check the opportunity for rank reduction at the conclusion of each sweep.)
[bind (SCOUNT ← 0)
(EPSSQUARED ←(FTIMES EPS EPS))
RCOUNT repeatwhile (IGREATERP RCOUNT 0)
eachtime (SETQ RCOUNT (IQUOTIENT (ITIMES NT (SUB1 NT))
2))
(SETQ SCOUNT (ADD1 SCOUNT))
do (if (IGREATERP SCOUNT SLIMIT)
then (HELP "Number of sweeps exceeds sweep limit." SCOUNT))
[for J from 0 to (IDIFFERENCE NT 2)
do (for K from (IPLUS J 1) to (SUB1 NT)
bind P Q R C S V
do (SETQ P (BLAS.DOTPROD UMATRIX J N UMATRIX K N M))
(SETQ Q (BLAS.DOTPROD UMATRIX J N UMATRIX J N M))
(SETQ R (BLAS.DOTPROD UMATRIX K N UMATRIX K N M))
(LASET Q SVECTOR J)
(LASET R SVECTOR K)
(LET ((FP P)
(FQ Q)
(FR R)
FV FC FS)
(DECLARE (TYPE FLOATP FP FQ FR FV FC FS))
(if (FLESSP Q R)
then (SETQ FQ (FDIFFERENCE (FQUOTIENT FQ FR)
1.0))
(SETQ FP (FQUOTIENT FP FR))
[SETQ FV (SQRT (SETQ FV
(FPLUS (FTIMES 4.0 FP FP)
(FTIMES FQ FQ]
[SETQ FS (SQRT (SETQ FS
(FTIMES .5
(FDIFFERENCE
1.0
(FQUOTIENT
FQ FV]
(if (FLESSP FP 0.0)
then (SETQ FS (FDIFFERENCE 0.0 FS)))
[SETQ C (SETQ FC (FQUOTIENT FP
(FTIMES FV FS]
(SETQ S FS)
(BLAS.ROT C S UMATRIX J N UMATRIX K N M)
(BLAS.ROT C S VMATRIX J N VMATRIX K N N)
elseif (OR (LEQ (FTIMES Q R)
EPSSQUARED)
(LEQ (FTIMES (FQUOTIENT P Q)
(FQUOTIENT P R))
EPS))
then (SETQ RCOUNT (SUB1 RCOUNT))
else (SETQ FR (FDIFFERENCE 1.0 (FQUOTIENT FR FQ)))
(SETQ FP (FQUOTIENT FP FQ))
[SETQ FV (SQRT (SETQ FV
(FPLUS (FTIMES 4.0 FP FP)
(FTIMES FR FR]
[SETQ FC (SQRT (SETQ FC
(FTIMES .5
(FPLUS 1.0
(FQUOTIENT
FR FV]
[SETQ S (SETQ FS (FQUOTIENT FP
(FTIMES FV FC]
(SETQ C FC)
(BLAS.ROT C S UMATRIX J N UMATRIX K N M)
(BLAS.ROT C S VMATRIX J N VMATRIX K N N]
(while (AND (GEQ NT 3)
(LEQ (FQUOTIENT (LAREF SVECTOR (SUB1 NT))
(FPLUS (LAREF SVECTOR 0)
TOL))
TOL))
do (SETQ NT (SUB1 NT]
(* * Finish the decomposition by returning all N singular values, and by normalizing those columns of UMATRIX
judged to be non-zero.)
(for J from 0 to (SUB1 N) bind Q
do (SETQ Q (SQRT (LAREF SVECTOR J)))
(LASET Q SVECTOR J)
(if (LEQ J NT)
then (BLAS.SCAL (FQUOTIENT 1.0 Q)
UMATRIX J N M])
(SVDTEST
[LAMBDA NIL (* jop: "29-Jan-86 14:06")
(* * comment)
(LET ((UU (MAKE-ARRAY (QUOTE (24 19))
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
(SS (MAKE-ARRAY 19 (QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
(VV (MAKE-ARRAY (QUOTE (19 19))
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
(L 0))
[for TT from 1.0 to 3.0 do (for PP from 1.0 to 2.0
do (for CC from 1.0 to 4.0
do (ASET TT UU L 0)
(ASET PP UU L 1)
(ASET CC UU L 2)
(ASET (FTIMES TT PP)
UU L 3)
(ASET (FTIMES TT CC)
UU L 4)
(ASET (FTIMES PP CC)
UU L 5)
(ASET (FTIMES PP PP)
UU L 6)
(ASET (FTIMES CC CC)
UU L 7)
(ASET (FTIMES TT TT)
UU L 8)
(ASET (FTIMES TT PP PP)
UU L 9)
(ASET (FTIMES TT CC CC)
UU L 10)
(ASET (FTIMES PP TT TT)
UU L 11)
(ASET (FTIMES PP CC CC)
UU L 12)
(ASET (FTIMES CC TT TT)
UU L 13)
(ASET (FTIMES CC PP PP)
UU L 14)
(ASET (FTIMES TT TT TT)
UU L 15)
(ASET (FTIMES PP PP PP)
UU L 16)
(ASET (FTIMES CC CC CC)
UU L 17)
(ASET (FTIMES TT PP CC)
UU L 18)
(SETQ L (ADD1 L]
(TIMEALL (SVDNASH UU SS VV])
)
(RPAQQ STACK ((80 27 89)
(80 27 88)
(75 25 90)
(62 24 87)
(62 22 87)
(62 23 87)
(62 24 93)
(62 24 93)
(58 23 87)
(58 18 80)
(58 18 89)
(58 17 88)
(58 18 82)
(58 19 93)
(50 18 89)
(50 18 86)
(50 19 72)
(50 19 79)
(50 20 80)
(56 20 82)
(70 20 91)))
(DECLARE: DONTCOPY
(FILEMAP (NIL (473 8169 (MAKESAVINGS 483 . 726) (MAKESTACK 728 . 969) (RSOLV 971 . 2085) (SVDNASH 2087
. 6399) (SVDTEST 6401 . 8167)))))
STOP