(FILECREATED "29-May-86 12:41:55" {QV}<PEDERSEN>LISP>OTHERMATRIXOPS.;2
 19199  

      changes to:  (VARS OTHERMATRIXOPSCOMS)
                   (FNS SVDNASH)

      previous date: "27-May-86 22:00:09" {QV}<PEDERSEN>LISP>OTHERMATRIXOPS.;1)


(* Copyright (c) 1986 by Xerox Corporation. All rights reserved.)

(PRETTYCOMPRINT OTHERMATRIXOPSCOMS)

(RPAQQ OTHERMATRIXOPSCOMS
 ((* File created by Coms Manager.)
                           (FNS DIAGXTXINVERSE OLSMGS QRMGS SVDNASH)))



(* File created by Coms Manager.)


(DEFINEQ

(DIAGXTXINVERSE

  [LAMBDA (R D)
                                              (* jop: " 3-Feb-86 15:41")

          
          (* * Return (as the P vector D) the diagonal elements of the matrix

          (R-transpose R) Inverse, where R is P x P upper triangular with non-zero 

          diagonal elements. When the decomposition X=Q*R has been formed, these diagonal 

          elements are just those of (X-transpose X) Inverse

          (mathematically)%. Notice that the vector D must be created by the CALLER!)


    (BLAS.CHECKARRAY R)
    (BLAS.CHECKARRAY D)
    (PROG [(P (ARRAY-DIMENSION R 0))
           (C (MAKE-ARRAY (ARRAY-DIMENSION R 0)
                     (QUOTE :ELEMENT-TYPE)
                     (QUOTE FLOAT]
          
          (* * Check args)


          (if
 (NOT (IEQP P (ARRAY-DIMENSION R 1)))
              then
 (HELP "R not P x P" R))
          (if
 (NOT (IEQP P (ARRAY-DIMENSION D 0)))
              then
 (HELP "D not a P vector" D))
          
          (* * The solution proceeds as a succession of back substitutions against 

          columns of the identity matrix.)


          (BLAS.ARRAYFILL 0.0 D)
          (bind
 FTEMP declare
 (TYPE FLOATP FTEMP) for
 K from
 (SUB1 P) to
 0 by
 -1
             do
 (\FLOATASET (FQUOTIENT 1.0 (\FLOATAREF R K K))
                       C K)
                (SETQ FTEMP (\FLOATAREF C K))
                (\FLOATASET (FPLUS (\FLOATAREF D K)
                                   (FTIMES FTEMP FTEMP))
                       D K)
                (for
 J from
 (SUB1 K) to
 0 by
 -1
                   do
 (SETQ FTEMP 0.0)
                      (for
 I from
 (ADD1 J) to
 K do
 [SETQ FTEMP (FDIFFERENCE FTEMP
                                                                      (FTIMES (\FLOATAREF R J I)
                                                                             (\FLOATAREF C I]
                         finally
 (\FLOATASET (FQUOTIENT FTEMP (\FLOATAREF R J J))
                                        C J))
                      (SETQ FTEMP (\FLOATAREF C J))
                      (\FLOATASET (FPLUS (\FLOATAREF D J)
                                         (FTIMES FTEMP FTEMP))
                             D J)))
          (RETURN NIL])

(OLSMGS

  [LAMBDA (Y X Q R B)
                                        (* jop: "27-May-86 17:08")

          
          (* * OLSMGS calculates the least squares

          (multiple) regression of Y on X. An N vector Y, and a N x P matrix X are 

          required. Q is a N x P matrix containing an orthogonal basis for the column 

          spce of X. R is an upper triangular P x P array.

          B is a P vector of coefficients)


    (PROG ((N (ARRAY-DIMENSION X 0))
           (P (ARRAY-DIMENSION X 1))
           C)
          
          (* * Arg checks)


          (if
 (NOT (IEQP (ARRAY-TOTAL-SIZE Y)
                         N))
              then
 (HELP "Y not an N vector" Y))
          (if
 (NULL Q)
              then
 (SETQ Q (MAKE-ARRAY (ARRAY-DIMENSIONS X)
                                  (QUOTE :ELEMENT-TYPE)
                                  (QUOTE FLOAT)))
            elseif
 (NOT (EQUAL (ARRAY-DIMENSIONS Q)
                               (ARRAY-DIMENSIONS X)))
              then
 (HELP "Q not N x P" Q))
          (if
 (NULL R)
              then
 (SETQ R (MAKE-ARRAY (LIST P P)
                                  (QUOTE :ELEMENT-TYPE)
                                  (QUOTE FLOAT)))
            elseif
 (NOT (EQUAL (ARRAY-DIMENSIONS R)
                               (LIST P P)))
              then
 (HELP "R not P x P" R))
          (if
 (NULL B)
              then
 (SETQ B (MAKE-ARRAY P (QUOTE :ELEMENT-TYPE)
                                  (QUOTE FLOAT)))
            elseif
 (NOT (IEQP (ARRAY-TOTAL-SIZE R)
                              P))
              then
 (HELP "B not a P vector" B))
          (SETQ C (MAKE-ARRAY P (QUOTE :ELEMENT-TYPE)
                         (QUOTE FLOAT)))
          
          (* * Compute Q R decomposition)


          (QRMGS
 X Q R)
          (for
 J from
 0 to
 (SUB1 P) do
 (bind
 (FTEMP _ 0.0) declare
 (TYPE FLOATP FTEMP) for
 I
                                          from
 0 to
 (SUB1 N)
                                          do
 [SETQ FTEMP (FPLUS FTEMP (FTIMES (\FLOATAREF Q I J)
                                                                             (\FLOATAREF Y I]
                                          finally
 (\FLOATASET FTEMP C J)))
          (RSOLV R C B)
          (RETURN B])

(QRMGS

  [LAMBDA (X Q R)
                                            (* jop: " 3-Feb-86 16:04")

          
          (* * This function determines a decomposition X = Q*R for the rectangular 

          CMLArray matrix X, where Q is an orthonormal basis for the column space of X 

          and R is a square upper triangular matrix.)

          
          (* * On entry: X is an n x p CMLARRAY, i.e.

          a matrix, whose QR decomposition is required;

          Q is an n x p array; R is a p x p array)

          
          (* * On return: the value of the function is NIL for a normal return.

          Should this routine be unable to complete the QR decomposition of X because of 

          near singularity of its columns, the function will have a nonNIL value;

          the matrix Q has been OVERWRITTEN by the orthonormal basis Q;

          R contains in its upper triangle the factor R from the decomposition.)

          
          (* * Algorithm notes: The method employed for the decomposition is 

          modified-Gram-Schmidt with reorthogonalization of the columns of Q.

          A test relative to the machine arithmetic precision is made to recognize 

          possible singularity of X. The decision to reorthoganalize is based on the 

          relative decrease in column norm obtained in the Gram-Schmidt process.

          These ideas are due to David M. Gay. For the sake of time and space economy, 

          this function does not perform column pivoting;

          this will exclude successful decomposition for only the most severely 

          ill-conditioned problems. For the sake of time and space economy, no special 

          precautions have been taken in the computation of inner products and vector 

          2-norms. Only those problems where intermediate quantities lie outside the 

          range 10^-19 to 10^19 will notice difficulties, but values of these magnitudes 

          are surely indication of very poor scaling or severe ill-conditioning.)


    (BLAS.CHECKARRAY X)
    (BLAS.CHECKARRAY Q)
    (BLAS.CHECKARRAY R)
    (PROG [(N (ARRAY-DIMENSION X 0))
           (P (ARRAY-DIMENSION X 1))
           (C (MAKE-ARRAY (ARRAY-DIMENSION X 1)
                     (QUOTE :ELEMENT-TYPE)
                     (QUOTE FLOAT]
          
          (* * Check args)


          (if
 (NOT (AND (IEQP (ARRAY-DIMENSION Q 0)
                              N)
                        (IEQP (ARRAY-DIMENSION Q 1)
                              P)))
              then
 (HELP "Q not n x p" Q))
          (if
 (NOT (AND (IEQP (ARRAY-DIMENSION R 0)
                              P)
                        (IEQP (ARRAY-DIMENSION R 1)
                              P)))
              then
 (HELP "R not p x p" R))
          
          (* * Copy X to Q)


          (BLAS.ARRAYBLT X 0 1 Q 0 1)
          
          (* * Initialize diagonal of R to 2-norms of columns of Q)


          (bind
 NORMXJ for
 J from
 0 to
 (SUB1 P) do
 (SETQ NORMXJ (BLAS.NRM2 Q J P N))
                                                   (\FLOATASET NORMXJ R J J)
                                                   (\FLOATASET (if
 (FGREATERP NORMXJ 0.0)
                                                                   then
 0.0
                                                                 else
 1.0)
                                                          C J))
                                                             (* Loop over each column of Q in turn.)

          [bind
 RKK CK V8 V9 declare
 (TYPE FLOATP RKK CK V8) for
 K from
 0 to
 (SUB1 P)
             do
 (SETQ RKK (\FLOATAREF R K K))
                (\FLOATASET 1.0 R K K)
                (SETQ V9 RKK)
                (SETQ CK (\FLOATAREF C K))
                (bind
 NO.ORTHOFLG repeatuntil
 NO.ORTHOFLG
                   do
 [if
 (UFGREATERP CK .75)
                          then
                               (* Weight loss has been substantial, 

                                                             recompute 2-norm.)

                               (SETQ V8 (BLAS.NRM2 Q K P N))
                        else
                                 (* Otherwise, fast update of length.)

                             (SETQ V8 (FTIMES V9 (SQRT (SETQ V8 (FDIFFERENCE 1.0 CK]
                      (\FLOATASET (FTIMES V8 (\FLOATAREF R K K))
                             R K K) 
          
          (* Test if new column length is so much smaller than original length that we 

          must give up and declare singularity. The test relies implicitly on the finite 

          precision of machine addition.)


                      (if
 [NOT (UFLESSP RKK (FPLUS RKK (\FLOATAREF R K K]
                          then
 (RETURN K)
                        else
                                 (* Otherwise, accept column length.)

                             (for
 I from
 0 to
 (SUB1 N)
                                do
 (\FLOATASET (FQUOTIENT (\FLOATAREF Q I K)
                                                      V8)
                                          Q I K)))           (* Column K should be nearly 

                                                             orthogonal to the previous ones, and 

                                                             we check if it can benefit from 

                                                             reorthoganalization.)

                      (if
 (UFLESSP (FQUOTIENT V8 V9)
                                 .25)
                          then
 (SETQ NO.ORTHOFLG NIL)        (* Reorthoganalize with respect 

                                                             preceding columns.)

                               (bind
 DOTXJXK declare
 (TYPE FLOATP DOTXJXK) for
 J from
 0
                                  to
 (SUB1 K)
                                  do
 (SETQ DOTXJXK (BLAS.DOTPROD Q J P Q K P N))
                                     (\FLOATASET (FPLUS (\FLOATAREF R J K)
                                                        (FTIMES DOTXJXK (\FLOATAREF R K K)))
                                            R J K)
                                     (for
 I from
 0 to
 (SUB1 N)
                                        do
 (\FLOATASET (FDIFFERENCE (\FLOATAREF Q I K)
                                                              (FTIMES DOTXJXK (\FLOATAREF Q I J)))
                                                  Q I K)))
                               (SETQ V9 1.0)
                        else
 (SETQ NO.ORTHOFLG T)))          (* Orthogonalize remaining columns 

                                                             with respect to current one.)

                (bind
 DOTXJXK declare
 (TYPE FLOATP DOTXJXK) for
 J from
 (ADD1 K)
                   to
 (SUB1 P) do
 (SETQ DOTXJXK (BLAS.DOTPROD Q J P Q K P N))
                                  (\FLOATASET DOTXJXK R K J)
                                  (for
 I from
 0 to
 (SUB1 N)
                                     do
 (\FLOATASET (FDIFFERENCE (\FLOATAREF Q I J)
                                                           (FTIMES DOTXJXK (\FLOATAREF Q I K)))
                                               Q I J))
                                  (SETQ V9 (\FLOATAREF R J J))
                                  (if
 (UFGREATERP V9 0.0)
                                      then
 (\FLOATASET (FPLUS (\FLOATAREF C J)
                                                              (FTIMES (FQUOTIENT DOTXJXK V9)
                                                                     (FQUOTIENT DOTXJXK V9)))
                                                  C J]
          (RETURN NIL])

(SVDNASH

  [LAMBDA (UMATRIX SVECTOR VMATRIX)
                          (* jop: "31-Jan-86 13:50")

          
          (* * 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))
          
          (* * 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)
                RCOUNT
                (FEPS _ EPS)
                (FEPSSQUARED _ (FTIMES EPS EPS)) declare
 (TYPE FLOATP FEPS FEPSSQUARED)
             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
 (bind
 P Q R C FC S FS V for
 K from
 (IPLUS J 1) to
 (SUB1 NT)
                         declare
 (TYPE FLOATP P Q R FC FS 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)
                            (if
 (UFLESSP 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 FS (SQRT (SETQ FS (FTIMES .5 (FDIFFERENCE 1.0
                                                                               (FQUOTIENT Q V]
                                     (if
 (UFLESSP P 0.0)
                                         then
 (SETQ FS (FDIFFERENCE 0.0 FS)))
                                     (SETQ FC (FQUOTIENT P (FTIMES V FS))) 
                                                             (* box before the COLROT calls)

                                     (SETQ C FC)
                                     (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 (UFLEQ (FTIMES Q R)
                                                FEPSSQUARED)
                                         (UFLEQ (FTIMES (FQUOTIENT P Q)
                                                       (FQUOTIENT P R))
                                                FEPS))
                                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 FC (SQRT (SETQ FC (FTIMES .5 (FPLUS 1.0 (FQUOTIENT R V]
                                   (SETQ FS (FQUOTIENT P (FTIMES V FC))) 
                                                             (* box before the COLROT calls)

                                   (SETQ C FC)
                                   (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]
                (bind
 (TOL _ EPS) declare
 (TYPE FLOATP TOL)
                   while
 (AND (IGEQ NT 3)
                              (UFLEQ (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.)


          (bind
 Q for
 J from
 0 to
 (SUB1 N) do
 (SETQ Q (SQRT (LAREF SVECTOR J)))
                                              (LASET Q SVECTOR J)
                                              (if
 (ILEQ J NT)
                                                  then
 (BLAS.SCAL (FQUOTIENT 1.0 Q)
                                                              UMATRIX J N M])
)
(PUTPROPS OTHERMATRIXOPS COPYRIGHT ("Xerox Corporation" 1986))
(DECLARE: DONTCOPY
  (FILEMAP (NIL (529 19114 (DIAGXTXINVERSE 539 . 2873) (OLSMGS 2875 . 5265) (QRMGS 5267 . 13239) (
SVDNASH 13241 . 19112)))))
STOP