acknowledgement to Dylan Simon, updated readme and fixed mapCAux

2 files changed, 33 insertions, 13 deletions

diff --git a/README b/README
index 8d7775b..3a260a4 100644
--- a/README
+++ b/README
@@ -10,15 +10,12 @@ REQUIREMENTS ----------------------------
 2) BLAS/LAPACK (http://www.netlib.org/lapack).
    An optimized implementation is recommended. I have tested:
 
-   - Intel's MKL (http://www.intel.com/cd/software/products).
-     There is a free noncommercial download of MKL for Linux.
-
    - ATLAS (http://math-atlas.sourceforge.net).
      In Ubuntu the required packages are "refblas3-dev", "lapack3-dev",
      and "atlas3-base-dev" (or a version tuned for your machine).
-     Please note that ATLAS currently requires compilation -fviaC in 32bit
-     machines. Otherwise many functions fail, producing strange NaN's.
-     Even with -fvia-C we may get wrong behavior in some cases.
+
+   - Intel's MKL (http://www.intel.com/cd/software/products).
+     There is a free noncommercial download of MKL for Linux.
 
 For ghc-6.8.x you may also need:
 
@@ -59,6 +56,7 @@ Using Intel's MKL:
            $ runhaskell Setup.lhs build
            $ runhaskell Setup.lhs install
 
+More information: http://www.hmatrix.googlepages.com/installation
 
 See below for installation on Windows.
 
@@ -70,7 +68,7 @@ Prelude> Numeric.LinearAlgebra.Tests.runTests 20
 
 Additional tests with big matrices (taking a few minutes):
 
-$ runhaskell examples/experiments bigtests
+$ runhaskell examples/experiments/bigtests
 
 EXAMPLES ------------------------------------------------------
 
@@ -99,10 +97,8 @@ 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 blas/lapack calls. In these machines
-  the library is automatically compiled -fvia-C, which apparently
-  solves the problem.
-  On 64-bit, or using MKL, the default and faster -fasm seems to work well.
+  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
@@ -228,3 +224,6 @@ in the Haskell mailing lists for their help.
 - Don Stewart fixed the implementation of the internal data structures
   to achieve excellent, C-like performance in Haskell functions which
   explicitly work with the elements of vectors and matrices.
+
+- Dylan Alex Simon improved the numeric instances to allow optimized
+  implementations of signum and abs on Vectors.
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 052cafd..77e793e 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -95,6 +95,26 @@ inline double sign(double x) {
     }
 }
 
+inline gsl_complex complex_abs(gsl_complex z) {
+    gsl_complex r;
+    r.dat[0] = gsl_complex_abs(z);
+    r.dat[1] = 0;
+    return r;
+}
+
+inline gsl_complex complex_signum(gsl_complex z) {
+    gsl_complex r;
+    double mag;
+    if (z.dat[0] == 0 && z.dat[1] == 0) {
+        r.dat[0] = 0;
+        r.dat[1] = 0;
+    } else {
+        mag = gsl_complex_abs(z);
+        r.dat[0] = z.dat[0]/mag;
+        r.dat[1] = z.dat[1]/mag;
+    }
+    return r;
+}
 
 #define OP(C,F) case C: { for(k=0;k<xn;k++) rp[k] = F(xp[k]); OK }
 #define OPV(C,E) case C: { for(k=0;k<xn;k++) rp[k] = E; OK }
@@ -124,6 +144,7 @@ int mapR(int code, KRVEC(x), RVEC(r)) {
     }
 }
 
+
 int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
     int k;
     REQUIRES(xn == rn,BAD_SIZE);
@@ -132,7 +153,7 @@ int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
         OP(0,gsl_complex_sin)
         OP(1,gsl_complex_cos)
         OP(2,gsl_complex_tan)
-
+        OP(3,complex_abs)
         OP(4,gsl_complex_arcsin)
         OP(5,gsl_complex_arccos)
         OP(6,gsl_complex_arctan)
@@ -144,7 +165,7 @@ int mapCAux(int code, KGCVEC(x), GCVEC(r)) {
         OP(12,gsl_complex_arctanh)
         OP(13,gsl_complex_exp)
         OP(14,gsl_complex_log)
-
+        OP(15,complex_signum)
         OP(16,gsl_complex_sqrt)
 
         // gsl_complex_arg