From 3d15dffe52629fd621ad548b671c3043caefb0d0 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Mon, 28 Dec 2009 17:59:18 +0000
Subject: additional eig functions

---
 lib/Numeric/LinearAlgebra/Algorithms.hs       | 24 +++++++-
 lib/Numeric/LinearAlgebra/LAPACK.hs           | 85 +++++++++++++++++----------
 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 41 +++++++------
 lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h |  4 +-
 lib/Numeric/LinearAlgebra/Tests.hs            |  4 ++
 lib/Numeric/LinearAlgebra/Tests/Properties.hs |  8 ++-
 6 files changed, 112 insertions(+), 54 deletions(-)

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 1544d7c..3d4a209 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -40,6 +40,7 @@ module Numeric.LinearAlgebra.Algorithms (
     leftSV, rightSV,
 -- ** Eigensystems
     eig, eigSH, eigSH',
+    eigenvalues, eigenvaluesSH, eigenvaluesSH',
 -- ** QR
     qr,
 -- ** Cholesky
@@ -92,6 +93,8 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH''      :: Matrix t -> (Vector Double, Matrix t)
+    eigOnly      :: Matrix t -> Vector (Complex Double)
+    eigOnlySH    :: Matrix t -> Vector Double
     cholSH'      :: Matrix t -> Matrix t
     qr'          :: Matrix t -> (Matrix t, Matrix t)
     hess'        :: Matrix t -> (Matrix t, Matrix t)
@@ -129,16 +132,24 @@ linearSolve = linearSolve'
 linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolveSVD = linearSolveSVD'
 
--- | Eigenvalues and eigenvectors of a general square matrix using lapack's dgeev or zgeev.
+-- | Eigenvalues and eigenvectors of a general square matrix.
 --
 -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
 eig         :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
 eig = eig'
 
+-- | Eigenvalues of a general square matrix.
+eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
+eigenvalues = eigOnly
+
 -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric.
 eigSH'      :: Field t => Matrix t -> (Vector Double, Matrix t)
 eigSH' = eigSH''
 
+-- | Similar to 'eigenvaluesSH' without checking that the input matrix is hermitian or symmetric.
+eigenvaluesSH' :: Field t => Matrix t -> Vector Double
+eigenvaluesSH' = eigOnlySH
+
 -- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
 cholSH      :: Field t => Matrix t -> Matrix t
 cholSH = cholSH'
@@ -187,6 +198,8 @@ instance Field Double where
     ctrans' = trans
     eig' = eigR
     eigSH'' = eigS
+    eigOnly = eigOnlyR
+    eigOnlySH = eigOnlyS
     cholSH' = cholS
     qr' = unpackQR . qrR
     hess' = unpackHess hessR
@@ -203,7 +216,9 @@ instance Field (Complex Double) where
     linearSolveSVD' = linearSolveSVDC Nothing
     ctrans' = conj . trans
     eig' = eigC
+    eigOnly = eigOnlyC
     eigSH'' = eigH
+    eigOnlySH = eigOnlyH
     cholSH' = cholH
     qr' = unpackQR . qrC
     hess' = unpackHess hessC
@@ -211,13 +226,18 @@ instance Field (Complex Double) where
     multiply' = multiplyC
 
 
--- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
+-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
 --
 -- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
 eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
 eigSH m | m `equal` ctrans m = eigSH' m
         | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
 
+-- | Eigenvalues of a complex hermitian or real symmetric matrix.
+eigenvaluesSH :: Field t => Matrix t -> Vector Double
+eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
+                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
+
 -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
 --
 -- If @c = chol m@ then @m == ctrans c \<> c@.
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs
index 2eae9b7..a1ac1cf 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK.hs
+++ b/lib/Numeric/LinearAlgebra/LAPACK.hs
@@ -27,6 +27,7 @@ module Numeric.LinearAlgebra.LAPACK (
     rightSVR, rightSVC, leftSVR, leftSVC,
     -- * Eigensystems
     eigR, eigC, eigS, eigS', eigH, eigH',
+    eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH,
     -- * LU
     luR, luC,
     -- * Cholesky
@@ -198,18 +199,16 @@ leftSVAux f st x = unsafePerformIO $ do
 
 foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM
 foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM
-foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM
-foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: CInt -> TMVM
+foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM
 
-eigAux f st m
-    | r == 1 = (fromList [flatten m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
+eigAux f st m = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        dummy <- createMatrix ColumnMajor 1 1
-        app4 f mat m mat dummy vec l mat v st
+        app3 g mat m vec l mat v st
         return (l,v)
   where r = rows m
+        g ra ca pa = f ra ca pa 0 0 nullPtr
 
 
 -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
@@ -217,26 +216,39 @@ eigAux f st m
 eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
 eigC = eigAux zgeev "eigC" . fmat
 
+eigOnlyAux f st m = unsafePerformIO $ do
+        l <- createVector r
+        app2 g mat m vec l st
+        return l
+  where r = rows m
+        g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr
+
+-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
+-- The eigenvalues are not sorted.
+eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
+eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat
+
 -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
 -- The eigenvectors are the columns of v. The eigenvalues are not sorted.
 eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
 eigR m = (s', v'')
     where (s,v) = eigRaux (fmat m)
-          s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
+          s' = fixeig1 s
           v' = toRows $ trans v
           v'' = fromColumns $ fixeig (toList s') v'
           r = rows m
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
-eigRaux m
-    | r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1)
-    | otherwise = unsafePerformIO $ do
+eigRaux m = unsafePerformIO $ do
         l <- createVector r
         v <- createMatrix ColumnMajor r r
-        dummy <- createMatrix ColumnMajor 1 1
-        app4 dgeev mat m mat dummy vec l mat v "eigR"
+        app3 g mat m vec l mat v "eigR"
         return (l,v)
   where r = rows m
+        g ra ca pa = dgeev ra ca pa 0 0 nullPtr
+
+fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s))
+    where r = dim s
 
 fixeig  []  _ =  []
 fixeig [_] [v] = [comp v]
@@ -246,8 +258,22 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
   where scale = vectorMapValR Scale
 fixeig _ _ = error "fixeig with impossible inputs"
 
+
+-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
+-- The eigenvalues are not sorted.
+eigOnlyR :: Matrix Double -> Vector (Complex Double)
+eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
+
+
 -----------------------------------------------------------------------------
 
+eigSHAux f st m = unsafePerformIO $ do
+        l <- createVector r
+        v <- createMatrix ColumnMajor r r
+        app3 f mat m vec l mat v st
+        return (l,v)
+  where r = rows m
+
 -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
 -- The eigenvectors are the columns of v.
 -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).
@@ -258,16 +284,7 @@ eigS m = (s', fliprl v)
 
 -- | 'eigS' in ascending order
 eigS' :: Matrix Double -> (Vector Double, Matrix Double)
-eigS' m
-    | r == 1 = (fromList [flatten m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        app3 dsyev mat m vec l mat v "eigS"
-        return (l,v)
-  where r = rows m
-
------------------------------------------------------------------------------
+eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
 
 -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
 -- The eigenvectors are the columns of v.
@@ -279,14 +296,20 @@ eigH m = (s', fliprl v)
 
 -- | 'eigH' in ascending order
 eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
-eigH' m
-    | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix ColumnMajor r r
-        app3 zheev mat m vec l mat v "eigH"
-        return (l,v)
-  where r = rows m
+eigH' = eigSHAux (zheev 1) "eigH'" . fmat
+
+
+-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
+-- The eigenvalues are sorted in descending order.
+eigOnlyS :: Matrix Double -> Vector Double
+eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat
+
+-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
+-- The eigenvalues are sorted in descending order.
+eigOnlyH :: Matrix (Complex Double) -> Vector Double
+eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat
+
+vrev = flatten . flipud . reshape 1
 
 -----------------------------------------------------------------------------
 foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index db94cc6..06c2479 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -274,19 +274,20 @@ int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) {
 
 //////////////////// general complex eigensystem ////////////
 
-int eig_l_C(KCMAT(a),CMAT(u), CVEC(s),CMAT(v)) {
+int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac==n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_C");
-    double *B = (double*)malloc(2*n*n*sizeof(double));
+    doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
     CHECK(!B,MEM);
-    memcpy(B,ap,n*n*2*sizeof(double));
-
+    memcpy(B,ap,n*n*sizeof(doublecomplex));
     double *rwork = (double*) malloc(2*n*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
-    char jobvl = ur==1?'N':'V';
-    char jobvr = vr==1?'N':'V';
     integer res;
     // ask for optimal lwork
     doublecomplex ans;
@@ -301,7 +302,7 @@ int eig_l_C(KCMAT(a),CMAT(u), CVEC(s),CMAT(v)) {
              &res);
     lwork = ceil(ans.r);
     //printf("ans = %d\n",lwork);
-    doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     //printf("zgeev\n");
     zgeev_  (&jobvl,&jobvr,
@@ -325,14 +326,16 @@ int eig_l_C(KCMAT(a),CMAT(u), CVEC(s),CMAT(v)) {
 
 int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac == n && (ur==1 || (ur==n && uc==n)) && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(up==NULL || ur==n && uc==n, BAD_SIZE);
+    char jobvl = up==NULL?'N':'V';
+    REQUIRES(vp==NULL || vr==n && vc==n, BAD_SIZE);
+    char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_R");
     double *B = (double*)malloc(n*n*sizeof(double));
     CHECK(!B,MEM);
     memcpy(B,ap,n*n*sizeof(double));
     integer lwork = -1;
-    char jobvl = ur==1?'N':'V';
-    char jobvr = vr==1?'N':'V';
     integer res;
     // ask for optimal lwork
     double ans;
@@ -366,13 +369,14 @@ int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) {
 //////////////////// symmetric real eigensystem ////////////
 
 
-int eig_l_S(KDMAT(a),DVEC(s),DMAT(v)) {
+int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac == n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_S");
     memcpy(vp,ap,n*n*sizeof(double));
     integer lwork = -1;
-    char jobz = vr==1?'N':'V';
     char uplo = 'U';
     integer res;
     // ask for optimal lwork
@@ -399,15 +403,16 @@ int eig_l_S(KDMAT(a),DVEC(s),DMAT(v)) {
 
 //////////////////// hermitian complex eigensystem ////////////
 
-int eig_l_H(KCMAT(a),DVEC(s),CMAT(v)) {
+int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) {
     integer n = ar;
-    REQUIRES(n>=2 && ac==n && sn==n && (vr==1 || (vr==n && vc==n)),BAD_SIZE);
+    REQUIRES(ac==n && sn==n, BAD_SIZE);
+    REQUIRES(vr==n && vc==n, BAD_SIZE);
+    char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_H");
     memcpy(vp,ap,2*n*n*sizeof(double));
     double *rwork = (double*) malloc((3*n-2)*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
-    char jobz = vr==1?'N':'V';
     char uplo = 'U';
     integer res;
     // ask for optimal lwork
@@ -421,7 +426,7 @@ int eig_l_H(KCMAT(a),DVEC(s),CMAT(v)) {
              &res);
     lwork = ceil(ans.r);
     //printf("ans = %d\n",lwork);
-    doublecomplex * work = (doublecomplex*)malloc(lwork*2*sizeof(double));
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     zheev_  (&jobz,&uplo,
              &n,(doublecomplex*)vp,&n,
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index dc7a98f..adf096e 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -68,8 +68,8 @@ int svd_l_C(KCMAT(a),CMAT(u),DVEC(s),CMAT(v));
 int eig_l_C(KCMAT(a),CMAT(u),CVEC(s),CMAT(v));
 int eig_l_R(KDMAT(a),DMAT(u),CVEC(s),DMAT(v));
 
-int eig_l_S(KDMAT(a),DVEC(s),DMAT(v));
-int eig_l_H(KCMAT(a),DVEC(s),CMAT(v));
+int eig_l_S(int,KDMAT(a),DVEC(s),DMAT(v));
+int eig_l_H(int,KCMAT(a),DVEC(s),CMAT(v));
 
 int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x));
 int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x));
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index efc8c40..7284cd9 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -210,6 +210,10 @@ runTests n = do
     test (eigSHProp . cHer)
     test (eigProp   . rSq)
     test (eigProp   . cSq)
+    test (eigSHProp2 . rHer)
+    test (eigSHProp2 . cHer)
+    test (eigProp2   . rSq)
+    test (eigProp2   . cSq)
     putStrLn "------ nullSpace"
     test (nullspaceProp . rM)
     test (nullspaceProp . cM)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index d4dff34..7bb914f 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -31,7 +31,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     nullspaceProp,
     svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,
     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
-    eigProp, eigSHProp,
+    eigProp, eigSHProp, eigProp2, eigSHProp2,
     qrProp,
     hessProp,
     schurProp1, schurProp2,
@@ -192,6 +192,12 @@ eigSHProp m = m <> v |~| v <> real (diag s)
               && m |~| v <> real (diag s) <> ctrans v
     where (s, v) = eigSH m
 
+eigProp2 m = fst (eig m) |~| eigenvalues m
+
+eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
+
+------------------------------------------------------------------
+
 qrProp m = q <> r |~| m && unitary q && upperTriang r
     where (q,r) = qr m
 
-- 
cgit v1.2.3