(DEFINE-FILE-INFO READTABLE "INTERLISP" PACKAGE "INTERLISP")
(FILECREATED " 5-Dec-88 19:07:01" {QV}<IDL>KOTO>MATRIXOPS.;6 29914
changes to%: (FNS SVDFACTOR QRFACTOR QRQTY MREGRESS MSOLVE LUINVERSE LUSOLV LUFACTOR LSOLV \FLOATAREFMACRO \FLOATASETMACRO MTRANSPOSE MTIMES)
previous date%: "23-Jun-86 14:22:28" {QV}<IDL>KOTO>MATRIXOPS.;1)
(* "
Copyright (c) 1986, 1988 by Xerox Corporation. All rights reserved.
")
(PRETTYCOMPRINT MATRIXOPSCOMS)
(RPAQQ MATRIXOPSCOMS ((FNS CHOLESKYFACTOR MTRANSPOSE MINVERT LSOLV LUFACTOR LUINVERSE LUSOLV MTIMES QRFACTOR QROLS QRQTY QRQY QRSOLV MREGRESS RSOLV MSOLVE SVDFACTOR SVDTEST \FLOATAREFMACRO \FLOATASETMACRO) (VARS STACK) (MACROS \FLOATAREF \FLOATASET) (FILES BLAS)))
(DEFINEQ
(CHOLESKYFACTOR
(LAMBDA (MATRIX FACTORMATRIX) (* jop%: "23-Jun-86 13:41") (* * Lifted from LINPACK algorithm SCHDC) (BLAS.CHECKARRAY MATRIX) (LET ((P (CL:ARRAY-DIMENSION MATRIX 0))) (* * Arg Checks) (if (NOT (AND (EQ 2 (CL:ARRAY-RANK MATRIX)) (EQ P (CL:ARRAY-DIMENSION MATRIX 1)))) then (HELP "Matrix not sqaure" MATRIX)) (if (NULL FACTORMATRIX) then (SETQ FACTORMATRIX (CL:MAKE-ARRAY (CL:ARRAY-DIMENSIONS MATRIX) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY FACTORMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS FACTORMATRIX) (CL:ARRAY-DIMENSIONS MATRIX))) then (HELP "Illegal FACTORMATRIX" FACTORMATRIX))) (* Copy MATRIX to FACTORMATRIX) (BLAS.ARRAYBLT MATRIX 0 1 FACTORMATRIX 0 1) (* * Compute the cholesky decomposition of FACTORMATRIX) (bind (WORK ← (CL:MAKE-ARRAY P (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) TEMP for K from 0 to (SUB1 P) do (if (LEQ (\FLOATAREF FACTORMATRIX K K) 0.0) then (HELP "Zero pivot element")) (\FLOATASET (SQRT (\FLOATAREF FACTORMATRIX K K)) WORK K) (\FLOATASET (\FLOATAREF WORK K) FACTORMATRIX K K) (if (NOT (EQ K (SUB1 P))) then (for J from (ADD1 K) to (SUB1 P) do (\FLOATASET (FQUOTIENT (\FLOATAREF FACTORMATRIX K J) (\FLOATAREF WORK K)) FACTORMATRIX K J) (\FLOATASET (\FLOATAREF FACTORMATRIX K J) WORK J) (SETQ TEMP (FMINUS (\FLOATAREF FACTORMATRIX K J))) (BLAS.AXPY TEMP WORK (ADD1 K) 1 FACTORMATRIX (IPLUS J (ITIMES P (ADD1 K))) P (IDIFFERENCE J K))))) FACTORMATRIX))
)
(MTRANSPOSE
(LAMBDA (SOURCEMATRIX DESTMATRIX) (* jop%: " 4-Jun-86 14:07") (* * Transpose the M x N matrix SOURCEMATRIX. DESTMATRIX should be N x M. Returns DESTMATRIX) (BLAS.CHECKARRAY SOURCEMATRIX) (PROG ((M (CL:ARRAY-DIMENSION SOURCEMATRIX 0)) (N (CL:ARRAY-DIMENSION SOURCEMATRIX 1))) (if (NULL DESTMATRIX) then (SETQ DESTMATRIX (CL:MAKE-ARRAY (LIST N M) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY DESTMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS DESTMATRIX) (LIST N M))) then (HELP "DESTMATRIX of incorrect size" DESTMATRIX))) (if (ILESSP M N) then (bind (SOURCEBASE ← (%%ARRAY-BASE SOURCEMATRIX)) (DESTBASE ← (%%ARRAY-BASE DESTMATRIX)) for I from 0 to (SUB1 M) do (\FLOATARRAYBLT SOURCEBASE (ITIMES N I) 1 DESTBASE I M N)) else (bind (SOURCEBASE ← (%%ARRAY-BASE SOURCEMATRIX)) (DESTBASE ← (%%ARRAY-BASE DESTMATRIX)) for J from 0 to (SUB1 N) do (\FLOATARRAYBLT SOURCEBASE J N DESTBASE (ITIMES J M) 1 M))) (RETURN DESTMATRIX)))
)
(MINVERT
(LAMBDA (MATRIX SOLUTION) (* jop%: "26-May-86 18:35") (* * Solves to system A x = b. BVECTOR should to the RHS of the system. Returns SOLUTION) (LET ((PIVOTVECTOR (CL:MAKE-ARRAY (CL:ARRAY-DIMENSION MATRIX 0)))) (LUINVERSE (LUFACTOR MATRIX PIVOTVECTOR) PIVOTVECTOR SOLUTION)))
)
(LSOLV
(LAMBDA (LMATRIX CVECTOR BVECTOR) (* ; "Edited 4-Dec-88 18:21 by jop") (* ;; "Calcluate the solution vector BVECTOR for the system of linear equations R*B=C, where LMATRIX is lower triangular M X N with non-zero diagonal elements. BVECTOR and CVECTOR must be of size N. Always returns BVECTOR") (BLAS.CHECKARRAY LMATRIX) (BLAS.CHECKARRAY CVECTOR) (PROG ((M (CL:ARRAY-DIMENSION LMATRIX 0)) (N (CL:ARRAY-DIMENSION LMATRIX 1))) (* ;;; "Arg Checks") (if (ILESSP M N) then (HELP "Order of system less than" N)) (if (NOT (EQ (CL:ARRAY-TOTAL-SIZE CVECTOR) N)) then (HELP "CVECTOR not of size" N)) (if (NULL BVECTOR) then (SETQ BVECTOR (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY BVECTOR) (if (NOT (EQ (CL:ARRAY-TOTAL-SIZE BVECTOR) N)) then (HELP "BVECTOR not of size" N))) (* ;; "Check for zero diagonal elements") (if (for I from 0 to (SUB1 N) thereis (UFEQP 0.0 (\FLOATAREF LMATRIX I I))) then (HELP "LMATRIX has a zero diagonal element")) (* ;; "Solution by forward substitution") (* ;; "Copy CVECTOR to BVECTOR") (BLAS.ARRAYBLT CVECTOR 0 1 BVECTOR 0 1 N) (* ;; "Compute the first value") (\FLOATASET (FQUOTIENT (\FLOATAREF BVECTOR 0) (\FLOATAREF LMATRIX 0 0)) BVECTOR 0) (for J from 1 to (SUB1 N) do (BLAS.AXPY (FMINUS (\FLOATAREF BVECTOR (SUB1 J))) LMATRIX (IPLUS (SUB1 J) (ITIMES J N)) N BVECTOR J 1 (IDIFFERENCE N J)) (\FLOATASET (FQUOTIENT (\FLOATAREF BVECTOR J) (\FLOATAREF LMATRIX J J)) BVECTOR J)) (RETURN BVECTOR)))
)
(LUFACTOR
(LAMBDA (MATRIX PIVOTVECTOR FACTORMATRIX) (* ; "Edited 4-Dec-88 18:26 by jop") (* ;;; "Computes the LU decomposition of the N x N matrix MATRIX by Gauss elimination with row pivoting. FACTORMATRIX will be overwritten with the packed result. PIVOTVECTOR will be a vector of smallposp's, holding the pivot permutation, and must be supplied. Returns NIL in the normal case, else returns the row index") (* ;;; "Lifted from LINPACK algorithm SGESL") (BLAS.CHECKARRAY MATRIX) (if (NOT (AND (CL:ARRAYP PIVOTVECTOR) (EQ (CL:ARRAY-ELEMENT-TYPE PIVOTVECTOR) T))) then (HELP "Must be a pointer array" PIVOTVECTOR)) (LET ((N (CL:ARRAY-DIMENSION MATRIX 0))) (* ;;; "Arg Checks") (if (AND (EQ 2 (CL:ARRAY-RANK MATRIX)) (NOT (EQ N (CL:ARRAY-DIMENSION MATRIX 1)))) then (HELP "MATRIX not square" MATRIX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK PIVOTVECTOR)) (EQ N (CL:ARRAY-TOTAL-SIZE PIVOTVECTOR)))) then (HELP "PIVOTVECTOR not of size N" PIVOTVECTOR)) (if (NULL FACTORMATRIX) then (SETQ FACTORMATRIX (CL:MAKE-ARRAY (LIST N N) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY FACTORMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS FACTORMATRIX) (CL:ARRAY-DIMENSIONS MATRIX))) then (HELP "Illegal FACTORMATRIX" FACTORMATRIX))) (* ; "Copy MATRIX to FACTORMATRIX") (BLAS.ARRAYBLT MATRIX 0 1 FACTORMATRIX 0 1) (* ;;; "Compute the LU decomposition of FACTORMATRIX") (bind PIVOTINDEX TEMP for K from 0 to (IDIFFERENCE N 2) do (* ; "find pivot index") (SETQ PIVOTINDEX (IPLUS (BLAS.MAX FACTORMATRIX (IPLUS K (ITIMES N K)) N (IDIFFERENCE N K)) K)) (ASET PIVOTINDEX PIVOTVECTOR K) (if (NOT (FEQP (\FLOATAREF FACTORMATRIX PIVOTINDEX K) 0.0)) then (if (NOT (EQ PIVOTINDEX K)) then (* ; "Interchange") (SETQ TEMP (\FLOATAREF FACTORMATRIX PIVOTINDEX K)) (\FLOATASET (\FLOATAREF FACTORMATRIX K K) FACTORMATRIX PIVOTINDEX K) (\FLOATASET TEMP FACTORMATRIX K K)) (* ; "compute Multpliers") (BLAS.SCAL (FMINUS (FQUOTIENT 1.0 (\FLOATAREF FACTORMATRIX K K))) FACTORMATRIX (IPLUS K (ITIMES N (ADD1 K))) N (SUB1 (IDIFFERENCE N K))) (* ; "Row eliminate with column indexing") (bind (KPLUS1 ← (ADD1 K)) for J from (ADD1 K) to (SUB1 N) do (SETQ TEMP (\FLOATAREF FACTORMATRIX PIVOTINDEX J)) (if (NOT (EQ PIVOTINDEX K)) then (* ; "Interchange") (\FLOATASET (\FLOATAREF FACTORMATRIX K J) FACTORMATRIX PIVOTINDEX J) (\FLOATASET TEMP FACTORMATRIX K J)) (BLAS.AXPY TEMP FACTORMATRIX (IPLUS K (ITIMES N KPLUS1)) N FACTORMATRIX (IPLUS J (ITIMES N KPLUS1)) N (IDIFFERENCE N KPLUS1))))) (* ; "No row elimination on last column") (ASET (SUB1 N) PIVOTVECTOR (SUB1 N)) FACTORMATRIX))
)
(LUINVERSE
(LAMBDA (LUMATRIX PIVOTVECTOR SOLUTION) (* ; "Edited 4-Dec-88 18:26 by jop") (* ;;; "Forms MATRIX inverse where LUMATRIX and PIVOTVECTOR are the outputs of LUFACTOR.") (* ;;; "lifted from LINPACK SGEDI") (BLAS.CHECKARRAY LUMATRIX) (if (NOT (AND (CL:ARRAYP PIVOTVECTOR) (EQ (CL:ARRAY-ELEMENT-TYPE PIVOTVECTOR) T))) then (HELP "Must be an array of pointers" PIVOTVECTOR)) (PROG ((N (CL:ARRAY-DIMENSION LUMATRIX 0))) (* ;;; "Arg Checks") (if (AND (EQ 2 (CL:ARRAY-RANK LUMATRIX)) (NOT (EQ N (CL:ARRAY-DIMENSION LUMATRIX 1)))) then (HELP "MATRIX not square" LUMATRIX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK PIVOTVECTOR)) (EQ N (CL:ARRAY-TOTAL-SIZE PIVOTVECTOR)))) then (HELP "PIVOTVECTOR not a vector of size N" PIVOTVECTOR)) (if (NULL SOLUTION) then (SETQ SOLUTION (CL:MAKE-ARRAY (LIST N N) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY SOLUTION) (if (NOT (AND (EQ 2 (CL:ARRAY-RANK SOLUTION)) (EQUAL (CL:ARRAY-DIMENSIONS LUMATRIX) (CL:ARRAY-DIMENSIONS SOLUTION)))) then (HELP "SOLUTION not an N x N array" SOLUTION))) (* ; "copy LUMATRIX to SOLUTION") (BLAS.ARRAYBLT LUMATRIX 0 1 SOLUTION 0 1) (* ;;; "first compute INVERSE (U)") (bind TEMP for K from 0 to (SUB1 N) do (\FLOATASET (FQUOTIENT 1.0 (\FLOATAREF SOLUTION K K)) SOLUTION K K) (SETQ TEMP (FMINUS (\FLOATAREF SOLUTION K K))) (BLAS.SCAL TEMP SOLUTION K N K) (bind TEMP for J from (ADD1 K) to (SUB1 N) do (SETQ TEMP (\FLOATAREF SOLUTION K J)) (\FLOATASET 0.0 SOLUTION K J) (BLAS.AXPY TEMP SOLUTION K N SOLUTION J N (ADD1 K)))) (* ;;; "Form INVERSE (U) *INVERSE (L)") (bind (TEMPARRAY ← (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) L for K from (IDIFFERENCE N 2) to 0 by -1 do (for I from (ADD1 K) to (SUB1 N) do (\FLOATASET (\FLOATAREF SOLUTION I K) TEMPARRAY I) (\FLOATASET 0.0 SOLUTION I K)) (bind TEMP for J from (ADD1 K) to (SUB1 N) do (SETQ TEMP (\FLOATAREF TEMPARRAY J)) (BLAS.AXPY TEMP SOLUTION J N SOLUTION K N N)) (SETQ L (CL:AREF PIVOTVECTOR K)) (if (NEQ L K) then (BLAS.SWAP SOLUTION K N SOLUTION L N N))) (RETURN SOLUTION)))
)
(LUSOLV
(LAMBDA (LUMATRIX PIVOTVECTOR CVECTOR SOLUTION) (* ; "Edited 4-Dec-88 18:26 by jop") (* ;;; "Solves to system A x = b. LUMATRIX and PIVOTVECTOR should be the outputs of LUFACTOR. CVECTOR should to the RHS of the system. Returns SOLUTION") (* ;;; "lifted from LINPACK SGESL") (BLAS.CHECKARRAY LUMATRIX) (if (NOT (AND (CL:ARRAYP PIVOTVECTOR) (EQ (CL:ARRAY-ELEMENT-TYPE PIVOTVECTOR) T))) then (HELP "Must be an array of pointers" PIVOTVECTOR)) (BLAS.CHECKARRAY CVECTOR) (PROG ((N (CL:ARRAY-DIMENSION LUMATRIX 0))) (* ;;; "Arg Checks") (if (AND (EQ 2 (CL:ARRAY-RANK LUMATRIX)) (NOT (EQ N (CL:ARRAY-DIMENSION LUMATRIX 1)))) then (HELP "MATRIX not square" LUMATRIX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK PIVOTVECTOR)) (EQ N (CL:ARRAY-TOTAL-SIZE PIVOTVECTOR)))) then (HELP "PIVOTVECTOR not a vector of size N" PIVOTVECTOR)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK CVECTOR)) (EQ N (CL:ARRAY-TOTAL-SIZE CVECTOR)))) then (HELP "CVECTOR not a vector of size N" CVECTOR)) (if (NULL SOLUTION) then (SETQ SOLUTION (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY SOLUTION) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK SOLUTION)) (EQ N (CL:ARRAY-TOTAL-SIZE SOLUTION)))) then (HELP "SOLUTION not avector of size N" SOLUTION))) (* ; "Copy CVECTOR to SOLUTION") (BLAS.ARRAYBLT CVECTOR 0 1 SOLUTION 0 1 N) (* ;;; "First solve L*y = b") (bind PIVOTINDEX TEMP for K from 0 to (IDIFFERENCE N 2) do (SETQ PIVOTINDEX (CL:AREF PIVOTVECTOR K)) (SETQ TEMP (\FLOATAREF SOLUTION PIVOTINDEX)) (if (NOT (EQ PIVOTINDEX K)) then (* ; "interchange") (\FLOATASET (\FLOATAREF SOLUTION K) SOLUTION PIVOTINDEX) (\FLOATASET TEMP SOLUTION K)) (BLAS.AXPY TEMP LUMATRIX (IPLUS K (ITIMES N (ADD1 K))) N SOLUTION (ADD1 K) 1 (IDIFFERENCE N (ADD1 K)))) (* ;;; "Then solve U*x = y") (bind TEMP for K from (SUB1 N) to 0 by -1 do (SETQ TEMP (FMINUS (\FLOATASET (FQUOTIENT (\FLOATAREF SOLUTION K) (\FLOATAREF LUMATRIX K K)) SOLUTION K))) (BLAS.AXPY TEMP LUMATRIX K N SOLUTION 0 1 K)) (RETURN SOLUTION)))
)
(MTIMES
(LAMBDA (A B PRODUCT) (* jop%: "17-Jun-86 18:14") (* * Matrix multiply. A may be an N vector or a (M x N) matrix and B may be a N vector or a N x P matrix. PRODUCT defaults to a M x P array. RETURNS PRODUCT) (BLAS.CHECKARRAY A) (BLAS.CHECKARRAY B) (LET* ((RANKA (CL:ARRAY-RANK A)) (ADIMS (CL:ARRAY-DIMENSIONS A)) (RANKB (CL:ARRAY-RANK B)) (BDIMS (CL:ARRAY-DIMENSIONS B)) (M (if (EQ RANKA 1) then 1 else (CAR ADIMS))) (N (if (EQ RANKA 1) then (CAR ADIMS) else (CADR ADIMS))) (P (if (EQ RANKB 1) then 1 else (CADR BDIMS)))) (* * Args checks) (if (NOT (OR (EQ RANKA 1) (EQ RANKA 2))) then (HELP "A not a one-d or two-d array" A)) (if (NOT (OR (EQ RANKB 1) (EQ RANKB 2))) then (HELP "B not a one-d or two-d array" B)) (if (NOT (EQ (CAR BDIMS) N)) then (HELP "Args not conformable")) (* * Do the multiply) (if (AND (EQ M 1) (EQ P 1)) then (bind (ABASE ← (%%ARRAY-BASE A)) (BBASE ← (%%ARRAY-BASE B)) (FTEMP ← 0.0) declare (TYPE FLOATP FTEMP) for I from 0 to (MUL2 (SUB1 N)) by 2 do (SETQ FTEMP (FPLUS FTEMP (FTIMES (\GETBASEFLOATP ABASE I) (\GETBASEFLOATP BBASE I)))) finally (RETURN FTEMP)) else (if (NULL PRODUCT) then (SETQ PRODUCT (CL:MAKE-ARRAY (APPEND (LDIFF ADIMS (LAST ADIMS)) (CDR BDIMS)) (QUOTE :ELEMENT-TYPE) (QUOTE CL:SINGLE-FLOAT))) else (BLAS.CHECKARRAY PRODUCT) (if (NOT (AND (EQ (CL:ARRAY-RANK PRODUCT) (if (OR (EQ M 1) (EQ P 1)) then 1 else 2)) (if (EQ M 1) then (EQ P (CL:ARRAY-DIMENSION PRODUCT 0)) else (if (EQ P 1) then (EQ M (CL:ARRAY-DIMENSION PRODUCT 0)) else (AND (EQ M (CL:ARRAY-DIMENSION PRODUCT 0)) (EQ P (CL:ARRAY-DIMENSION PRODUCT 1))))))) then (HELP "Product of incorrect size" PRODUCT))) (bind (ABASE ← (%%ARRAY-BASE A)) (BBASE ← (%%ARRAY-BASE B)) (CBASE ← (%%ARRAY-BASE PRODUCT)) for I from 0 to (SUB1 M) do (for J from 0 to (SUB1 P) as COFFSET from (MUL2 (ITIMES P I)) by 2 do (bind (FTEMP ← 0.0) declare (TYPE FLOATP FTEMP) for K from 0 to (SUB1 N) as AOFFSET from (MUL2 (ITIMES N I)) by 2 as BOFFSET from (MUL2 J) by (MUL2 P) do (SETQ FTEMP (FPLUS FTEMP (FTIMES (\GETBASEFLOATP ABASE AOFFSET) (\GETBASEFLOATP BBASE BOFFSET)))) finally (\PUTBASEFLOATP CBASE COFFSET FTEMP)))) PRODUCT)))
)
(QRFACTOR
(LAMBDA (MATRIX QRAUX FACTORMATRIX) (* ; "Edited 4-Dec-88 18:31 by jop") (* ;;; "Computes the LU decomposition of the N x N matrix MATRIX by Gauss elimination with row pivoting. FACTORMATRIX will be overwritten with the packed result. QRAUX will be a vector of smallposp's, holding the pivot permutation, and must be supplied. Returns NIL in the normal case, else returns the row index") (* ;;; "Lifted from LINPACK algorithm SGESL") (BLAS.CHECKARRAY MATRIX) (BLAS.CHECKARRAY QRAUX) (LET ((N (CL:ARRAY-DIMENSION MATRIX 0)) (P (CL:ARRAY-DIMENSION MATRIX 1))) (* ;;; "Arg Checks") (if (NOT (AND (EQ 1 (CL:ARRAY-RANK QRAUX)) (EQ P (CL:ARRAY-TOTAL-SIZE QRAUX)))) then (HELP "QRAUX not of size P" QRAUX)) (if (NULL FACTORMATRIX) then (SETQ FACTORMATRIX (CL:MAKE-ARRAY (CL:ARRAY-DIMENSIONS MATRIX) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY FACTORMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS FACTORMATRIX) (CL:ARRAY-DIMENSIONS MATRIX))) then (HELP "Illegal FACTORMATRIX" FACTORMATRIX))) (* ; "Copy MATRIX to FACTORMATRIX") (BLAS.ARRAYBLT MATRIX 0 1 FACTORMATRIX 0 1) (* ;;; "Compute the QR decomposition of FACTORMATRIX") (for I from 0 to (SUB1 P) do (\FLOATASET 0.0 QRAUX I)) (bind NRMXL for L from 0 to (SUB1 (IMIN N P)) unless (EQ L (SUB1 N)) do (* ; "Compute the Householder transformation for column L") (SETQ NRMXL (BLAS.NRM2 FACTORMATRIX (IPLUS L (ITIMES P L)) P (IDIFFERENCE N L))) (if (FGREATERP NRMXL 0.0) then (if (FLESSP (\FLOATAREF FACTORMATRIX L L) 0.0) then (SETQ NRMXL (FMINUS NRMXL))) (BLAS.SCAL (FQUOTIENT 1.0 NRMXL) FACTORMATRIX (IPLUS L (ITIMES P L)) P (IDIFFERENCE N L)) (\FLOATASET (FPLUS 1.0 (\FLOATAREF FACTORMATRIX L L)) FACTORMATRIX L L) (* ; "apply the transform to the remaining columns") (bind TEMP for J from (ADD1 L) to (SUB1 P) do (SETQ TEMP (FMINUS (FQUOTIENT (BLAS.DOTPROD FACTORMATRIX (IPLUS L (ITIMES P L)) P FACTORMATRIX (IPLUS J (ITIMES P L)) P (IDIFFERENCE N L)) (\FLOATAREF FACTORMATRIX L L)))) (BLAS.AXPY TEMP FACTORMATRIX (IPLUS L (ITIMES P L)) P FACTORMATRIX (IPLUS J (ITIMES P L)) P (IDIFFERENCE N L))) (\FLOATASET (\FLOATAREF FACTORMATRIX L L) QRAUX L) (\FLOATASET (FMINUS NRMXL) FACTORMATRIX L L))) FACTORMATRIX))
)
(QROLS
(LAMBDA (QRMATRIX QRAUX Y QTY B RSD YHAT) (* jop%: "23-Jun-86 13:45") (* * Lifted from LINPACK algorithm SQRSL) (BLAS.CHECKARRAY QRMATRIX) (BLAS.CHECKARRAY QRAUX) (BLAS.CHECKARRAY Y) (LET ((N (CL:ARRAY-DIMENSION QRMATRIX 0)) (P (CL:ARRAY-DIMENSION QRMATRIX 1))) (* * Arg Checks) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK QRAUX)) (EQ P (CL:ARRAY-TOTAL-SIZE QRAUX)))) then (HELP "QRAUX not of size P" QRAUX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK Y)) (EQ N (CL:ARRAY-TOTAL-SIZE Y)))) then (HELP "Y not of size N" Y)) (if (NULL QTY) then (SETQ QTY (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY QTY) (if (NOT (EQ N (CL:ARRAY-TOTAL-SIZE QTY))) then (HELP "QTY not of size N" QTY))) (if (NULL B) then (SETQ B (CL:MAKE-ARRAY P (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY B) (if (NOT (EQ P (CL:ARRAY-TOTAL-SIZE B))) then (HELP "B not of size P" B))) (if RSD then (BLAS.CHECKARRAY RSD) (if (NOT (EQ N (CL:ARRAY-TOTAL-SIZE RSD))) then (HELP "RSD not of size N" RSD))) (if YHAT then (BLAS.CHECKARRAY YHAT) (if (NOT (EQ N (CL:ARRAY-TOTAL-SIZE YHAT))) then (HELP "XB not of size N" YHAT))) (* Compute TRANS (Q) * Y) (QRQTY QRMATRIX QRAUX Y QTY) (* * Compute B) (* Set up computation of B) (BLAS.ARRAYBLT QTY 0 1 B 0 1 P) (for J from (SUB1 P) to 0 by -1 do (if (UFEQP (\FLOATAREF QRMATRIX J J) 0.0) then (HELP "Singular Matrix" QRMATRIX)) (\FLOATASET (FQUOTIENT (\FLOATAREF B J) (\FLOATAREF QRMATRIX J J)) B J) (if (NOT (EQ J 0)) then (BLAS.AXPY (FMINUS (\FLOATAREF B J)) QRMATRIX J P B 0 1 J))) (* * Compute RSD) (if RSD then (* Set up computation of RSD) (if (ILESSP P N) then (BLAS.ARRAYBLT QTY P 1 RSD P 1)) (BLAS.ARRAYFILL 0.0 RSD 0 1 P) (bind TEMP for J from (SUB1 (IMIN P (SUB1 N))) to 0 by -1 do (if (NOT (UFEQP (\FLOATAREF QRAUX J) 0.0)) then (SETQ TEMP (\FLOATAREF QRMATRIX J J)) (\FLOATASET (\FLOATAREF QRAUX J) QRMATRIX J J) (BLAS.AXPY (FMINUS (FQUOTIENT (BLAS.DOTPROD QRMATRIX (IPLUS J (ITIMES P J)) P RSD J 1 (IDIFFERENCE N J)) (\FLOATAREF QRMATRIX J J))) QRMATRIX (IPLUS J (ITIMES P J)) P RSD J 1 (IDIFFERENCE N J)) (\FLOATASET TEMP QRMATRIX J J))) (* * Compute YHAT) (if YHAT then (* Set up computation of YHAT) (BLAS.ARRAYBLT QTY 0 1 YHAT 0 1 P) (BLAS.ARRAYFILL 0.0 YHAT P 1) (bind TEMP for J from (SUB1 (IMIN P (SUB1 N))) to 0 by -1 do (if (NOT (UFEQP (\FLOATAREF QRAUX J) 0.0)) then (SETQ TEMP (\FLOATAREF QRMATRIX J J)) (\FLOATASET (\FLOATAREF QRAUX J) QRMATRIX J J) (BLAS.AXPY (FMINUS (FQUOTIENT (BLAS.DOTPROD QRMATRIX (IPLUS J (ITIMES P J)) P YHAT J 1 (IDIFFERENCE N J)) (\FLOATAREF QRMATRIX J J))) QRMATRIX (IPLUS J (ITIMES P J)) P YHAT J 1 (IDIFFERENCE N J)) (\FLOATASET TEMP QRMATRIX J J))))) B))
)
(QRQTY
(LAMBDA (QRMATRIX QRAUX Y PRODUCT) (* ; "Edited 4-Dec-88 18:40 by jop") (* ;;; "COMPUTE (TRANS Q) * Y given a QR factorization described by QRMATRIX and QRAUX where Y is an N vector") (* ;;; "Lifted from LINPACK algorithm SQRSL") (BLAS.CHECKARRAY QRMATRIX) (BLAS.CHECKARRAY QRAUX) (BLAS.CHECKARRAY Y) (LET ((N (CL:ARRAY-DIMENSION QRMATRIX 0)) (P (CL:ARRAY-DIMENSION QRMATRIX 1))) (* ;;; "Arg Checks") (if (NOT (AND (EQ 1 (CL:ARRAY-RANK QRAUX)) (EQ P (CL:ARRAY-TOTAL-SIZE QRAUX)))) then (HELP "QRAUX not of size P" QRAUX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK Y)) (EQ N (CL:ARRAY-TOTAL-SIZE Y)))) then (HELP "Y not of size N" Y)) (if (NULL PRODUCT) then (SETQ PRODUCT (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY PRODUCT) (if (NOT (EQ N (CL:ARRAY-TOTAL-SIZE PRODUCT))) then (HELP "PRODUCT not of size N" PRODUCT))) (BLAS.ARRAYBLT Y 0 1 PRODUCT 0 1 N) (bind TEMP for J from 0 to (SUB1 (IMIN P N)) do (if (NOT (UFEQP (\FLOATAREF QRAUX J) 0.0)) then (SETQ TEMP (\FLOATAREF QRMATRIX J J)) (\FLOATASET (\FLOATAREF QRAUX J) QRMATRIX J J) (BLAS.AXPY (FMINUS (FQUOTIENT (BLAS.DOTPROD QRMATRIX (IPLUS J (ITIMES P J)) P PRODUCT J 1 (IDIFFERENCE N J)) (\FLOATAREF QRMATRIX J J))) QRMATRIX (IPLUS J (ITIMES P J)) P PRODUCT J 1 (IDIFFERENCE N J)) (\FLOATASET TEMP QRMATRIX J J))) PRODUCT))
)
(QRQY
(LAMBDA (QRMATRIX QRAUX Y PRODUCT) (* jop%: "23-Jun-86 13:46") (* * COMPUTE QX given a QR factorization described by QRMATRIX and QRAUX where Y is an N vector) (* * Lifted from LINPACK algorithm SQRSL) (BLAS.CHECKARRAY QRMATRIX) (BLAS.CHECKARRAY QRAUX) (BLAS.CHECKARRAY Y) (LET ((N (CL:ARRAY-DIMENSION QRMATRIX 0)) (P (CL:ARRAY-DIMENSION QRMATRIX 1))) (* * Arg Checks) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK QRAUX)) (EQ P (CL:ARRAY-TOTAL-SIZE QRAUX)))) then (HELP "QRAUX not of size P" QRAUX)) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK Y)) (EQ N (CL:ARRAY-TOTAL-SIZE Y)))) then (HELP "Y not of size N" Y)) (if (NULL PRODUCT) then (SETQ PRODUCT (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY PRODUCT) (if (NOT (EQ N (CL:ARRAY-TOTAL-SIZE PRODUCT))) then (HELP "PRODUCT not of size N" PRODUCT))) (BLAS.ARRAYBLT Y 0 1 PRODUCT 0 1 N) (bind TEMP for J from (SUB1 (IMIN P (SUB1 N))) to 0 by -1 do (if (NOT (UFEQP (\FLOATAREF QRAUX J) 0.0)) then (SETQ TEMP (\FLOATAREF QRMATRIX J J)) (\FLOATASET (\FLOATAREF QRAUX J) QRMATRIX J J) (BLAS.AXPY (FMINUS (FQUOTIENT (BLAS.DOTPROD QRMATRIX (IPLUS J (ITIMES P J)) P PRODUCT J 1 (IDIFFERENCE N J)) (\FLOATAREF QRMATRIX J J))) QRMATRIX (IPLUS J (ITIMES P J)) P PRODUCT J 1 (IDIFFERENCE N J)) (\FLOATASET TEMP QRMATRIX J J))) PRODUCT))
)
(QRSOLV
(LAMBDA (QRMATRIX QRAUX BVECTOR SOLUTION) (* jop%: "27-May-86 20:38") (* * Solves to system A x = b. BVECTOR should to the RHS of the system. Returns SOLUTION) (RSOLV QRMATRIX (QRQTY QRMATRIX QRAUX BVECTOR SOLUTION) SOLUTION))
)
(MREGRESS
(LAMBDA (Y X B RSD YHAT) (* ; "Edited 4-Dec-88 18:36 by jop") (* ;;; "MREGRESS calculates the least squares (multiple) regression of Y on X. An N vector Y.") (LET* ((QRAUX (CL:MAKE-ARRAY (CL:ARRAY-DIMENSION X 1) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) (QRMATRIX (QRFACTOR X QRAUX))) (QROLS QRMATRIX QRAUX Y NIL B RSD YHAT)))
)
(RSOLV
(LAMBDA (RMATRIX CVECTOR BVECTOR) (* jop%: "23-Jun-86 13:46") (* * Calcluate the solution vector BVECTOR for the system of linear equations R*B=C, where RMATRIX is upper triangular M X N with non-zero diagonal elements. BVECTOR and CVECTOR must be of size N. Always returns BVECTOR) (BLAS.CHECKARRAY RMATRIX) (BLAS.CHECKARRAY CVECTOR) (PROG ((M (CL:ARRAY-DIMENSION RMATRIX 0)) (N (CL:ARRAY-DIMENSION RMATRIX 1))) (* * Arg Checks) (if (ILESSP M N) then (HELP "Order of system less than" N)) (if (NOT (EQ (CL:ARRAY-TOTAL-SIZE CVECTOR) N)) then (HELP "CVECTOR not of size" N)) (if (NULL BVECTOR) then (SETQ BVECTOR (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY BVECTOR) (if (NOT (EQ (CL:ARRAY-TOTAL-SIZE BVECTOR) N)) then (HELP "BVECTOR not of size" N))) (* Check for zero diagonal elements) (if (for I from 0 to (SUB1 N) thereis (UFEQP 0.0 (\FLOATAREF RMATRIX I I))) then (HELP "RMATRIX has a zero diagonal element")) (* * Solution by backsubstitution.) (BLAS.ARRAYBLT CVECTOR 0 1 BVECTOR 0 1 N) (LET ((INDEXLIMIT (SUB1 N))) (* Compute the last value) (\FLOATASET (FQUOTIENT (\FLOATAREF BVECTOR INDEXLIMIT) (\FLOATAREF RMATRIX INDEXLIMIT INDEXLIMIT)) BVECTOR INDEXLIMIT) (bind J JLESS1 for JJ from 1 to INDEXLIMIT do (SETQ J (IDIFFERENCE N JJ)) (SETQ JLESS1 (SUB1 J)) (BLAS.AXPY (FMINUS (\FLOATAREF BVECTOR J)) RMATRIX J N BVECTOR 0 1 J) (\FLOATASET (FQUOTIENT (\FLOATAREF BVECTOR JLESS1) (\FLOATAREF RMATRIX JLESS1 JLESS1)) BVECTOR JLESS1))) (RETURN BVECTOR)))
)
(MSOLVE
(LAMBDA (MATRIX CVECTOR SOLUTION) (* ; "Edited 4-Dec-88 18:32 by jop") (* ;;; "Solves to system A x = b. CVECTOR should to the RHS of the system. Returns SOLUTION") (LET ((PIVOTVECTOR (CL:MAKE-ARRAY (CL:ARRAY-DIMENSION MATRIX 0)))) (LUSOLV (LUFACTOR MATRIX PIVOTVECTOR) PIVOTVECTOR CVECTOR SOLUTION)))
)
(SVDFACTOR
(LAMBDA (XMATRIX SVECTOR UMATRIX VMATRIX) (* ; "Edited 5-Dec-88 19:02 by jop") (* ;;; "Singular-value decomposition by means of orthogonalization by plane rotations. Taken from Nash and Shlien: 'Partial svd algorithms.' On entry X 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.") (BLAS.CHECKARRAY XMATRIX) (LET ((M (CL:ARRAY-DIMENSION XMATRIX 0)) (N (CL:ARRAY-DIMENSION XMATRIX 1))) (* ;; "Args checks") (if (NOT (EQ 2 (CL:ARRAY-RANK XMATRIX))) then (HELP "XMATRIX not a matrix" XMATRIX)) (if (NULL SVECTOR) then (SETQ SVECTOR (CL:MAKE-ARRAY N (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY SVECTOR) (if (NOT (AND (EQ 1 (CL:ARRAY-RANK SVECTOR)) (EQ N (CL:ARRAY-TOTAL-SIZE SVECTOR)))) then (HELP "Illegal SVECTOR" SVECTOR))) (if (NULL UMATRIX) then (SETQ UMATRIX (CL:MAKE-ARRAY (LIST M N) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY UMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS UMATRIX) (CL:ARRAY-DIMENSIONS XMATRIX))) then (HELP "Illegal UMATRIX" UMATRIX))) (if (NULL VMATRIX) then (SETQ VMATRIX (CL:MAKE-ARRAY (LIST N N) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) else (BLAS.CHECKARRAY VMATRIX) (if (NOT (EQUAL (CL:ARRAY-DIMENSIONS VMATRIX) (LIST N N))) then (HELP "Illegal VMATRIX" VMATRIX))) (* ;; "Copy XMATRIX to UMATRIX") (BLAS.ARRAYBLT XMATRIX NIL NIL UMATRIX) (* ;; "Initialize VMATRIX to identity matrix.") (BLAS.ARRAYFILL 0.0 VMATRIX) (for I from 0 to (SUB1 N) do (\FLOATASET 1.0 VMATRIX I I)) (* ;; "Start the computation") (LET ((NT N)) (* ;; "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 (EPS ← 1.0E-6) (SLIMIT ← (IMAX (IQUOTIENT N 4) 6)) (SCOUNT ← 0) RCOUNT eachtime (SETQ RCOUNT (IQUOTIENT (ITIMES NT (SUB1 NT)) 2)) (SETQ SCOUNT (ADD1 SCOUNT)) repeatwhile (IGREATERP RCOUNT 0) do (if (IGREATERP SCOUNT SLIMIT) then (HELP "Number of sweeps exceeds sweep limit." SCOUNT)) (for J from 0 to (IDIFFERENCE NT 2) do (bind P Q R C S V for K from (ADD1 J) to (SUB1 NT) 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)) (\FLOATASET Q SVECTOR J) (\FLOATASET R SVECTOR K) (if (FLESSP Q R) then (SETQ Q (FDIFFERENCE (FQUOTIENT Q R) 1.0)) (SETQ P (FQUOTIENT P R)) (SETQ V (SQRT (SETQ V (FPLUS (FTIMES 4.0 P P) (FTIMES Q Q))))) (SETQ S (SQRT (FTIMES 0.5 (FDIFFERENCE 1.0 (FQUOTIENT Q V))))) (if (FLESSP P 0.0) then (SETQ S (FDIFFERENCE 0.0 S))) (SETQ C (FQUOTIENT P (FTIMES V S))) (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) (FTIMES EPS EPS)) (LEQ (FTIMES (FQUOTIENT P Q) (FQUOTIENT P R)) EPS)) then (SETQ RCOUNT (SUB1 RCOUNT)) else (SETQ R (FDIFFERENCE 1.0 (FQUOTIENT R Q))) (SETQ P (FQUOTIENT P Q)) (SETQ V (SQRT (SETQ V (FPLUS (FTIMES 4.0 P P) (FTIMES R R))))) (SETQ C (SQRT (FTIMES 0.5 (FPLUS 1.0 (FQUOTIENT R V))))) (SETQ S (FQUOTIENT P (FTIMES V C))) (* ;; "box before the COLROT calls") (BLAS.ROT C S UMATRIX J N UMATRIX K N M) (BLAS.ROT C S VMATRIX J N VMATRIX K N N)))) (while (AND (IGEQ NT 3) (LEQ (FQUOTIENT (\FLOATAREF SVECTOR (SUB1 NT)) (FPLUS (\FLOATAREF SVECTOR 0) EPS)) EPS)) 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.") (bind Q for J from 0 to (SUB1 N) do (SETQ Q (SQRT (\FLOATAREF SVECTOR J))) (\FLOATASET Q SVECTOR J) (if (ILEQ J NT) then (BLAS.SCAL (FQUOTIENT 1.0 Q) UMATRIX J N M))) SVECTOR)))
)
(SVDTEST
(LAMBDA NIL (* jop%: "30-Jan-86 17:37") (* * comment) (LET ((UU (CL:MAKE-ARRAY (QUOTE (24 19)) (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) (SS (CL:MAKE-ARRAY 19 (QUOTE :ELEMENT-TYPE) (QUOTE FLOAT))) (VV (CL: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))))
)
(\FLOATAREFMACRO
(LAMBDA (ARGS) (* ; "Edited 4-Dec-88 18:22 by jop") (* * macro expander for \FLOATAREF) (if (IGREATERP (LENGTH ARGS) 3) then (HELP "\FLOATAREF takes no more than three args" ARGS)) (PROG ((BARRAY (CAR ARGS)) (BINDICES (CDR ARGS)) INDEXFORM) (if (EQLENGTH BINDICES 1) then (SETQ INDEXFORM (CAR BINDICES)) else (SETQ INDEXFORM (BQUOTE (IPLUS %, (CADR BINDICES) (ITIMES %, (CAR BINDICES) (CL:ARRAY-DIMENSION %, BARRAY 1)))))) (RETURN (BQUOTE (\GETBASEFLOATP (%%ARRAY-BASE %, BARRAY) (LLSH %, INDEXFORM 1))))))
)
(\FLOATASETMACRO
(LAMBDA (ARGS) (* ; "Edited 4-Dec-88 18:21 by jop") (* * macro expander for \FLOATASET) (if (IGREATERP (LENGTH ARGS) 4) then (HELP "\FLOATASET takes no more than four args" ARGS)) (PROG ((BNEWVALUE (CAR ARGS)) (BARRAY (CADR ARGS)) (BINDICES (CDDR ARGS)) INDEXFORM) (if (EQLENGTH BINDICES 1) then (SETQ INDEXFORM (CAR BINDICES)) else (SETQ INDEXFORM (BQUOTE (IPLUS %, (CADR BINDICES) (ITIMES %, (CAR BINDICES) (CL:ARRAY-DIMENSION %, BARRAY 1)))))) (RETURN (BQUOTE (\PUTBASEFLOATP (%%ARRAY-BASE %, BARRAY) (LLSH %, INDEXFORM 1) %, BNEWVALUE)))))
)
)
(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%: EVAL@COMPILE
(PUTPROPS \FLOATAREF MACRO (ARGS (* *) (\FLOATAREFMACRO ARGS)))
(PUTPROPS \FLOATASET MACRO (ARGS (* *) (\FLOATASETMACRO ARGS)))
)
(FILESLOAD BLAS)
(PUTPROPS MATRIXOPS COPYRIGHT ("Xerox Corporation" 1986 1988))
(DECLARE%: DONTCOPY
(FILEMAP (NIL (737 29405 (CHOLESKYFACTOR 747 . 2166) (MTRANSPOSE 2168 . 3120) (MINVERT 3122 . 3412) (
LSOLV 3414 . 4883) (LUFACTOR 4885 . 7439) (LUINVERSE 7441 . 9470) (LUSOLV 9472 . 11455) (MTIMES 11457
. 13586) (QRFACTOR 13588 . 15778) (QROLS 15780 . 18437) (QRQTY 18439 . 19758) (QRQY 19760 . 21058) (
QRSOLV 21060 . 21300) (MREGRESS 21302 . 21642) (RSOLV 21644 . 23147) (MSOLVE 23149 . 23467) (SVDFACTOR
23469 . 27326) (SVDTEST 27328 . 28302) (\FLOATAREFMACRO 28304 . 28834) (\FLOATASETMACRO 28836 . 29403
)))))
STOP