module Numeric.LinearAlgebra.EigenSolvers.Experimental where
import Control.Exception.Common
import Control.Iterative
import Data.Sparse.Common
import Control.Monad.Catch
import Control.Monad.State.Strict
import Data.VectorSpace
gklBidiag aa q1nn | dim q1nn == n = return (pp, bb, qq)
| otherwise = throwM (MatVecSizeMismatchException "hhBidiag" (dim aa) (dim q1nn))
where
(m,n) = (nrows aa, ncols aa)
aat = transpose aa
bb = fromListSM (n,n) bl
(ql, _, pl, _, pp, bl, qq) = execState (modifyUntil tf bidiagStep) bidiagInit
tf (_, _, _, i, _, _, _) = i == n
bidiagInit = (q2n, beta1, p1n, 1 :: Int, pp, bb', qq)
where
q1 = normalize2' q1nn
p1 = aa #> q1
alpha1 = norm2' p1
p1n = p1 ./ alpha1
q2 = (aat #> p1) ^-^ (alpha1 .* q1)
beta1 = norm2' q2
q2n = q2 ./ beta1
pp = insertCol (zeroSM m n) p1n 0
qq = insertCol (zeroSM n n) q2n 0
bb' = [(0, 0, alpha1)]
bidiagStep (qj , betajm, pjm , j , pp, bb, qq ) =
(qjp, betaj , pj, succ j, pp', bb', qq') where
u = (aa #> qj) ^-^ (betajm .* pjm)
alphaj = norm2' u
pj = u ./ alphaj
v = (aat #> pj) ^-^ (alphaj .* qj)
betaj = norm2' v
qjp = v ./ betaj
pp' = insertCol pp pj j
bb' = [(j1, j, betaj),
(j ,j, alphaj)] ++ bb
qq' = insertCol qq qjp j