[comp.sys.handhelds] QR and SVD

dominop@en.ecn.purdue.edu (Philippos A. Peleties) (12/02/90)

Hi there!

Here are two HP-28 programs for you Linear Algebra buffs: QR and SVD.

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

QR: 
---

Decompose a real matrix A into a real orthogonal matrix Q and a real
upper triangular matrix R such that A=Q*R.

Input:
------

Level 1: Matrix A (real)
Level 2: Q (real orthogonal)
Level 3: R (real upper triangular)

Checksum [AEB]

<< -> A
   << A SIZE LIST->
DROP -> n m
     << n IDN A 0 0 ->
Q R s c
       << 1 m
         FOR j 1 j +
n DUP2
           IF >
           THEN DROP2
           ELSE
             FOR i 'R
(j,j)' EVAL SQ 'R(i,
j)' EVAL SQ + sqrt DUP
              IF DUP
0 ==
              THEN 1
'c' STO 0 's' STO
DROP2
              ELSE '
R(i,j)' EVAL SWAP /
's' STO 'R(j,j)'
EVAL SWAP / 'c' STO
              END Q
n IDN i i 2 ->LIST c 
PUT i j 2 ->LIST s
NEG PUT j j 2 ->LIST
c PUT j i 2 ->LIST s
PUT TRN DUP TRN R *
'R' STO * 'Q' STO
            NEXT
          END
        NEXT Q R
       >>
     >>
   >>
>>


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

Singular Value Decomposition (SVD)
----------------------------------

A=U*S*V'

Input:
------

Level 2: Matrix A (real)
Level 1: Number of QR itterations


Output:
-------

Level 1: U (real orthogonal)
Level 2: S (singular values)
Level 3: V (real orthogonal)
Level 4: ABS(U*S*V'-A)


Checksum [2FCC]

<< -> A r
   << A SIZE LIST->
DROP A 0 0 -> n m s u
v
    << n IDN 'u' STO
m IDN 'v' STO 1 r
      START u s QR
TRN 's' STO * 'u'
STO v s QR TRN 's'
STO * 'v' STO
      NEXT s
      IF s n m MIN
DUP 2 ->LIST GET 0 <
      THEN n IDN n n
2 ->LIST 1 NEG PUT
DUP u SWAP * 'u' STO
SWAP * 's' STO s
      END u SWAP v u
s v TRN * * A - ABS 
    >>
  >>
>>


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


Enjoy!


Philip Peleties
-- 


I speak for myself, I think for myself, I work for myself ... but I don't want
to play by myself ... so bring your toys and let's share ...