Simple matrix tridiagonalization algorithm in Forth

[Krishna Myneni]

The following module in Forth computes the "similar" matrix, in
tridiagonal form, for a real square matrix, i.e. it tridiagonalizes the
input matrix using a series of similarity transformations. Similarity
transformations preserve many properties of the original matrix, such as
eigenvalues, determinant, trace, etc. -- a nice video explaining this
may be viewed here: https://www.youtube.com/watch?v=L1a6yb3vjIk

The Forth code here implements the procedure given on the web page,

https://en.wikipedia.org/wiki/Householder_transformation

and reduces the original N by N matrix into tridiagonal form through a
series of Householder transformations. The example on the page above is
used within the test code.

Tridiagonalization may be used as a first step in computing the
eigenvalues and eigenvectors of a matrix, so this algorithm may be
useful as the first step, in conjunction with other well-known
algorithms for obtaining the eigenvalues/eigenvectors. The present
implementation does not compute ad save the accumulated similarity
transform, which would be needed for obtaining eigenvectors of the
original matrix. I also don't know what types of numerical instabilities
it may suffer, relative to similar routines from EISPACK or LAPACK, and
I don't know its relative efficiency in comparison with those routines
(e.g. tred2 from EISPACK), but it's a starting point.

The code is designed to work with the Forth Scientific Library, with an
extra matrix multiplication module (mmul.4th), and using the portable
modules framework provided with kForth (kForth-64). It should be easy to
use on modern Forth systems.

--
Krishna Myneni


tridiag.4th
--- start of file ---
\ tridiag.4th
\
\ Tridiagonalize a square matrix using a series of Householder
\ transforms. The input matrix is replaced by its similar
\ tridiagonal matrix.
\
\ Copyright 2022 Krishna Myneni. Permission is granted by the
\ author to use this software for any application provided this
\ copyright notice is preserved.
\
\ Glossary:
\
\ }}TRIDIAG ( N amatrix -- )
\
\ See example(s) of use in test code.
\
\ References:
\ 1. Wikipedia entry on Householder transformation,
\ https://en.wikipedia.org/wiki/Householder_transformation
\
\ 2. Wikipedia entry on matrix similarity,
\ https://en.wikipedia.org/wiki/Matrix_similarity
\
\ Requires:
\ fsl/fsl-util
\ fsl/dynmem
\ fsl/extras/mmul
\

BEGIN-MODULE
BASE @

[undefined] fsquare [IF] : fsquare fdup f* ; [THEN]

FLOAT DMATRIX Ap{{
FLOAT DMATRIX P{{
FLOAT DARRAY v{

fvariable alpha
fvariable 2r
0 value N
0 ptr AA{{

Public:

: }}tridiag ( N a -- )
\ Initialization
to AA{{ to N
N 2- 0 <= IF EXIT THEN

& v{ N }malloc
& P{{ N N }}malloc
& Ap{{ N N }}malloc
malloc-fail? IF 1 throw THEN

N 2- 0 DO
0.0e0 N I 1+ DO AA{{ I J }} f@ fsquare f+ LOOP fsqrt
AA{{ I 1+ I }} f@ f0< IF fnegate THEN fnegate fdup alpha f!
fdup AA{{ I 1+ I }} f@ f- f* 2.0e0 f* fsqrt 2r f!
\ Find reflection plane normal vector, v{
I 1+ 0 DO 0.0e0 v{ I } f! LOOP
AA{{ I 1+ I }} f@ alpha f@ f- 2r f@ f/ v{ I 1+ } f!
N I 2+ DO AA{{ I J }} f@ 2r f@ f/ v{ I } f! LOOP
\ Find similarity transform matrix, P{{, for reflection about plane
N 0 DO
N 0 DO
v{ J } f@ v{ I } f@ f* 2.0e0 f* fnegate
I J = IF 1.0e0 f+ THEN P{{ J I }} f!
LOOP
LOOP
\ Apply similarity transformation to matrix, AA{{
AA{{ P{{ Ap{{ N N N df_mmul
P{{ Ap{{ AA{{ N N N df_mmul
LOOP

\ cleanup
& Ap{{ }}free
& P{{ }}free
& v{ }free
;

BASE !
END-MODULE


TEST-CODE? [IF]
[UNDEFINED] T{ [IF] include ttester.4th [THEN]
BASE @ DECIMAL

1e-14 rel-near f!
1e-14 abs-near f!
set-near

\ Test case: 4 x 4 real, symmetric matrix (see ref. [1]).
4 4 FLOAT MATRIX A{{
4e 1e -2e 2e
1e 2e 0e 1e
-2e 0e 3e -2e
2e 1e -2e -1e
4 4 A{{ }}fput


CR
TESTING }}TRIDIAG
4 A{{ }}tridiag
t{ A{{ 0 0 }} f@ -> 4.0e0 r}t
t{ A{{ 0 1 }} f@ -> -3.0e0 r}t
t{ A{{ 0 2 }} f@ -> 0.0e0 r}t
t{ A{{ 0 3 }} f@ -> 0.0e0 r}t
t{ A{{ 1 0 }} f@ -> A{{ 0 1 }} f@ r}t
t{ A{{ 1 1 }} f@ -> 10.0e0 3.0e0 f/ r}t
t{ A{{ 1 2 }} f@ -> -5.0e0 3.0e0 f/ r}t
t{ A{{ 1 3 }} f@ -> 0.0e0 r}t
t{ A{{ 2 0 }} f@ -> 0.0e0 r}t
t{ A{{ 2 1 }} f@ -> A{{ 1 2 }} f@ r}t
t{ A{{ 2 2 }} f@ -> -33.0e0 25.0e0 f/ r}t
t{ A{{ 2 3 }} f@ -> 68.0e0 75.0e0 f/ r}t
t{ A{{ 3 0 }} f@ -> 0.0e0 r}t
t{ A{{ 3 1 }} f@ -> 0.0e0 r}t
t{ A{{ 3 2 }} f@ -> A{{ 2 3 }} f@ r}t
t{ A{{ 3 3 }} f@ -> 149.0e0 75.0e0 f/ r}t

BASE !
[THEN]

--- end of file ---

test-tridiag.4th
--- start of file ---
\ test-tridiag.4th

include ans-words
include modules
include fsl/fsl-util
include fsl/dynmem
include fsl/extras/mmul
true to test-code?
include fsl/extras/tridiag

[Krishna Myneni]

On 6/2/22 22:05, Krishna Myneni wrote:
> The following module in Forth computes the "similar" matrix, in
> tridiagonal form, for a real square matrix, i.e. it tridiagonalizes the

> input matrix using a series of similarity transformations. ...

> I don't know its relative efficiency in comparison with those routines
> (e.g. tred2 from EISPACK), but it's a starting point.
>

...

I translated the EISPACK subroutine tred1 from Fortran to Forth, and
verified for one matrix case that both Fortran and Forth versions gave
the same output.

TRED1 uses a different set of similarity transforms for
tridiagonalization than used by tridiag.4th; hence, for a given input
double precision symmetric matrix, TRED1 and }}TRIDIAG will give
different tridiagonal matrices. But the output matrices should be
similar, i.e. connected by a similarity transformation.

Differences between TRED1 and }}TRIDIAG include the following:

1. TRED1 accepts a matrix which is assumed to be symmetric, while
}}TRIDIAG can accept a general matrix. In fact TRED1 only uses the lower
triangular part of the input matrix, so it enforces the symmetry.

2. TRED1 takes array arguments for receiving the diagonal and the
subdiagonal elements. }}TRIDIAG replaces the input matrix elements with
the tridiagonlized version. Only two arrays are needed by TRED1 to store
the tridiagonal matrix since the two subdiagonals are the same (symmetry
is preserved).

Although LAPACK replaced EISPACK as the standard for numerical linear
algebra some years ago, I think it is worthwhile to be able to run the
EISPACK routines, particularly in Forth. TRED1 maps to the LAPACK
routine SSYTRD. The EISPACK archive on netlib may be found at

http://www.netlib.org/eispack/

The correspondence and differences between EISPACK routines and their
LAPACK counterparts are described at

https://www.netlib.org/lapack/lug/node147.html

My Forth translation of tred1 may be found at

https://github.com/mynenik/kForth-64/tree/master/forth-src/fsl/extras

Note that FSL arrays, as defined in the fsl-util.4th file, are assumed
by both tred1.4th and tridiag.4th.

--
Krishna Myneni

[Krishna Myneni]

On 6/8/22 09:38, Krishna Myneni wrote:
...
> I translated the EISPACK subroutine tred1 from Fortran to Forth, and
> verified for one matrix case that both Fortran and Forth versions gave
> the same output.
>

...

I completed translation of the EISPACK subroutine IMTQL1 from Fortran to
Forth. The word IMTQL1 computes the eigenvalues of a real symmetric
tridiagonal matrix. The initial Forth version also includes test code
for a 10x10 symmetrized Clement matrix, a tridiagonal matrix with known
integer eigenvalues.

The original Fortran version of the imtql1 EISPACK routine may be found
here:

http://www.netlib.org/eispack/imtql1.f

My Forth version may be found here:

https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/imtql1.4th


The above Forth module, together with the already translated tred1.4th,
and using FSL arrays and matrices, provide a way to compute eigenvalues
of real symmetric matrices in Forth. Thus, finding eigenvalues for one
major category of matrices may now be handled simply in Forth, without
the need for external library interfaces.

Of course, IMTQL1 can only provide eigenvalues. For solving the complete
eigenvalue problem for this category of matrices, we also need the
eigenvectors. For this we need Forth implementations of EISPACK's
tred2() and imtql2() as well.

To run the test code for imtql1.4th under kForth-64/32, the following
statements may be used:

include ans-words
include modules
include fsl/fsl-util

true to test-code?
include fsl/extras/imtql1

--
Krishna Myneni

[Marcel Hendrix]

On Saturday, June 11, 2022 at 6:05:56 PM UTC+2, Krishna Myneni wrote:
> On 6/8/22 09:38, Krishna Myneni wrote:

[..]

> The above Forth module, together with the already translated tred1.4th,
> and using FSL arrays and matrices, provide a way to compute eigenvalues
> of real symmetric matrices in Forth. Thus, finding eigenvalues for one
> major category of matrices may now be handled simply in Forth, without
> the need for external library interfaces.

In 2006, I sent Skip Carter eigenvals.frt (358 lines) for the FSL.
It has been under review for 16 years by now :--)

: .about
CR ." HOUSEHOLDER-REDUCE ( a{{ d{ e{ -- ) must be real and symmetric"
CR ." TRIDIAGONAL-QL-IMPLICIT ( d{ e{ -- bool ) must be real and symmetric"
CR ." TQLI-EIGENVECTORS ( z{{ d{ e{ -- flag ) must be real and symmetric"
CR ." EXTREME-EVS ( s{{ -- ) ( F: -- ev_max ev_min ) any matrix"
CR ." LARGEST-EV ( s{{ -- ) ( F: -- ev ) any matrix"
CR ." SMALLEST-EV ( s{{ -- ) ( F: -- ev ) any matrix"
CR ." SYMMETRIC? ( x{{ -- TRUE=yes ) " ;

-marcel

[Krishna Myneni]

Good to know. Have you posted this code?

There are problems with the FSL review process, and I use the practice
of putting all of my additional FSL-related code in the extras/
subfolder so that they are not taken to be part of the FSL proper. The
review process suffers from

1) lack of many users for the FSL modules

2) lack of reviewers with both competence in Forth and in numerical methods

3) difficulty in performing review to the FSL standards

I have been both a contributor and reviewer. As a reviewer I found the
amount of work in performing a review of a proposed module to be
comparable to or exceeding the effort to review an academic paper. When
I submitted the Faddeeva function (zwofz.4th) for inclusion into the
FSL, no one was willing to do the review. I had submitted the same
function for inclusion into an R package, and it was adopted into R's
NORMT3 package. A somewhat similar experience occurred with my Numerov
integrator package (numerov.4th).

In any case, my goal here is to have a set of Forth modules which map
one to one to the well-tested and well-documented EISPACK routines, for
my own use. The Forth words will be demonstrated to have identical
behavior to their Fortran subroutine counterparts, and thus will be
traceable to EISPACK. The comprehensive book, Matrix Eigensystem
Routines -- Eispack Guide, by B. T. Smith, J. M. Boyle, et al., Springer
2013, appears to be still in print. Other guides also exist on the internet.

--
Krishna

[Marcel Hendrix]

On Sunday, June 12, 2022 at 4:13:36 AM UTC+2, Krishna Myneni wrote:
[..]

> In any case, my goal here is to have a set of Forth modules which map
> one to one to the well-tested and well-documented EISPACK routines, for
> my own use.

Up to last week I thought that those routines could just as well be
accessed with iForth's foreign library API / Flag's EXTERN.
Unfortunately, a few of them work with callbacks, notably root solvers
and optimizers (in my case Nelder and Mead's Amoeba routine, also
used in MATLAB's fmins before release 12).
Although iForth of course supports callbacks (C calls to Forth code),
there is a hairy problem: in such a callback Forth apparently can not
access ( open ) files. In Windows it simply returned immediately and
in the C debugger there was an uncaught exception stating the
process had switched stream IDs.
Instead of struggling with this mess (in C) any further, I decided to
integrate one of the many Forth versions of Amoeba. In the process
I found a very easy way to limit the search to strictly positive vertices,
giving Amoeba the desirable properties of the nonlinear fmincon solver
(inequality constraints). The Forth version used absolute tolerance
for its stopping criterium (a problem when the result is 0 and having
nasty scaling properties). I'm afraid this Amoeba+ is not similar to the
"well-tested and well-documented EISPACK routines" anymore.

> The Forth words will be demonstrated to have identical
> behavior to their Fortran subroutine counterparts, and thus will be
> traceable to EISPACK. The comprehensive book, Matrix Eigensystem
> Routines -- Eispack Guide, by B. T. Smith, J. M. Boyle, et al., Springer
> 2013, appears to be still in print. Other guides also exist on the internet.

Eispack's eigenvalue solvers for *non*-symmetric matrices were too
much effort to code in Forth, even for me :-) Luckily they don't use
callbacks.

-marcel

[Krishna Myneni]

On 6/11/22 11:05, Krishna Myneni wrote:

> https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/imtql1.4th
>
...

> The above Forth module, together with the already translated tred1.4th,
> and using FSL arrays and matrices, provide a way to compute eigenvalues
> of real symmetric matrices in Forth. Thus, finding eigenvalues for one
> major category of matrices may now be handled simply in Forth, without
> the need for external library interfaces.
>

...

An example of using the two modules, tred1 and imtql1, to solve for the
eigenvalues of a real symmetric matrix, using the 4x4 matrix from the
original post, is given below.

--- start of code --

include ans-words ( kForth-only )
include modules
include fsl/fsl-util
include fsl/extras/tred1
include fsl/extras/imtql1

\ 4 x 4 real symmetric matrix example

4 4 FLOAT MATRIX A{{

4.0e0 1.0e0 -2.0e0 2.0e0
1.0e0 2.0e0 0.0e0 1.0e0
-2.0e0 0.0e0 3.0e0 -2.0e0
2.0e0 1.0e0 -2.0e0 -1.0e0
4 4 A{{ }}fput

4 FLOAT ARRAY diag{
4 FLOAT ARRAY subdiag{
4 FLOAT ARRAY subdiag2{

4 4 A{{ diag{ subdiag{ subdiag2{ tred1
4 diag{ subdiag{ imtql1 . cr \ print error code from IMTQL1
4 diag{ }fprint \ print the eigenvalues

--- end of code --

--- output of above using kforth64/kforth32 ---
$ kforth64 tred1-test

/home/krishna/kforth/ans-words.4th

/home/krishna/kforth/modules.4th

/home/krishna/kforth/fsl/fsl-util.4th

FSL-UTIL V1.3c 11 Sep 2021 EFC, KM
/home/krishna/kforth/fsl/extras/tred1.4th

/home/krishna/kforth/fsl/extras/imtql1.4th
0
-2.19752 1.08436 2.26853 6.84462 ok
--- end of output ---

The eigenvalues of the matrix, sorted from smallest to largest, are
shown in the last line of the output. One may use R or matlab to verify
the eigenvalues.

Note that after specifying the input matrix and creating arrays to
receive the diagonal and subdiagonal of the tridiagonalized matrix,
TRED1 is called and then IMTQL1 is called.

--
Krishna Myneni

[Marcel Hendrix]

On Monday, June 13, 2022 at 4:00:40 PM UTC+2, Krishna Myneni wrote:
> On 6/11/22 11:05, Krishna Myneni wrote:
>
> > https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/imtql1.4th
> >
> ...
> > The above Forth module, together with the already translated tred1.4th,
> > and using FSL arrays and matrices, provide a way to compute eigenvalues
> > of real symmetric matrices in Forth. Thus, finding eigenvalues for one
> > major category of matrices may now be handled simply in Forth, without
> > the need for external library interfaces.
> >
> ...
>
> An example of using the two modules, tred1 and imtql1, to solve for the
> eigenvalues of a real symmetric matrix, using the 4x4 matrix from the
> original post, is given below.

[..]

> \ 4 x 4 real symmetric matrix example
> 4 4 FLOAT MATRIX A{{
> 4.0e0 1.0e0 -2.0e0 2.0e0
> 1.0e0 2.0e0 0.0e0 1.0e0
> -2.0e0 0.0e0 3.0e0 -2.0e0
> 2.0e0 1.0e0 -2.0e0 -1.0e0
> 4 4 A{{ }}fput

[..]

> -2.19752 1.08436 2.26853 6.84462 ok

> The eigenvalues of the matrix, sorted from smallest to largest, are
> shown in the last line of the output. One may use R or matlab to verify
> the eigenvalues.

FORTH> a{{ =eigv ok
[2]FORTH> swap }}print ( complex eigenvalues on diagonal )
( 6.844621e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 )
( 0.000000e+0000 0.000000e+0000 ) ( -2.197517e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 )
( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 1.084364e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 )
( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 0.000000e+0000 0.000000e+0000 ) ( 2.268531e+0000 0.000000e+0000 ) ok

[1]FORTH> }}print ( complex eigenvectors in columns )
( 7.180460e-0001 0.000000e+0000 ) ( 1.767052e-0001 0.000000e+0000 ) ( -6.422600e-0001 0.000000e+0000 ) ( -2.017111e-0001 0.000000e+0000 )
( 2.211530e-0001 0.000000e+0000 ) ( 1.781005e-0001 0.000000e+0000 ) ( 5.441878e-0001 0.000000e+0000 ) ( -7.894499e-0001 0.000000e+0000 )
( -5.573514e-0001 0.000000e+0000 ) ( -2.876680e-0001 0.000000e+0000 ) ( -5.202219e-0001 0.000000e+0000 ) ( -5.796342e-0001 0.000000e+0000 )
( 3.533565e-0001 0.000000e+0000 ) ( -9.242849e-0001 0.000000e+0000 ) ( 1.439823e-0001 0.000000e+0000 ) ( -1.028100e-0002 0.000000e+0000 ) ok

The eigenvalues appear to be correct.

-marcel

[Krishna Myneni]

The eigenvalues are real, as they should be for a real symmetric matrix.
If the matrix is real, but not symmetric, then there is the
possibility/likelihood of complex eigenvalues.

--
Krishna

[Marcel Hendrix]

On Monday, June 13, 2022 at 7:43:10 PM UTC+2, Krishna Myneni wrote:
> On 6/13/22 12:30, Marcel Hendrix wrote:

[..]

> >> The eigenvalues of the matrix, sorted from smallest to largest, are

Sorting the eigenvalues (and eigenvectors) is a nice touch.

> >> shown in the last line of the output. One may use R or matlab to verify
> >> the eigenvalues.
> >
> > FORTH> a{{ =eigv ok

[..]

> > [1]FORTH> }}print ( complex eigenvectors in columns )

[..]

> The eigenvalues are real, as they should be for a real symmetric matrix.
> If the matrix is real, but not symmetric, then there is the
> possibility/likelihood of complex eigenvalues.

For electrical engineering problems the matrices are normally non-symmetric (eigenvalues almost always complex).

-marcel

[Krishna Myneni]

On 6/14/22 03:52, Marcel Hendrix wrote:
> On Monday, June 13, 2022 at 7:43:10 PM UTC+2, Krishna Myneni wrote:
>> On 6/13/22 12:30, Marcel Hendrix wrote:
> [..]
>>>> The eigenvalues of the matrix, sorted from smallest to largest, are
>
> Sorting the eigenvalues (and eigenvectors) is a nice touch.
>

The EISPACK routine, IMTQL1, does the sorting as it finds each
additional eigenvalue.

>>>> shown in the last line of the output. One may use R or matlab to verify
>>>> the eigenvalues.
>>>
>>> FORTH> a{{ =eigv ok
> [..]
>>> [1]FORTH> }}print ( complex eigenvectors in columns )
> [..]
>> The eigenvalues are real, as they should be for a real symmetric matrix.
>> If the matrix is real, but not symmetric, then there is the
>> possibility/likelihood of complex eigenvalues.
>
> For electrical engineering problems the matrices are normally non-symmetric (eigenvalues almost always complex).

In physics, I've mostly encountered symmetric real matrices and their
complex counterpart, Hermitian matrices (self-adjoint). Both have real
eigenvalues, with the significance being that the eigenvalues correspond
to physically measurable properties, for example resonant frequencies.

The Harwell-Boeing set of matrices[1], used for testing eigensolvers,
come from actual problems in different disciplines (astrophysics,
quantum chemistry, mechanical structures such as suspension bridge and
aircraft design, power networks, economics, fluid flow, nuclear reactor
modeling, etc), and they include both symmetric and non-symmetric
matrices, often large and sparse.

The EISPACK routines are generally considered useful for matrix order
about N ~ 100. Typical large matrices from applications tend to be
sparse, and specialized eigensolvers for sparse matrices are needed in
those cases. I expect EISPACK routines such as IMTQL1 will work well for
large sparse matrices when the sparse matrix is in symmetric tridiagonal
form. There are also routines present for some relaxed constraints on
tridiagonal matrices.

--
Krishna

1. https://math.nist.gov/MatrixMarket/data/Harwell-Boeing/

[Krishna Myneni]

On 6/11/22 11:05, Krishna Myneni wrote:

> On 6/8/22 09:38, Krishna Myneni wrote:
> ...
>> I translated the EISPACK subroutine tred1 from Fortran to Forth, and
>> verified for one matrix case that both Fortran and Forth versions gave
>> the same output.
>>
> ...
>
> I completed translation of the EISPACK subroutine IMTQL1 from Fortran to
> Forth. The word IMTQL1 computes the eigenvalues of a real symmetric

> tridiagonal matrix. ...

> Of course, IMTQL1 can only provide eigenvalues. For solving the complete
> eigenvalue problem for this category of matrices, we also need the
> eigenvectors. For this we need Forth implementations of EISPACK's
> tred2() and imtql2() as well.
>

...

Initial Forth versions of EISPACK routines, tred2 and imtql2, are now
available. These two routines provide the complete eigensolutions for
real symmetric matrices.

https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/tred2.4th

https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/extras/imtql2.4th

The combination of the two procedures, tred2 and imtql2, provide both
the complete eigenvalues and eigenvectors for a real symmetric matrix.
Example code and its output is provided below to illustrate its use.

Caution: TRED2 and IMTQL2 have only been verified to provide correct
solutions for one test case -- they have not undergone extensive testing.

The eigenvectors agree, to within the expected sign ambiguity, with the
output of Marcel Hendrix's code, for the given example.

--
Krishna Myneni

--- start of code ---
include ans-words ( needed for kForth only )
include modules
include fsl/fsl-util
include fsl/extras/tred2
include fsl/extras/imtql2

\ 4 x 4 real symmetric matrix example

4 4 FLOAT MATRIX A{{

4.0e0 1.0e0 -2.0e0 2.0e0
1.0e0 2.0e0 0.0e0 1.0e0
-2.0e0 0.0e0 3.0e0 -2.0e0
2.0e0 1.0e0 -2.0e0 -1.0e0
4 4 A{{ }}fput

4 FLOAT ARRAY diag{
4 FLOAT ARRAY subdiag{

4 4 FLOAT MATRIX ot{{

4 4 a{{ diag{ subdiag{ ot{{ tred2
4 4 diag{ subdiag{ ot{{ imtql2 . cr

cr .( Eigenvalues: ) cr
4 diag{ }fprint

cr .( Eigenvectors: ) cr
4 4 ot{{ }}fprint \ Each column is an eigenvector.

--- end of code ---


--- start of output ---
$ kforth64
kForth-64 v 0.3.0 (Build: 2022-03-17)
Copyright (c) 1998--2022 Krishna Myneni
Contributions by: dpw gd mu bk abs tn cmb bg dnw imss
Provided under the GNU Affero General Public License, v3.0 or later


Ready!
include test-imtql2

/home/krishna/kforth/ans-words.4th

/home/krishna/kforth/modules.4th

/home/krishna/kforth/fsl/fsl-util.4th

FSL-UTIL V1.3c 11 Sep 2021 EFC, KM

/home/krishna/kforth/fsl/extras/tred2.4th

/home/krishna/kforth/fsl/extras/imtql2.4th


Eigenvalues:
-2.19752 1.08436 2.26853 6.84462
Eigenvectors:
-0.176705 -0.64226 0.201711 0.718046
-0.1781 0.544188 0.78945 0.221153
0.287668 -0.520222 0.579634 -0.557351
0.924285 0.143982 0.010281 0.353356

ok
--- end of output ---

The first eigenvector is column zero of the eigenvector matrix and
corresponds to the first eigenvalue (-2.19752), column 1 to the second
eigenvalue (1.08436), etc.

[Krishna Myneni]

On 6/8/22 09:38, Krishna Myneni wrote:

...
> I translated the EISPACK subroutine tred1 from Fortran to Forth, ...

> My Forth translation of tred1 may be found at
>
> https://github.com/mynenik/kForth-64/tree/master/forth-src/fsl/extras
>

...

Please note that the EISPACK routines translated to Forth so far have
been moved within the kForth repos (kForth-64, kForth-32, and
kForth-Win32). Previous links to the Forth EISPACK files, given in this
thread, will not work.

Within the kForth repo folder structure, the files have been moved from
"forth-src/fsl/extras/" to "forth-src/eispack/". Demo Forth files
illustrating the use of the translated EISPACK routines have been moved
from "forth-src/fsl/demo" to "forth-src/eispack/demo/".

To date, I have translated the following EISPACK routines to Forth:

tred1
tred2
imtql1
imtql2

The two existing demo files are

tred1-ex01.4th
tred2-ex01.4th

A Forth translation of the tridiagonalization routine for complex
Hermitian matrices, htridi, is in preliminary testing right now. I hope
to add additional translations of the EISPACK routines in the future, to
be able to do the full eigensolutions for real symmetric and real
general matrices, and for at least complex Hermitian matrices. That's
the goal anyway, for now. Extensive testing of the Forth translations
will be needed.

--
Krishna Myneni