svdR, some quickCheck

6 files changed, 142 insertions, 19 deletions

diff --git a/examples/pru.hs b/examples/pru.hs
index 8b25780..10789d2 100644
--- a/examples/pru.hs
+++ b/examples/pru.hs
@@ -38,7 +38,7 @@ bf = (3>|<4) [7,11,15,8,12,16,9,13,17,10,14,18::Double]
 
 a |=| b = rows a == rows b &&
           cols a == cols b &&
-          toList (dat a) == toList (dat b)
+          toList (cdat a) == toList (cdat b)
 
 mulC a b = multiply RowMajor a b
 mulF a b = multiply ColumnMajor a b
@@ -75,15 +75,14 @@ delta i j | i==j      = 1
 
 e i n = fromList [ delta k i | k <- [1..n]]
 
-ident n = fromRows [ e i n | i <- [1..n]]
+diagl = diag.fromList
 
-diag l = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..c]]
-    where c = length l
+ident n = diag (constant n 1)
 
 tensorFromVector idx v = T {dims = [(dim v,idx)], ten = v}
 tensorFromMatrix idxr idxc m = T {dims = [(rows m,idxr),(cols m,idxc)], ten = cdat m}
 
-td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diag [1..4] :: Tensor Double
+td = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ diagl [1..4] :: Tensor Double
 
 tn = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
 tt = tensorFromMatrix (Contravariant,"i") (Covariant,"j") $ (2><3) [1..6] :: Tensor Double
@@ -114,3 +113,4 @@ names t = sort $ map (snd.snd) (dims t)
 normal t = tridx (names t) t
 
 contractions t1 t2 = [ contraction t1 n1 t2 n2 | n1 <- names t1, n2 <- names t2, compatIdx t1 n1 t2 n2 ]
+
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index db53cd1..bd333d4 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -74,8 +74,7 @@ common f = commonval . map f where
     commonval (a:b:xs) = if a==b then commonval (b:xs) else Nothing
 
 
-toLists m | fortran m = transpose $ partit (rows m) . toList . dat $ m
-          | otherwise = partit (cols m) . toList . dat $ m
+toLists m = partit (cols m) . toList . cdat $ m
 
 dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp
     where
