3 files changed, 46 insertions, 9 deletions

diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 01d2ccf..68547bd 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -23,6 +23,7 @@ import Data.Packed.Internal.Vector
 import Foreign hiding (xor)
 import Complex
 import Foreign.C.Types
+import Foreign.C.String
 
 -----------------------------------------------------------------
 
@@ -350,6 +351,19 @@ fromComplex :: Vector (Complex Double) -> (Vector Double, Vector Double)
 fromComplex z = (r,i) where
     [r,i] = toColumns $ reshape 2 $ asReal z
 
--- | loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
-fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
-fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+--------------------------------------------------------------------------
+
+-- | Saves a matrix as 2D ASCII table.
+saveMatrix :: FilePath
+           -> String     -- ^ format (%f, %g, %e)
+           -> Matrix Double
+           -> IO ()
+saveMatrix filename fmt m = do
+    charname <- newCString filename
+    charfmt <- newCString fmt
+    let o = if orderOf m == RowMajor then 1 else 0
+    app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
+    free charname
+    free charfmt
+
+foreign import ccall "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 32eb991..ad0776e 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -30,7 +30,7 @@ module Data.Packed.Matrix (
     ident, diag, diagRect, takeDiag,
     liftMatrix, liftMatrix2,
     format,
-    loadMatrix, fromFile, fileDimensions,
+    loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix, fromArray2D
 ) where
 
@@ -247,11 +247,15 @@ fileDimensions fname = do
         then return (tot `div` c, c)
         else return (0,0)
 
-{- | loads a matrix from a formatted ASCII file.
--}
+-- | Loads a matrix from an ASCII file formatted as a 2D table.
 loadMatrix :: FilePath -> IO (Matrix Double)
 loadMatrix file = fromFile file =<< fileDimensions file
 
+-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
+fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
+fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
+
+
 -- | rearranges the rows of a matrix according to the order given in a list of integers. 
 extractRows :: Element t => [Int] -> Matrix t -> Matrix t
 extractRows l m = fromRows $ extract (toRows $ m) l
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index b3b524d..d129aeb 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -351,7 +351,7 @@ int polySolve(KRVEC(a), CVEC(z)) {
 }
 
 int vector_fscanf(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_matrix_fscanf");
+    DEBUGMSG("gsl_vector_fscanf");
     DVVIEW(a);
     FILE * f = fopen(filename,"r");
     CHECK(!f,BAD_FILE);
@@ -373,7 +373,7 @@ int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
 }
 
 int vector_fread(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_matrix_fscanf");
+    DEBUGMSG("gsl_vector_fread");
     DVVIEW(a);
     FILE * f = fopen(filename,"r");
     CHECK(!f,BAD_FILE);
@@ -384,7 +384,7 @@ int vector_fread(char*filename, RVEC(a)) {
 }
 
 int vector_fwrite(char*filename, RVEC(a)) {
-    DEBUGMSG("gsl_vector_fprintf");
+    DEBUGMSG("gsl_vector_fwrite");
     DVVIEW(a);
     FILE * f = fopen(filename,"w");
     CHECK(!f,BAD_FILE);
@@ -394,6 +394,25 @@ int vector_fwrite(char*filename, RVEC(a)) {
     OK
 }
 
+int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
+    DEBUGMSG("matrix_fprintf");
+    FILE * f = fopen(filename,"w");
+    CHECK(!f,BAD_FILE);
+    int i,j,sr,sc;
+    if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
+    #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
+    for (i=0; i<mr; i++) {
+        for (j=0; j<mc-1; j++) {
+            fprintf(f,fmt,AT(m,i,j));
+            fprintf(f," ");
+        }
+        fprintf(f,fmt,AT(m,i,j));
+        fprintf(f,"\n");
+    }
+    fclose (f);
+    OK
+}
+
 //---------------------------------------------------------------
 
 typedef double Trawfun(int, double*);