8 files changed, 19 insertions, 165 deletions

diff --git a/README b/README
index 3a260a4..efb32da 100644
--- a/README
+++ b/README
@@ -96,10 +96,6 @@ A number of illustrative programs are included in the examples folder.
 
 KNOWN PROBLEMS / BUGS -------------------------------
 
-- Compilation with -O -fasm on 32-bit machines produces strange
-  NaN's results on certain foreign calls. More info at
-  http://www.hmatrix.googlepages.com/bugs
-
 - On 64-bit machines the example "minimize.hs", when run from ghci,
   produces a segmentation fault. It happens in the call to
   gsl_multimin_fdfminimizer_alloc, inside the C wrapper.
@@ -108,17 +104,8 @@ KNOWN PROBLEMS / BUGS -------------------------------
   program seems to work perfectly well.
 
 - On Ubuntu 6.06 LTS (Dapper) atlas3-sse2-dev (3.6.0-20)
-  produces segmentation faults when working with big matrices 
-  on compiled programs. To expose the problem:
-
-  $ cd examples
-  $ ghc --make -O -fvia-C tests.hs
-  $ ./tests --big
-
-  If this crashes, just uninstall atlas3-sse2 and use atlas3-base-dev instead.
-  Fortunately, atlas3-sse2-dev seems to work well on Ubuntu 7.10 Gutsy.
-  A similar problem was reported at:
-      http://article.gmane.org/gmane.linux.debian.devel.bugs.general/323065
+  produced segmentation faults when working with big matrices
+  on compiled programs.
 
 - On distributions with old GSL versions you should comment out a couple of functions
   in the export lists of Ellint.hs and Debye.hs
@@ -227,3 +214,6 @@ in the Haskell mailing lists for their help.
 
 - Dylan Alex Simon improved the numeric instances to allow optimized
   implementations of signum and abs on Vectors.
+
+- Pedro E. López de Teruel discovered the need of asm("finit") to
+  avoid the wrong NaNs produced by foreign functions.
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 1ffcd4d..6bffaf9 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -29,13 +29,6 @@ flag unsafe
     description: Compile the library with bound checking disabled.
     default: False
 
-flag no-workaround
-    description:    Expose the NaN problem
-    default: False
-
-flag expose-bug
-    description:    Expose the NaN problem with debug information
-    default: False
 
 library
     if flag(splitBase)
@@ -113,12 +106,6 @@ library
     if flag(unsafe)
         cpp-options: -DUNSAFE
 
-    if arch(i386) && !flag(no-workaround) && !flag(expose-bug) && !flag(mkl)
-        cpp-options: -DWORKAROUND
-
-    if flag(expose-bug)
-        cpp-options: -DEXPOSEBUG
-
 
     if flag(mkl)
       if arch(x86_64)
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
index 3609bc2..bce9922 100644
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -80,9 +80,14 @@ errorCode 2006 = "the input matrix is not positive definite"
 errorCode 2007 = "not yet supported in this OS"
 errorCode n    = "code "++show n
 
+
+-- | clear the fpu
+foreign import ccall "auxi.h asm_finit" finit :: IO ()
+
 -- | check the error code
 check :: String -> IO CInt -> IO ()
 check msg f = do
+    finit
     err <- f
     when (err/=0) $ if err > 1024
                       then (error (msg++": "++errorCode err)) -- our errors
diff --git a/lib/Data/Packed/Internal/auxi.c b/lib/Data/Packed/Internal/auxi.c
index 04dc7ad..bbb8cd4 100644
--- a/lib/Data/Packed/Internal/auxi.c
+++ b/lib/Data/Packed/Internal/auxi.c
@@ -149,3 +149,9 @@ int conjugate(KCVEC(x),CVEC(t)) {
     }
     OK
 }
+
+//---------------------------------------
+void asm_finit() {
+    asm("finit");
+}
+//---------------------------------------
diff --git a/lib/Data/Packed/Internal/auxi.h b/lib/Data/Packed/Internal/auxi.h
index 377a4a1..507a30d 100644
--- a/lib/Data/Packed/Internal/auxi.h
+++ b/lib/Data/Packed/Internal/auxi.h
@@ -26,3 +26,5 @@ const char * gsl_strerror (const int gsl_errno);
 int matrix_fscanf(char*filename, RMAT(a));
 
 int conjugate(KCVEC(x),CVEC(t));
+
+void asm_finit();
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 45298b5..e2fec9d 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -139,16 +139,8 @@ instance Field (Complex Double) where
     qr = unpackQR . qrC
     hess = unpackHess hessC
     schur = schurC
-
-#if defined(WORKAROUND)
-    multiply = mulCW
-#else
-#if defined(EXPOSEBUG)
-    multiply = multiplyC
-#else
     multiply = multiplyC3
-#endif
-#endif
+
 
 -- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix using lapack's dsyev or zheev.
 --
