rodatz@cat.de (Christoph Rodatz) (02/24/91)
Morgan R LaVake (akcs.lavake@hpcvbbs.UUCP) schrieb am 23. 2.: > Several days ago a set of programs to compute Eigen values and vectors > was posted. I can't seem to find it now. Can someone help me out? > > E-mail to c-lmr@math.utah.edu > Thanks, > > Morgan I'd also be interested. So I'd appreciate it if someone could help me out just as well. Thanks, Chris --- Christoph Rodatz (rodatz@cat.de) C.A.T. Kommunikations-System, Frankfurt, Germany
dominop@en.ecn.purdue.edu (Philippos A. Peleties) (02/25/91)
Well, here are my versions of eigenvalue routines.
The first one is through the Leverierre algorithm. Not very good
(roots of polynomials are sensitive to coefficient perturbations) but
quite adequate for quick checks.
The second routine converts a matrix to a schur form. Use 0 or 1 to
indicate whether you want shift or not. Shift results in quadratic
convergence rate. Sometimes you might not want to use shift though.
As far as eigenvectors are concerned, that's not so easy. I don't
have any eigenvector routines ( null(A-lambda*I) is not the way to
go). Does anybody have any intelligent algorithms he/she wants to
provide?
Have fun.
----------------------------------------------------------------------
Name: EIGEN
----
Status: In hp28s.
------
Author: Philip Peleties.
------
Function: Compute the eigenvalues of a matrix via the
-------- Leverierre algorithm.
Input: A Matrix.
-----
Output: List with eigenvalues.
------
Calls: TRACE, BAIR
-----
Checksum: DFBF
--------
<< -> A
<< A SIZE LIST->
DROP IDN -> n N
<< 1 1 n
FOR i N A *
TRACE NEG i / DUP n
IDN * N A * + 'N'
STO
NEXT 1 n +
->LIST BAIR
>>
>>
>>
------------------------------------------------------------------
Name: TRACE
----
Status: In hp28s.
------
Author: Philip Peleties.
------
Function: Computes the trace of a square matrix.
--------
Input: A Square matrix.
-----
Output: Trace of A.
------
Calls: None.
-----
Checksum: A803
--------
<< -> A
<< 0 A SIZE LIST->
DROP2 1 SWAP
FOR i A i i 2
->LIST GET +
NEXT
>>
--------------------------------------------------------------------
Name: BAIR
----
Status: In hp28s.
------
Author: Mark Wicks.
------
Function: Computes the roots, both real and complex, of a
-------- polynomial.
Input: List containing coefficients of polynomial.
-----
Output: List containing roots of the polynomial.
------
Calls: NR.
-----
Checksum: 9AF7.
--------
<< { } SWAP 0 0 0 0 0
0 0 0 0 0 0
[ 0 0 ]
-> a b r s p q b1 b2
bp1 bp2 bq1 bq2 opq
<<
IF a SIZE DUP 2
/ IP 2 * ==
THEN a NR SWAP
'a' STO +
END
WHILE a SIZE 3 >
REPEAT 'a' 2 GET
'a' 1 GET / 'p' STO
'a' 3 GET 'a' 1 GET
/ 'q' STO
IF p 0 ==
THEN .1 'p'
STO
END
IF q 0 ==
THEN .2 'q'
STO
END
DO 'a' 1 GET
'b2' STO 'a' 2 GET p
b2 * - 'b1' STO 0
'bp2' STO 0 'bq2'
STO b2 NEG 'bp1' STO
0 'bq1' STO b2 b1 2
->LIST 'b' STO 3 a
SIZE 2 - DUP2
IF <=
THEN
FOR i 'a'
i GET p b1 * - q b2
* - b1 NEG p bp1 * -
q bp2 * - bp1 'bp2'
STO 'bp1' STO b2 NEG
p bq1 * - q bq2 * -
bq1 'bq2' STO 'bq1'
STO b1 'b2' STO 'b1'
STO b b1 + 'b' STO
NEXT
ELSE DROP2
END a DUP
SIZE 1 - GET p b1 *
- q b2 * - 'r' STO a
DUP SIZE GET q b1 *
- 's' STO p q 2
->ARRY r s 2 ->ARRY b1
NEG p bp1 * - q bp2
* - b2 NEG p bq1 * -
q bq2 * - q NEG bp1
* b1 NEG q bq1 * - {
2 2 } ->ARRY / -
ARRY-> DROP 'q' STO
'p' STO
UNTIL p q 2
->ARRY DUP opq - ABS
SWAP DUP 'opq' STO
ABS .0000000001 * <
END p NEG 2 /
DUP SQ q - $sqrt$ DUP2 +
3 ROLLD - 2 ->LIST +
b 'a' STO
END 'a' 2 GET
NEG 'a' 1 GET / 2 /
DUP SQ 'a' 3 GET 'a'
1 GET / - $sqrt$ DUP2 + 3
ROLLD - 2 ->LIST +
>>
>>
--------------------------------------------------------------------
Name: NR
----
Status: In hp28s.
------
Author: Mark Wicks.
------
Function: Computes a single real root of a polynomial and returns
-------- the deflated polynomial and the root.
Input: List containing coefficients of polynomial.
-----
Output: A list containing the coefficients of the deflated
------ polynomial.
A single real root of the polynomial.
Calls: None.
-----
Checksum: B7C2
--------
<< 0 0 0 0 0 -> a c r
p b db
<< 'a' 2 GET 'a' 1
GET / 'p' STO
DO 0 'db' STO a
1 GET 'b' STO b 1
->LIST 2 a SIZE 1 -
FOR i 'a' i
GET p b * - db p * b
+ NEG 'db' STO DUP
'b' STO +
NEXT 'c' STO
'a' a SIZE GET p b *
- b p db * + NEG / p
SWAP DUP2 - 'p' STO
UNTIL ABS SWAP
ABS .0000000001 * <
END c p NEG
>>
>>
----------------------------------------------------------------------
Name: SCHUR
----
Status: In hp28s.
------
Author: Philip Peleties.
------
Function: Compute the Schur form of a matrix.
--------
Input: A Matrix.
-----
r Number of QR iterations.
s Shift indicator (0 if no shift, 1 if shift).
Output: Q Orthogonal matrix.
------
T Schur form (A=Q*T*Q')
number ABS(Q*T*Q'-A)
Calls: QR
-----
Checksum: 5E5E
--------
<< -> A r s
<< A SIZE LIST->
DROP IDN DUP
IF s 0 ==
THEN 0 *
END A -> n Q IDNT
T
<< 1 r
START Q T 'T(n
,n)' EVAL IDNT * -
QR DROP DUP DUP TRN
T ROT * * 'T' STO *
'Q' STO
NEXT Q T Q T Q
TRN * * A - ABS
>>
>>
>>
----------------------------------------------------------------
Name: QR
----
Status: In hp28s.
------
Author: Philip Peleties.
------
Function: Compute the QR factorization of a matrix.
--------
Input: A Matrix.
-----
Output: Q Orthogonal matrix.
------ R Upper triangular matrix (A=Q*R).
Calls: None.
-----
Checksum: 7937
--------
<< -> A
<< A SIZE LIST->
DROP -> n m
<< n IDN DUP A
TYPE
IF 4 ==
THEN (1,0) *
END A 0 0 -> Q
IDNT R s c
<< 1 m
FOR j 1 j +
n DUP2
IF >
THEN DROP2
ELSE
FOR i 'R
(j,j)' EVAL 'R(i,j)'
EVAL DUP ABS
IF 0
==
THEN 1
'c' STO 0 's' STO
DROP2
ELSE
DUP2 ABS SWAP ABS
IF >
THEN ABS / DUP ABS
SQ 1 + $sqrt$ INV DUP 'R(
i,j)' EVAL DUP ABS /
* 's' STO * 'c' STO
ELSE SWAP ABS / DUP
ABS SQ 1 + $sqrt$ INV DUP
'R(j,j)' EVAL DUP
ABS / * 'c' STO *
's' STO
END
END Q
IDNT i i 2 ->LIST c
PUT i j 2 ->LIST s
NEG PUT j j 2 ->LIST
c CONJ PUT j i 2
->LIST s CONJ PUT TRN
DUP TRN R * 'R' STO
* 'Q' STO
NEXT
END
NEXT Q R
IF TYPE 4 ==
THEN IDNT n
m MIN DUP 2 ->LIST
DUP R SWAP GET DUP
ABS
IF
.000000000001 <=
THEN 3
DROPN
ELSE DUP
ABS / PUT DUP 3
ROLLD * SWAP TRN R *
'R' STO
END
END R
>>
>>
>>
>>
---------------------------------------------------------------------------------
--
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 ...