Skip to content

Instantly share code, notes, and snippets.

@abrudz
Created December 20, 2017 22:19
Show Gist options
  • Save abrudz/ab6f4e2f4ac7da8a97b8dcbd0174ea9a to your computer and use it in GitHub Desktop.
Save abrudz/ab6f4e2f4ac7da8a97b8dcbd0174ea9a to your computer and use it in GitHub Desktop.
Eigen function in as it used to be in old versions of Dyalog APL
Eigen←{⎕IO ⎕ML←1 ⍝ Eigen values⍪vectors
⍝ graeme@dyadic.com
dsyev←{ ⍝ real symmetric:
args←'dsyev_'assoc↓⍉↑⌽{ ⍝ associate external fn.
⍵,⊂' <T1 ' 'V'}{ ⍝ JOBZ
⍵,⊂' <T1 ' 'L'}{ ⍝ UPLO
⍵,⊂' <I4 'n}{ ⍝ N
⍵,⊂' =F8[] '(,mat)}{ ⍝ A
⍵,⊂' <I4 'n}{ ⍝ LDA
⍵,⊂' >F8[] 'n}{ ⍝ W
⍵,⊂' >F8[] '(¯1+3×n)}{ ⍝ WORK
⍵,⊂' <I4 '(¯1+3×n)}{ ⍝ LWORK
⍵,⊂' >I4 ' 0}⍬ ⍝ INFO
vecs vals _ code←dsyev_ args ⍝ call external fn.
code=0:vals join vecs ⍝ good rslt: values⍪vectors.
⎕SIGNAL 11 ⍝ else: DOMAIN ERROR.
}
dgeev←{ ⍝ real general:
args←'dgeev_'assoc↓⍉↑⌽{ ⍝ associate external fn.
⍵,⊂' <C1 ' 'N'}{ ⍝ JOBVL
⍵,⊂' <C1 ' 'V'}{ ⍝ JOBVR
⍵,⊂' <I4 'n}{ ⍝ N
⍵,⊂' =F8[] '(,⍉mat)}{ ⍝ A
⍵,⊂' <I4 'n}{ ⍝ LDA
⍵,⊂' >F8[] 'n}{ ⍝ WR
⍵,⊂' >F8[] 'n}{ ⍝ WI
⍵,⊂' >I4 ' 0}{ ⍝ VL
⍵,⊂' <I4 ' 1}{ ⍝ LDVL
⍵,⊂' >F8[] '(n×n)}{ ⍝ VR
⍵,⊂' <I4 'n}{ ⍝ LDVR
⍵,⊂' >F8[] '(4×n)}{ ⍝ WORK
⍵,⊂' <I4 '(4×n)}{ ⍝ LWORK
⍵,⊂' >I4 ' 0}⍬ ⍝ INFO
_ vals cmpx _ vecs _ code←dgeev_ args ⍝ call external fn.
code≠0:⎕SIGNAL 11 ⍝ DOMAIN ERROR.
0∧.=cmpx:vals join vecs ⍝ real
VR←⍉n n⍴vecs
Ei←vals,¨cmpx
Eii←VR×[2]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],¨VR[;1+pair]
Eii[;1+pair]←VR[;pair],¨-VR[;1+pair]
Ei⍪Eii
}
zheev←{ ⍝ complex hermitian:
args←'zheev_'assoc↓⍉↑⌽{ ⍝ associate external fn.
⍵,⊂' <C1 ' 'V'}{ ⍝ JOBZ
⍵,⊂' <C1 ' 'L'}{ ⍝ UPLO
⍵,⊂' <I4 'n}{ ⍝ N
⍵,⊂' =F8[] '(∊⍉mat)}{ ⍝ A
⍵,⊂' <I4 'n}{ ⍝ LDA
⍵,⊂' >F8[] 'n}{ ⍝ W
⍵,⊂' >F8[] '(¯2+4×n)}{ ⍝ WORK
⍵,⊂' <I4 '(¯1+2×n)}{ ⍝ LWORK
⍵,⊂' >F8[] '(¯2+3×n)}{ ⍝ RWORK
⍵,⊂' >I4 ' 0}⍬ ⍝ INFO
vecs vals _ _ code←zheev_ args ⍝ call external fn.
code=0:vals join cmpx vecs ⍝ good rslt: values⍪vectors.
⎕SIGNAL 11 ⍝ else: DOMAIN ERROR.
}
zgeev←{ ⍝ complex general:
args←'zgeev_'assoc↓⍉↑⌽{ ⍝ associate external fn.
⍵,⊂' <C1 ' 'N'}{ ⍝ JOBVL
⍵,⊂' <C1 ' 'V'}{ ⍝ JOBVR
⍵,⊂' <I4 'n}{ ⍝ N
⍵,⊂' =F8[] '(∊⍉mat)}{ ⍝ A
⍵,⊂' <I4 'n}{ ⍝ LDA
⍵,⊂' >F8[] '(2×n)}{ ⍝ W
⍵,⊂' <I4 ' 0}{ ⍝ VL
⍵,⊂' <I4 ' 1}{ ⍝ LDVL
⍵,⊂' >F8[] '(2×n×n)}{ ⍝ VR
⍵,⊂' <I4 'n}{ ⍝ LDVR
⍵,⊂' >F8[] '(4×n)}{ ⍝ WORK
⍵,⊂' <I4 '(2×n)}{ ⍝ LWORK
⍵,⊂' >F8[] '(2×n)}{ ⍝ RWORK
⍵,⊂' >I4 ' 0}⍬ ⍝ INFO
_ vals vecs _ _ code←zgeev_ args ⍝ call external fn.
code=0:(cmpx vals)join cmpx vecs ⍝ good rslt: values⍪vectors.
⎕SIGNAL 11 ⍝ else: DOMAIN ERROR.
}
assoc←{ ⍝ associate external function.
types args←⍵ ⍝ argument types and values.
3=⎕NC ⍺:args ⍝ function already exists.
{args}⎕NA'lapack|',⍺,∊types ⍝ associate external function.
}
fuzz←{ ⍝ fuzzy comparison.
⎕CT←2*¯32 ⍝ coarse tolerance for max fuzz.
noise←(○1)+*1 ⍝ offset for comparison.
(⍺+noise)⍺⍺ ⍵+noise ⍝ tolerant comparison.
}
join←{⍺⍪⍉n n⍴⍵} ⍝ join real vals and vecs.
cmpx←{↓(0.5×⊃⍴⍵)2⍴⍵} ⍝ complex from real.
n←⊃⍴mat←⍵ ⍝ matrix size and content.
1=≡⍵:{ ⍝ real matrix.
(⍵≡fuzz⍉⍵)∧2≠⊃⍴⍵:dsyev ⍵ ⍝ symmetric (not 2×2).
dgeev ⍵ ⍝ non-symmetric.
}⍵
⍝ complex matrix.
⍵≡fuzz⍉⍵×⊂1 ¯1:zheev ⍵ ⍝ hermitian:
zgeev ⍵ ⍝ non-hermitian.
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment