Article 4537 of comp.sys.handhelds:
Newsgroups: comp.sys.handhelds
Path: en.ecn.purdue.edu!dominop
From: dominop@en.ecn.purdue.edu (Philippos A. Peleties)
Subject: Re: Eigen-values/vectors request
Message-ID: <1991Feb25.145302.9837@en.ecn.purdue.edu>
Keywords: Eigenvalues, Schur form
Organization: Purdue University Engineering Computer Network
References: <1991Feb24.115734.14198@cat.de>
Date: Mon, 25 Feb 91 14:53:02 GMT

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 ...