@@ -145,6 +144,8 @@ transdata c1 d c2 | isReal baseOf d = scast $ transdataR c1 (scast d) c2
 --{-# RULES "transdataR" transdata=transdataR #-}
 --{-# RULES "transdataC" transdata=transdataC #-}
 
+-----------------------------------------------------------------------------
+
 -- | creates a Matrix from a list of vectors
 fromRows :: Field t => [Vector t] -> Matrix t
 fromRows vs = case common dim vs of
@@ -160,6 +161,34 @@ toRows m = toRows' 0 where
     toRows' k | k == r*c  = []
               | otherwise = subVector k c v : toRows' (k+c)
 
+-- | Creates a matrix from a list of vectors, as columns
+fromColumns :: Field t => [Vector t] -> Matrix t
+fromColumns m = trans . fromRows $ m
+
+-- | Creates a list of vectors from the columns of a matrix
+toColumns :: Field t => Matrix t -> [Vector t]
+toColumns m = toRows . trans $ m
+
+-- | creates a matrix from a vertical list of matrices
+joinVert :: Field t => [Matrix t] -> Matrix t
+joinVert ms = case common cols ms of
+    Nothing -> error "joinVert on matrices with different number of columns"
+    Just c  -> reshape c $ join (map cdat ms)
+
+-- | creates a matrix from a horizontal list of matrices
+joinHoriz :: Field t => [Matrix t] -> Matrix t
+joinHoriz ms = trans. joinVert . map trans $ ms
+
+------------------------------------------------------------------------------
+
+-- | Reverse rows 
+flipud :: Field t => Matrix t -> Matrix t
+flipud m = fromRows . reverse . toRows $ m
+
+-- | Reverse columns
+fliprl :: Field t => Matrix t -> Matrix t
+fliprl m = fromColumns . reverse . toColumns $ m
+
 -----------------------------------------------------------------
 
 liftMatrix f m = m { dat = f dat, tdat = f tdat } -- check sizes
@@ -168,7 +197,11 @@ liftMatrix f m = m { dat = f dat, tdat = f tdat } -- check sizes
 
 dotL a b = sum (zipWith (*) a b)
 
-multiplyL a b = [[dotL x y | y <- transpose b] | x <- a]
+multiplyL a b | ok = [[dotL x y | y <- transpose b] | x <- a]
+              | otherwise = error "inconsistent dimensions in contraction "
+    where ok = case common length a of
+                   Nothing -> False
+                   Just c  -> c == length b
 
 transL m = matrixFromVector RowMajor (rows m) $ transdataG (cols m) (cdat m) (rows m)
 
@@ -201,9 +234,8 @@ foreign import ccall safe "aux.h multiplyC"
 
 multiply :: (Num a, Field a) => MatrixOrder -> Matrix a -> Matrix a -> Matrix a
 multiply RowMajor a b    = multiplyD RowMajor a b
-multiply ColumnMajor a b = trans $ multiplyT ColumnMajor a b
-
-multiplyT order a b = multiplyD order (trans b) (trans a)
+multiply ColumnMajor a b = m {rows = cols m, cols = rows m, order = ColumnMajor}
+    where m = multiplyD RowMajor (trans b) (trans a)
 
 multiplyD order a b
     | isReal (baseOf.dat) a = scast $ multiplyAux order cmultiplyR (scast a) (scast b)
@@ -253,3 +285,33 @@ subMatrix st sz m
 
 subMatrixG (r0,c0) (rt,ct) x = reshape ct $ fromList $ concat $ map (subList c0 ct) (subList r0 rt (toLists x))
     where subList s n = take n . drop s
+
+---------------------------------------------------------------------
+
+diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
+    m <- createMatrix RowMajor n n
+    fun // vec v // mat dat m // check msg [dat m]
+    return m
+
+-- | diagonal matrix from a real vector
+diagR :: Vector Double -> Matrix Double
+diagR = diagAux c_diagR "diagR"
+foreign import ccall "aux.h diagR" c_diagR :: Double :> Double ::> IO Int
+
+-- | diagonal matrix from a real vector
+diagC :: Vector (Complex Double) -> Matrix (Complex Double)
+diagC = diagAux c_diagC "diagC"
+foreign import ccall "aux.h diagC" c_diagC :: (Complex Double) :> (Complex Double) ::> IO Int
+
+-- | diagonal matrix from a vector
+diag :: (Num a, Field a) => Vector a -> Matrix a
+diag v
+    | isReal (baseOf) v = scast $ diagR (scast v)
+    | isComp (baseOf) v = scast $ diagC (scast v)
+    | otherwise             = diagG v
+
+diagG v = reshape c $ fromList $ [ l!!(i-1) * delta k i | k <- [1..c], i <- [1..c]]
+    where c = dim v
+          l = toList v
+          delta i j | i==j      = 1
+                    | otherwise = 0
diff --git a/lib/Data/Packed/Internal/aux.c b/lib/Data/Packed/Internal/aux.c
index 01a2bb3..fe611e2 100644
--- a/lib/Data/Packed/Internal/aux.c
+++ b/lib/Data/Packed/Internal/aux.c
@@ -18,19 +18,17 @@
 #define MACRO(B) do {B} while (0)
 #define ERROR(CODE) MACRO(return CODE;)
 #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
+#define OK return 0;
 
 #define MIN(A,B) ((A)<(B)?(A):(B))
 #define MAX(A,B) ((A)>(B)?(A):(B))
 
 #ifdef DBG
-#define DEBUGMSG(M) printf("GSL Wrapper "M": "); size_t t0 = time(NULL);
-#define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;);
+#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
 #else
 #define DEBUGMSG(M)
-#define OK return 0;
 #endif
 
-
 #define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
 
 #ifdef DBG
@@ -45,7 +43,6 @@
 #define DEBUGVEC(MSG,X)
 #endif
 
-
 #define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
 #define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
 #define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
@@ -66,8 +63,6 @@
 #define MEM      1002
 #define BAD_FILE 1003
 
-
-
 int transR(KRMAT(x),RMAT(t)) {
     REQUIRES(xr==tc && xc==tr,BAD_SIZE);
     DEBUGMSG("transR");
@@ -122,6 +117,7 @@ int constantC(gsl_complex* pval, CVEC(r)) {
     OK
 }
 
+
 int multiplyR(int ta, KRMAT(a), int tb, KRMAT(b),RMAT(r)) {
     //printf("%d %d %d %d %d %d\n",ar,ac,br,bc,rr,rc);
     //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
@@ -155,3 +151,30 @@ int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r)) {
     CHECK(res,res);
     OK
 }
+
+
+int diagR(KRVEC(d),RMAT(r)) {
+    REQUIRES(dn==rr && rr==rc,BAD_SIZE);
+    DEBUGMSG("diagR");
+    int i,j;
+    for (i=0;i<rr;i++) {
+        for(j=0;j<rc;j++) {
+            rp[i*rc+j] = i==j?dp[i]:0.;
+        }
+    }
+    OK
+}
+
+int diagC(KCVEC(d),CMAT(r)) {
+    REQUIRES(dn==rr && rr==rc,BAD_SIZE);
+    DEBUGMSG("diagC");
+    int i,j;
+    gsl_complex zero;
+    GSL_SET_COMPLEX(&zero,0.,0.);
+    for (i=0;i<rr;i++) {
+        for(j=0;j<rc;j++) {
+            rp[i*rc+j] = i==j?dp[i]:zero;
+        }
+    }
+    OK
+}
diff --git a/lib/Data/Packed/Internal/aux.h b/lib/Data/Packed/Internal/aux.h
index 59d90db..d055d35 100644
--- a/lib/Data/Packed/Internal/aux.h
+++ b/lib/Data/Packed/Internal/aux.h
@@ -21,3 +21,6 @@ int multiplyR(int ta, KRMAT(a), int tb,  KRMAT(b),RMAT(r));
 int multiplyC(int ta, KCMAT(a), int tb, KCMAT(b),CMAT(r));
 
 int submatrixR(int r1, int r2, int c1, int c2, KRMAT(x),RMAT(r));
+
+int diagR(KRVEC(d),RMAT(r));
+int diagC(KCVEC(d),CMAT(r));
diff --git a/lib/LAPACK.hs b/lib/LAPACK.hs
index cab437c..67e49af 100644
--- a/lib/LAPACK.hs
+++ b/lib/LAPACK.hs
@@ -14,6 +14,7 @@
 
 module LAPACK (
     --module LAPACK.Internal
+    svdR, svdR',
     eigC,
     linearSolveLSR
 ) where
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs
index 4c755bc..2569215 100644
--- a/lib/LAPACK/Internal.hs
+++ b/lib/LAPACK/Internal.hs
@@ -27,6 +27,24 @@ import Foreign.C.String
 foreign import ccall "lapack-aux.h svd_l_R"
     dgesvd :: Double ::> Double ::> (Double :> Double ::> IO Int)
 
+-- | Wrapper for LAPACK's /dgesvd/, which computes the full svd decomposition of a real matrix.
+--
+-- @(u,s,v)=svdR m@ so that @m=u \<\> s \<\> 'trans' v@.
+svdR :: Matrix Double -> (Matrix Double, Matrix Double , Matrix Double)
+svdR x@M {rows = r, cols = c} = (u, s, v)
+    where (u,s',v) = svdR' x
+          s | r == c    = diag s'
+            | r < c     = joinHoriz [diag s' , zeros (r,c-r)]
+            | otherwise = joinVert  [diag s' , zeros (r-c,c)]
+          zeros (r,c) = reshape c $ constant (r*c) 0
+
+svdR' x@M {rows = r, cols = c} = unsafePerformIO $ do
+    u <- createMatrix ColumnMajor r r
+    s <- createVector (min r c)
+    v <- createMatrix ColumnMajor c c
+    dgesvd // mat fdat x // mat dat u // vec s // mat dat v // check "svdR" [fdat x]
+    return (u,s,trans v)
+
 -----------------------------------------------------------------------------
 -- dgesdd
 foreign import ccall "lapack-aux.h svd_l_Rdd"
@@ -62,11 +80,27 @@ foreign import ccall "lapack-aux.h eig_l_R"
     dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int)
 
 -----------------------------------------------------------------------------
-
 -- dsyev
 foreign import ccall "lapack-aux.h eig_l_S"
     dsyev :: Double ::> (Double :> Double ::> IO Int)
 
+-- | Wrapper for LAPACK's /dsyev/, which computes the eigenvalues and right eigenvectors of a symmetric real matrix:
+--
+-- if @(l,v)=eigSl m@ then @m \<\> v = v \<\> diag l@.
+--
+-- The eigenvectors are the columns of v.
+-- The eigenvalues are sorted in descending order (use eigS' for ascending order).
+eigS :: Matrix Double -> (Vector Double, Matrix Double)
+eigS m = (s', fliprl v)
+    where (s,v) = eigS' m
+          s' = fromList . reverse . toList $  s
+
+eigS' (m@M {rows = r}) = unsafePerformIO $ do
+    l <- createVector r
+    v <- createMatrix ColumnMajor r r
+    dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m]
+    return (l,v)
+
 -----------------------------------------------------------------------------
 -- zheev
 foreign import ccall "lapack-aux.h eig_l_H"