diag using ST

5 files changed, 10 insertions, 63 deletions

diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index bce9922..1b63dd8 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -101,9 +101,6 @@ check msg f = do
 foreign import ccall "auxi.h gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
 
 ---------------------------------------------------
--- ugly, but my haddock version doesn't understand
--- yet infix type constructors
----------------------------------------------------
 ---------- signatures of the C functions ---------
 --------------------------------------------------
 type PD = Ptr Double                            --
@@ -141,8 +138,3 @@ type TCVM = CInt -> PC -> TM                    --
 type TMCVM = CInt -> CInt -> PD -> TCVM         --
 type TMMCVM = CInt -> CInt -> PD -> TMCVM       --
 --------------------------------------------------
-
-type TauxMul a = CInt -> CInt -> CInt -> Ptr a
-               -> CInt -> CInt -> CInt -> Ptr a
-               -> CInt -> CInt -> Ptr a
-               -> IO CInt
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 63f0a8d..51fb6f8 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -213,19 +213,16 @@ class (Storable a, Floating a) => Element a where
     subMatrixD :: (Int,Int) -- ^ (r0,c0) starting position 
                -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
                -> Matrix a -> Matrix a
-    diagD :: Vector a -> Matrix a
 
 instance Element Double where
     constantD  = constantR
     transdata = transdataR
     subMatrixD = subMatrixR
-    diagD      = diagR
 
 instance Element (Complex Double) where
     constantD  = constantC
     transdata  = transdataC
     subMatrixD = subMatrixC
-    diagD      = diagC
 
 ------------------------------------------------------------------
 
@@ -287,28 +284,6 @@ subMatrix :: Element a
           -> Matrix a -- ^ result
 subMatrix = subMatrixD
 
-
----------------------------------------------------------------------
-
-diagAux fun msg (v@V {dim = n}) = unsafePerformIO $ do
-    m <- createMatrix RowMajor n n
-    app2 fun vec v mat m msg
-    return m
-
--- | diagonal matrix from a real vector
-diagR :: Vector Double -> Matrix Double
-diagR = diagAux c_diagR "diagR"
-foreign import ccall "auxi.h diagR" c_diagR :: TVM
-
--- | diagonal matrix from a real vector
-diagC :: Vector (Complex Double) -> Matrix (Complex Double)
-diagC = diagAux c_diagC "diagC"
-foreign import ccall "auxi.h diagC" c_diagC :: TCVCM
-
--- | creates a square matrix with the given diagonal
-diag :: Element a => Vector a -> Matrix a
-diag = diagD
-
 ------------------------------------------------------------------------
 
 constantAux fun x n = unsafePerformIO $ do
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index bbb8cd4..c449b9a 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -113,32 +113,6 @@ int constantC(gsl_complex* pval, CVEC(r)) {
 }
 
 
-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
-}
-
 int conjugate(KCVEC(x),CVEC(t)) {
     REQUIRES(xn==tn,BAD_SIZE);
     DEBUGMSG("conjugate");
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index 507a30d..b57fdaf 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -18,9 +18,6 @@ int constantC(gsl_complex *val, CVEC(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));
-
 const char * gsl_strerror (const int gsl_errno);
 
 int matrix_fscanf(char*filename, RMAT(a));
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 7d2c564..c56bf3d 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -33,6 +33,7 @@ module Data.Packed.Matrix (
 ) where
 
 import Data.Packed.Internal
+import qualified Data.Packed.ST as ST
 import Data.Packed.Vector
 import Data.List(transpose,intersperse)
 import Data.Array
@@ -73,6 +74,14 @@ fliprl m = fromColumns . reverse . toColumns $ m
 
 ------------------------------------------------------------
 
+-- | Creates a square matrix with a given diagonal.
+diag :: Element a => Vector a -> Matrix a
+diag v = ST.runSTMatrix $ do
+    let d = dim v
+    m <- ST.newMatrix 0 d d
+    mapM_ (\k -> ST.writeMatrix m k k (v@>k)) [0..d-1]
+    return m
+
 {- | creates a rectangular diagonal matrix
 
 @> diagRect (constant 5 3) 3 4
@@ -87,7 +96,7 @@ diagRect s r c
     | r == c    = diag s
     | r < c     = trans $ diagRect s c r
     | otherwise = joinVert  [diag s , zeros (r-c,c)]
-    where zeros (r',c') = reshape c' $ constantD 0 (r'*c')
+    where zeros (r',c') = reshape c' $ constant 0 (r'*c')
 
 -- | extracts the diagonal from a rectangular matrix
 takeDiag :: (Element t) => Matrix t -> Vector t