NaN problem solved with asm(finit)

6 files changed, 14 insertions, 137 deletions

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));