@@ -639,47 +631,3 @@ multiplyC3 x y = unsafePerformIO $ multiply3 zgemm "cblas_zgemm" (fmat x) (fmat
                 return s
                 -- if toLists s== toLists s then return s else error $ "HORROR " ++ (show (toLists s))
             | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
-
------------------------------------------------------------------------------------
--- BLAS via auxiliary C
------------------------------------------------------------------------------------
-
--- foreign import ccall "multiply.h multiplyR2" dgemmc :: TMMM
--- foreign import ccall "multiply.h multiplyC2" zgemmc :: TCMCMCM
--- 
--- multiply2 f st a b
---     | cols a == rows b = unsafePerformIO $ do
---         s <- createMatrix ColumnMajor (rows a) (cols b)
---         app3 f mat a mat b mat s st 
---         if toLists s== toLists s then return s else error $ "AYYY " ++ (show (toLists s))
---     | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
--- 
--- multiplyR2 :: Matrix Double -> Matrix Double -> Matrix Double
--- multiplyR2 a b = multiply2 dgemmc "dgemmc" (fmat a) (fmat b)
--- 
--- multiplyC2 :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
--- multiplyC2 a b = multiply2 zgemmc "zgemmc" (fmat a) (fmat b)
-
------------------------------------------------------------------------------------
--- direct C multiplication, to expose the NaN bug
------------------------------------------------------------------------------------
-
--- foreign import ccall "multiply.h multiplyR" cmultiplyR :: TMMM
-foreign import ccall "multiply.h multiplyC" cmultiplyC :: TCMCMCM
-
-cmultiply f st a b
---    | cols a == rows b = 
-      = unsafePerformIO $ do
-        s <- createMatrix RowMajor (rows a) (cols b)
-        app3 f mat a mat b mat s st
-        if toLists s== toLists s
-            then return s
-            else error $ "NaN FOUND!! " ++ (show (toLists s))
-        -- return s
---    | otherwise = error $ st ++ " (matrix product) of nonconformant matrices"
-
--- multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
--- multiplyR a b = cmultiply cmultiplyR "cmultiplyR" (cmat a) (cmat b)
-
-multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-multiplyC a b = cmultiply cmultiplyC "cmultiplyR" (cmat a) (cmat b)
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
index 61cb002..310f6ee 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c
@@ -814,81 +814,3 @@ int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) {
     free(auxipiv);
     OK
 }
-
-////////////////////////////////////////////////////////////
-
-// int multiplyR(KDMAT(a),KDMAT(b),DMAT(r)) {
-//     REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-//     int i,j,k;
-//     for (i=0;i<ar;i++) {
-//         for(j=0;j<bc;j++) {
-//             double temp = 0;
-//             for(k=0;k<ac;k++) {
-//                 temp += ap[i*ac+k]*bp[k*bc+j];
-//             }
-//             rp[i*rc+j] = temp;
-//         }
-//     }
-//     OK
-// }
-
-
-int multiplyC(KCMAT(a),KCMAT(b),CMAT(r)) {
-    REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-    int i,j,k;
-    for (i=0;i<ar;i++) {
-        for(j=0;j<bc;j++) {
-            doublecomplex temp = {0,0};
-            for(k=0;k<ac;k++) {
-                doublecomplex A = ((doublecomplex*)ap)[i*ac+k];
-                doublecomplex B = ((doublecomplex*)bp)[k*bc+j];
-                double w = A.r * B.r - A.i * B.i;
-                double w1 = A.r * B.r;
-                double w2 = A.i * B.i;
-                if(w != w) {
-                    printf("at : %d %d %d\n",i,j,k);
-                    printf("%f %f %f\n",A.i,B.i, A.i * B.i);
-                    printf("%f %f %f\n",A.i,B.i, w2);
-                    exit(1);
-                }
-                temp.r += (w + w1-w2)/2;
-                //temp.r += w;
-                temp.i += A.r * B.i + A.i * B.r;
-            }
-            ((doublecomplex*)rp)[i*rc+j] = temp;
-        }
-    }
-    OK
-}
-
-// void dgemm_(char *, char *, integer *, integer *, integer *,
-//            double *, const double *, integer *, const double *,
-//            integer *, double *, double *, integer *);
-// 
-// void zgemm_(char *, char *, integer *, integer *, integer *,
-//            doublecomplex *, const doublecomplex *, integer *, const doublecomplex *,
-//            integer *, doublecomplex *, doublecomplex *, integer *);
-// 
-// 
-// int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r)) {
-//     REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-//     double alpha = 1;
-//     double beta = 0;
-//     integer m = ar;
-//     integer n = bc;
-//     integer k = ac;
-//     int i,j;
-//     dgemm_("N","N",&m,&n,&k,&alpha,ap,&m,bp,&k,&beta,rp,&m);
-//     OK
-// }
-// 
-// int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r)) {
-//     REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE);
-//     integer m = ar;
-//     integer n = bc;
-//     integer k = ac;
-//     doublecomplex alpha = {1,0};
-//     doublecomplex beta = {0,0};
-//     zgemm_("N","N",&m,&n,&k,&alpha,(doublecomplex*)ap,&m,(doublecomplex*)bp,&k,&beta,(doublecomplex*)rp,&m);
-//     OK
-// }
diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
index e8cba30..79e52be 100644
--- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
+++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h
@@ -84,9 +84,3 @@ int schur_l_C(KCMAT(a), CMAT(u), CMAT(s));
 
 int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r));
 int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r));
-
-// int multiplyR(KDMAT(a),KDMAT(b),DMAT(r));
-// int multiplyC(KCMAT(a),KCMAT(b),CMAT(r));
-
-// int multiplyR2(KDMAT(a),KDMAT(b),DMAT(r));
-int multiplyC2(KCMAT(a),KCMAT(b),CMAT(r));