Graham and J-LA-packers generally
FWIW, I've had a go at scripting dggev.ijs and zggev.ijs,
although you don't need them at the moment.
When I got round to testing them, I discovered that they
don't have entry points; ie routines dggev_ and zggev_
are missing from jlapack.dll. So, pro tem, these scripts
call the "deprecated" routines dgegv_ and zgegv_ .
I've modelled them on dgeev.ijs and zgeev.ijs currently
supplied with the LAPACK add-on. Just as in those paradigms,
there's no attempt to choose between none, one or both of
the sets of eigenvectors; both left and right sets of
generalised eigenvectors are returned, together with
"alpha" and "beta", the generalised eigenvalues. The
allotment of work-space is possibly overgenerous.
There are a couple of testing verbs, but on very limited
data, found at http://www.nag.co.uk/lapack-ex/node136.html
The results shown there (for right eigenvectors only)
appear consistent with those from the J scripts apart from
sign and rounding.
I don't know how to incorporate these in the library of
LAPACK scripts, and in any case, they should be properly
tested with suitable data by experts in linear algebra.
I'm attaching dggev.ijs and zggev.ijs as separate files,
though they might be displayed in line by some mail
readers.
Apologies in advance for any errors...
Mike
Graham Parkhouse wrote:
Thank you Raul and Mike.
I wrote:
Eigenvalues again. For structural dynamics it is the generalised
Eigenvalue problem that is needed, the solution to
(A - w B) r = 0
where A and B are real square matrices, w the eigenvalue and r is the
eigenvector.
Raul asked:
For your case are A and B known, and w and r unknown? Or something
different?
Yes, A and B are known, and w and r are unknown.
Mike suggested:
There's a dgegv.f in ~\addons\math\lapack . Its preamble comments
say it's deprecated in favour of dggev... so are you looking for an
implementation of dggev.ijs?
Yes, I read that too. Though dgegv does what I want, maybe dggev is
more accurate.
I have subsequently found that at the moment I don't need either of
them, because for the typical problems I am attempting B is always a
unit matrix, so all I need is dgeev, which I've got.
Thanks for your help. Sorry for the false alarm.
Graham
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
NB. dggev_ is not installed in jlapack.dll
NB. so will use dgegv_ instead, though "deprecated"
NB. dggev returns generalized eigenvalues, and generalized eigenvectors
NB. for a pair of N-by-N real nonsymmetric matrices (A,B)
NB. This script is based on
NB. a- fortran source listing for this problem: dggev.f
NB. b- fortran source listing for eigenvalues for square matrix: dgeev.f
NB. c- j script ~\addons\lapack\math\lapack\dgeev.ijs
NB. - as a model of how to call dgeev_
NB. MJL Day 10/07
coclass 'jlapack'
NB. =========================================================
NB.*dggev v generalised eigenvalues and eigenvectors for real matrices A B
NB.
NB. form: dggev mata;matb
NB.
NB. returns:
NB. left eigenvectors; alpha; beta; right eigenvectors
NB. (all generalised)
NB. The generalized eigenvalue is represented as the pair ( alpha,beta)
NB. Eigenvalue V = alpha%beta
NB. there is a reasonable interpretation for beta = 0, or both zero
NB.
NB. if:
NB. 'L ALFA BETA R' =. dggev A;B
NB. then
NB. V is ALFA % BETA NB. careful....
NB. A mp R is V *"1 B mp R
NB. (+|:L)mp A is V * (+|:L) mp B
dggev=: 3 : 0
'ma mb'=. y
NB. validation -----------------------------------
NB. vsquare each ma;mb
if. (iscomplex ma) +. iscomplex mb do.
need 'zggev'
zggev y
return.
end.
NB. define variables needed for DLL call ---------
s=. $ma
n=. lda=. ldvl=. ldvr=. #ma
a=. dzero + |: ma
b=. dzero + |: mb
ldb =. #mb
alphai =. alphar =. beta =. n$dzero
jobvl=. jobvr=. 'V'
vl=. vr=. a * dzero
lwork=. 20*n
work=. lwork$dzero
info=. izero
NB.
JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHAR,ALPHAI,BETA,VL,LDVL,VR,LDVR,WORK,LWORK,INFO
arg=.'jobvl;jobvr;n;a;lda;b;ldb;alphar;alphai;beta;vl;ldvl;vr;ldvr;work;lwork;info'
NB. call DLL, reassign variables -----------------
NB. SHOULD BE (cutarg arg)=. 'dggev' call ".arg
NB. for now, using dgegv_ ...
(cutarg arg)=. 'dgegv' call ".arg
NB. define result variables ----------------------
alf =. alphar
bet =. beta
lvec=. |:s$vl
rvec=. |:s$vr
if. 1 e. alphai ~: 0 do.
alf =. alf + j. alphai
cx =. I. alphai ~: 0
lvec=. cx cxpair lvec
rvec=. cx cxpair rvec
end.
lvec;alf;bet;rvec NB. ;info,{.work
)
NB. =========================================================
NB. tests
NB. =========================================================
NB.*tdggev v test dggev single test
tdggev=: 3 : 0
match=. matchclean;;
smoutput 'A B' =. y
smoutput clean each 'L ALFA BETA R'=. dggev y
V =. ALFA % BETA
smoutput a=. (clean A mp R) match clean V *"1 B mp R
smoutput b=. (clean (+|:L)mp A) match clean V * (+|:L) mp B
(0 pick a) *. 0 pick b
)
NB. =========================================================
NB.*testdggev v testdggev test sample
NB. cribbed from http://www.nag.co.uk/lapack-ex/node117.html
testdggev=: 3 : 0
A =: >3.9 12.5 _34.5 _0.5; 4.3 21.5 _47.5 7.5; 4.3 21.5 _43.5 3.5; 4.4 26 _46 6
B =: >1 2 _3 1; 1 3 _5 4; 1 3 _4 3; 1 3 _4 4
NB. DGGEV Example Program Results (from website)
NB. Fortran data conventions
NB. Eigenvalue( 1) = 2.0000E+00
NB.
NB. Eigenvector( 1)
NB. 1.0000E+00
NB. 5.7143E-03
NB. 6.2857E-02
NB. 6.2857E-02
NB.
NB. Eigenvalue( 2) = ( 3.0000E+00, 4.0000E+00)
NB.
NB. Eigenvector( 2)
NB. (-4.3979E-01,-5.6021E-01)
NB. (-8.7958E-02,-1.1204E-01)
NB. (-1.4241E-01, 3.1418E-03)
NB. (-1.4241E-01, 3.1418E-03)
NB.
NB. Eigenvalue( 3) = ( 3.0000E+00,-4.0000E+00)
NB.
NB. Eigenvector( 3)
NB. (-4.3979E-01, 5.6021E-01)
NB. (-8.7958E-02, 1.1204E-01)
NB. (-1.4241E-01,-3.1418E-03)
NB. (-1.4241E-01,-3.1418E-03)
NB.
NB. Eigenvalue( 4) = 4.0000E+00
NB.
NB. Eigenvector( 4)
NB. -1.0000E+00
NB. -1.1111E-02
NB. 3.3333E-02
NB. -1.5556E-01
tdggev A;B
)
NB. zggev_ is not installed in jlapack.dll
NB. so will use zgegv_ instead, though "deprecated"
NB. zggev returns generalized eigenvalues, and generalized eigenvectors
NB. for a pair of N-by-N complex nonsymmetric matrices (A,B)
NB. This script is based on
NB. a- fortran source listing for this problem: zggev.f
NB. b- fortran source listing for eigenvalues for square matrix: zgeev.f
NB. c- j script ~\addons\lapack\math\lapack\zgeev.ijs
NB. - as a model of how to call zgeev_
NB. MJL Day 10/07
coclass 'jlapack'
NB. =========================================================
NB.*zggev v generalised eigenvalues and eigenvectors for complex matrices A B
NB.
NB. form: zggev mata;matb
NB.
NB. returns:
NB. left eigenvectors; alpha; beta; right eigenvectors
NB. (all generalised)
NB. The generalized eigenvalue is represented as the pair ( alpha,beta)
NB. Eigenvalue V = alpha%beta
NB. there is a reasonable interpretation for beta = 0, or both zero
NB.
NB. if:
NB. 'L ALFA BETA R' =. zggev A;B
NB. then
NB. V is ALFA % BETA NB. careful....
NB. A mp R is V *"1 B mp R
NB. (+|:L)mp A is V * (+|:L) mp B
zggev=: 3 : 0
'ma mb'=. y
NB. define variables needed for DLL call ---------
s=. $ma
n=. lda=. ldvl=. ldvr=. #ma
a=. zzero + |: ma
b=. zzero + |: mb
ldb =. #mb
alpha =. beta =. n$zzero
jobvl=. jobvr=. 'V'
vl=. vr=. a * zzero
lwork=. 20*n
rwork=. work=. lwork$dzero
info=. izero
NB.
JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHA,BETA,VL,LDVL,VR,LDVR,WORK,LWORK,RWORK,INFO
arg=.'jobvl;jobvr;n;a;lda;b;ldb;alpha;beta;vl;ldvl;vr;ldvr;work;lwork;rwork;info'
NB. call DLL, reassign variables -----------------
NB. SHOULD BE (cutarg arg)=. 'zggev' call ".arg
NB. for now, using zgegv_ ...
(cutarg arg)=. 'zgegv' call ".arg
NB. define result variables ----------------------
lvec=. |:s$vl
rvec=. |:s$vr
lvec;alpha;beta;rvec NB. ;info,{.work
)
NB. =========================================================
NB. tests
NB. =========================================================
NB.*tzggev v test zggev single test
tzggev=: 3 : 0
match=. matchclean;;
smoutput 'A B' =. y
smoutput clean each'L ALFA BETA R'=. zggev y
V =. ALFA % BETA
smoutput a=. (clean A mp R) match clean V *"1 B mp R
smoutput b=. (clean (+|:L)mp A) match clean V * (+|:L) mp B
(0 pick a) *. 0 pick b
)
NB. =========================================================
NB.*testzggev v testzggev test (real) sample
NB. data cribbed from http://www.nag.co.uk/lapack-ex/node122.html
testzggevA=: >( _2 j./\ ".) each LF cut 0 : 0
_21.10 _22.50 53.50 _50.50 _34.50 127.50 7.50 0.50
_0.46 _7.78 _3.50 _37.50 _15.50 58.50 _10.50 _1.50
4.30 _5.50 39.70 _17.10 _68.50 12.50 _7.50 _3.50
5.50 4.40 14.40 43.30 _32.50 _46.00 _19.00 _32.50
)
testzggevB=: >( _2 j./\ ".) each LF cut 0 : 0
1.00 _5.00 1.60 1.20 _3.00 0.00 0.00 _1.00
0.80 _0.60 3.00 _5.00 _4.00 3.00 _2.40 _3.20
1.00 0.00 2.40 1.80 _4.00 _5.00 0.00 _3.00
0.00 1.00 _1.80 2.40 0.00 _4.00 4.00 _5.00
)
testzggev=: 3 : 0
NB. To find all the eigenvalues and right eigenvectors of the matrix pair {A
B} where
NB. data cribbed from http://www.nag.co.uk/lapack-ex/node122.html
NB. ZGGEV Example Program Results (from website)
NB. Fortran data conventions
NB. Eigenvalue( 1) = ( 3.0000E+00,-9.0000E+00)
NB.
NB. Eigenvector( 1)
NB. (-8.2377E-01,-1.7623E-01) (-1.5295E-01, 7.0655E-02)
(-7.0655E-02,-1.5295E-01)
NB. ( 1.5295E-01,-7.0655E-02)
NB.
NB. Eigenvalue( 2) = ( 2.0000E+00,-5.0000E+00)
NB.
NB. Eigenvector( 2)
NB. ( 6.3974E-01, 3.6026E-01) ( 4.1597E-03,-5.4650E-04) ( 4.0212E-02,
2.2645E-02)
NB. (-2.2645E-02, 4.0212E-02)
NB.
NB. Eigenvalue( 3) = ( 3.0000E+00,-1.0000E-00)
NB.
NB. Eigenvector( 3)
NB. ( 9.7754E-01, 2.2465E-02) ( 1.5910E-01,-1.1371E-01) (
1.2090E-01,-1.5371E-01)
NB. ( 1.5371E-01, 1.2090E-01)
NB.
NB. Eigenvalue( 4) = ( 4.0000E+00,-5.0000E+00)
NB.
NB. Eigenvector( 4)
NB. (-9.0623E-01, 9.3766E-02) (-7.4303E-03, 6.8750E-03) (
3.0208E-02,-3.1255E-03)
NB. (-1.4586E-02,-1.4097E-01)
tzggev testzggevA;testzggevB
)
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
