Skip to content

Instantly share code, notes, and snippets.

@abrudz
Created December 20, 2017 22:25
Show Gist options
  • Save abrudz/20b0d91b72c7b7c72f9db8e9849e8621 to your computer and use it in GitHub Desktop.
Save abrudz/20b0d91b72c7b7c72f9db8e9849e8621 to your computer and use it in GitHub Desktop.
Eigen function (and its helper functions) as currently shipped with Dyalog APL
∇ 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