138 views

Skip to first unread message

Jun 2, 2022, 11:05:41 PMJun 2

to

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

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

Jun 8, 2022, 10:38:07 AMJun 8

to

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. ...
> The following module in Forth computes the "similar" matrix, in

> tridiagonal form, for a real square matrix, i.e. it tridiagonalizes the

> I don't know its relative efficiency in comparison with those routines

> (e.g. tred2 from EISPACK), but it's a starting point.

>

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

Jun 11, 2022, 12:05:56 PMJun 11

to

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

include fsl/extras/imtql1

--

Krishna Myneni

Jun 11, 2022, 4:44:44 PMJun 11

to

On Saturday, June 11, 2022 at 6:05:56 PM UTC+2, Krishna Myneni wrote:

> On 6/8/22 09:38, 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.
> 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.

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

Jun 11, 2022, 10:13:36 PMJun 11

to

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

Jun 12, 2022, 7:49:48 AMJun 12

to

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

[..]

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

[..]

> 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
> one to one to the well-tested and well-documented EISPACK routines, for

> my own use.

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.

much effort to code in Forth, even for me :-) Luckily they don't use

callbacks.

-marcel

Jun 13, 2022, 10:00:40 AMJun 13

to

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

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

>

...

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

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

>

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

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

Jun 13, 2022, 1:30:53 PMJun 13

to

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.

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

[..]
> 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
> shown in the last line of the output. One may use R or matlab to verify

> the eigenvalues.

[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

Jun 13, 2022, 1:43:10 PMJun 13

to

If the matrix is real, but not symmetric, then there is the

possibility/likelihood of complex eigenvalues.

--

Krishna

Jun 14, 2022, 4:52:03 AMJun 14

to

On Monday, June 13, 2022 at 7:43:10 PM UTC+2, Krishna Myneni wrote:

> On 6/13/22 12:30, Marcel Hendrix 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).
> If the matrix is real, but not symmetric, then there is the

> possibility/likelihood of complex eigenvalues.

-marcel

Jun 14, 2022, 9:45:50 AMJun 14

to

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

>

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

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/

Jun 14, 2022, 5:56:32 PMJun 14

to

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

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

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

>

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

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 FLOAT ARRAY subdiag{

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/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
--- end of output ---

corresponds to the first eigenvalue (-2.19752), column 1 to the second

eigenvalue (1.08436), etc.

Jun 16, 2022, 10:30:25 AMJun 16

to

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

>

...
>

> 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

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu