{-# language TypeFamilies, FlexibleContexts #-}
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


-- | `eigRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm` and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair.

eigRayleigh nitermax debq prntf m = untilConvergedGM "eigRayleigh" config (const True) (rayStep m)
  where
    ii = eye (nrows m)
    config = IterConf nitermax debq fst prntf
    rayStep aa (b, mu) = do
      nom <- (m ^-^ (mu `matScale` ii)) <\> b
      let b' = normalize2' nom
          mu' = (b' <.> (aa #> b')) / (b' <.> b')
      return (b', mu')

      

-- | Golub-Kahan-Lanczos bidiagonalization (see "Restarted Lanczos Bidiagonalization for the SVD", SLEPc STR-8, http://slepc.upv.es/documentation/manual.htm )
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' = [(j-1, j, betaj),
           (j ,j, alphaj)] ++ bb
    qq' = insertCol qq qjp j