Created
December 20, 2017 22:25
-
-
Save abrudz/20b0d91b72c7b7c72f9db8e9849e8621 to your computer and use it in GitHub Desktop.
Eigen function (and its helper functions) as currently shipped with Dyalog APL
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
∇ vv←Eigen mat;vecs;vals;code;n | |
n←≢mat | |
:If Real mat | |
:If (Sym mat)∧2≠n ⍝ symmetric | |
(vals vecs code)←Dsyev mat | |
:Else ⍝ general | |
(vals vecs code)←Dgeev mat | |
:EndIf | |
:Else ⍝ complex | |
:If Sym mat ⍝ hermitian | |
(vals vecs code)←Zheev mat | |
:Else ⍝ general | |
(vals vecs code)←Zgeev mat | |
:EndIf | |
:EndIf | |
⎕SIGNAL code/11 | |
vv←vals Join vecs | |
∇ | |
Join←{⍺⍪⍉n n⍴⍵} ⍝ hang vectors below values | |
Real←1289≠⎕dr ⍝ is mat real? | |
Sym←{w≡+⍉w←⍵+(○+*)≢⎕CT←2*¯32} ⍝ offset coarse tolerance for max fuzz symmetric/hermetian check. | |
J←{⍺+¯11○⍵} ⍝ J's j. verb | |
Call←{(⍎⍺⍺ Assoc ⍵⍵)⍵} ⍝ Call external fn ⍺⍺ with arg ⍺, associate with types ⍵⍵ if necessary | |
Assoc←{3=⎕NC ⍺:⍺ | |
plat←⊃# ⎕WG'APLVersion' | |
bits←¯2↑'32',⎕D∩⍨plat | |
ext←'.so'/⍨'Linux'≡5↑plat | |
call←'lapack',bits,ext,'|',⍺,' ',⍵ | |
0::'*** ERROR ',(⍕⎕EN),': ⎕NA''',call,'''' | |
⎕NA call} | |
Dsyev←{ | |
⍝ JOBZ UPLO N A LDA W WORK LWORK INFO | |
types←'<C1 <C1 <I4 =F8[] <I4 >F8[] >F8[] <I4 >I4' | |
args←'V' 'L'n(,⍵)n n(¯1+3×n)(¯1+3×n)0 | |
2 1 3⊃¨⊂1 1 0 1/('dsyev_'Call types)args ⍝ call external fn. | |
} | |
Dgeev←{ | |
⍝ JOBVL JOBVR N A LDA WR WI VL LDVL VR LDVR WORK LWORK INFO | |
types←'<C1 <C1 <I4 =F8[] <I4 >F8[] >F8[] >I4 <I4 >F8[] <I4 >F8[] <I4 >I4' | |
args←'N' 'V'n(,⍉mat)n n n 0 1(n×n)n(4×n)(4×n)0 | |
(real cmpx vecs code)←0 1 1 0 1 0 1/('dgeev_'Call types)args ⍝ call external fn. | |
vals←real J cmpx ⍝ combine parts | |
∧/0=cmpx:real vecs code ⍝ all-real case | |
vr←⍉n n⍴vecs ⍝ build matrix | |
eii←vr×⍤1⊢cmpx=0 | |
wi←cmpx,0 ⍝ {(⍵≠0)^⍵=-1⌽⍵}WI,0 | |
pair←(wi≠0)∧wi=-1⌽wi ⍝ problem if 3 same in a row! | |
pair←pair/⍳⍴pair | |
eii[;pair]←vr[;pair]J vr[;1+pair] | |
eii[;1+pair]←vr[;pair]J-vr[;1+pair] | |
vecs←⍉eii | |
vals vecs code ⍝ return parts | |
} | |
Zheev←{ | |
⍝ JOBZ UPLO N A LDA W WORK LWORK RWORK INFO | |
types←'<C1 <C1 <I4 =J16[] <I4 >F8[] >F8[] <I4 >F8[] >I4' | |
args←'V' 'L'n(∊⍉mat)n n(¯2+4×n)(¯1+2×n)(¯2+3×n)0 | |
2 1 3⊃¨⊂1 1 0 0 1/('zheev_'Call types)args ⍝ call external fn. | |
} | |
Zgeev←{ | |
⍝ JOBVL JOBVR N A LDA W VL LDVL VR LDVR WORK LWORK RWORK INFO | |
types←'<C1 <C1 <I4 =J16[] <I4 >J16[] <I4 <I4 >J16[] <I4 >F8[] <I4 >F8[] >I4' | |
args←'N' 'V'n(∊⍉mat)n n 0 1(n×n)n(4×n)(2×n)(2×n)0 | |
0 1 1 0 0 1/('zgeev_'Call types)args ⍝ call external fn. | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment