From 192ac5f4b98517862c37ecf161505396ad223cd8 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Thu, 2 Oct 2008 15:53:10 +0000
Subject: alternative multiply versions

---
 lib/Numeric/GSL/Matrix.hs | 311 ----------------------------------------------
 lib/Numeric/GSL/gsl-aux.c | 286 ------------------------------------------
 lib/Numeric/GSL/gsl-aux.h |  19 ---
 3 files changed, 616 deletions(-)
 delete mode 100644 lib/Numeric/GSL/Matrix.hs

(limited to 'lib/Numeric/GSL')

diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
deleted file mode 100644
index f62bb82..0000000
--- a/lib/Numeric/GSL/Matrix.hs
+++ /dev/null
@@ -1,311 +0,0 @@
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.GSL.Matrix
--- Copyright   :  (c) Alberto Ruiz 2007
--- License     :  GPL-style
---
--- Maintainer  :  Alberto Ruiz <aruiz@um.es>
--- Stability   :  provisional
--- Portability :  portable (uses FFI)
---
--- A few linear algebra computations based on GSL.
---
------------------------------------------------------------------------------
--- #hide
-
-module Numeric.GSL.Matrix(
-    eigSg, eigHg,
-    svdg,
-    qr, qrPacked, unpackQR,
-    cholR, cholC,
-    luSolveR, luSolveC,
-    luR, luC
-) where
-
-import Data.Packed.Internal
-import Data.Packed.Matrix(ident)
-import Numeric.GSL.Vector
-import Foreign
-import Complex
-
-{- | eigendecomposition of a real symmetric matrix using /gsl_eigen_symmv/.
-
-> > let (l,v) = eigS $ 'fromLists' [[1,2],[2,1]]
-> > l
-> 3.000 -1.000
->
-> > v
-> 0.707 -0.707
-> 0.707  0.707
->
-> > v <> diag l <> trans v
-> 1.000 2.000
-> 2.000 1.000
-
--}
-eigSg :: Matrix Double -> (Vector Double, Matrix Double)
-eigSg = eigSg' . cmat
-
-eigSg' m
-    | r == 1 = (fromList [cdat m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix RowMajor r r
-        app3 c_eigS mat m vec l mat v "eigSg"
-        return (l,v)
-  where r = rows m
-foreign import ccall "gsl-aux.h eigensystemR" c_eigS :: TMVM
-
-------------------------------------------------------------------
-
-
-
-{- | eigendecomposition of a complex hermitian matrix using /gsl_eigen_hermv/
-
-> > let (l,v) = eigH $ 'fromLists' [[1,2+i],[2-i,3]]
->
-> > l
-> 4.449 -0.449
->
-> > v
->         -0.544          0.839
-> (-0.751,0.375) (-0.487,0.243)
->
-> > v <> diag l <> (conjTrans) v
->          1.000 (2.000,1.000)
-> (2.000,-1.000)         3.000
-
--}
-eigHg :: Matrix (Complex Double)-> (Vector Double, Matrix (Complex Double))
-eigHg = eigHg' . cmat
-
-eigHg' m
-    | r == 1 = (fromList [realPart $ cdat m `at` 0], singleton 1)
-    | otherwise = unsafePerformIO $ do
-        l <- createVector r
-        v <- createMatrix RowMajor r r
-        app3 c_eigH mat m vec l mat v "eigHg"
-        return (l,v)
-  where r = rows m
-foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
-
-
-{- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/:
-
-
--}
-svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdg x = if rows x >= cols x
-    then svd' (cmat x)
-    else (v, s, u) where (u,s,v) = svd' (cmat (trans x))
-
-svd' x = unsafePerformIO $ do
-    u <- createMatrix RowMajor r c
-    s <- createVector c
-    v <- createMatrix RowMajor c c
-    app4 c_svd mat x mat u vec s mat v "svdg"
-    return (u,s,v)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
-
-{- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
-
--}
-qr :: Matrix Double -> (Matrix Double, Matrix Double)
-qr = qr' . cmat
-
-qr' x = unsafePerformIO $ do
-    q <- createMatrix RowMajor r r
-    rot <- createMatrix RowMajor r c
-    app3 c_qr mat x mat q mat rot  "qr"
-    return (q,rot)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h QR" c_qr :: TMMM
-
-qrPacked :: Matrix Double -> (Matrix Double, Vector Double)
-qrPacked = qrPacked' . cmat
-
-qrPacked' x = unsafePerformIO $ do
-    qrp <- createMatrix RowMajor r c
-    tau <- createVector (min r c)
-    app3 c_qrPacked mat x mat qrp vec tau "qrUnpacked"
-    return (qrp,tau)
-  where r = rows x
-        c = cols x
-foreign import ccall "gsl-aux.h QRpacked" c_qrPacked :: TMMV
-
-unpackQR :: (Matrix Double, Vector Double) -> (Matrix Double, Matrix Double)
-unpackQR (qrp,tau) = unpackQR' (cmat qrp, tau)
-
-unpackQR' (qrp,tau) = unsafePerformIO $ do
-    q <- createMatrix RowMajor r r
-    res <- createMatrix RowMajor r c
-    app4 c_qrUnpack mat qrp vec tau mat q mat res "qrUnpack"
-    return (q,res)
-  where r = rows qrp
-        c = cols qrp
-foreign import ccall "gsl-aux.h QRunpack" c_qrUnpack :: TMVMM
-
-{- | Cholesky decomposition of a symmetric positive definite real matrix using /gsl_linalg_cholesky_decomp/.
-
-@\> chol $ (2><2) [1,2,
-                   2,9::Double]
-(2><2)
- [ 1.0,              0.0
- , 2.0, 2.23606797749979 ]@
-
--}
-cholR :: Matrix Double -> Matrix Double
-cholR = cholR' . cmat
-
-cholR' x = unsafePerformIO $ do
-    r <- createMatrix RowMajor n n
-    app2 c_cholR mat x mat r "cholR"
-    return r
-  where n = rows x
-foreign import ccall "gsl-aux.h cholR" c_cholR :: TMM
-
-cholC :: Matrix (Complex Double) -> Matrix (Complex Double)
-cholC = cholC' . cmat
-
-cholC' x = unsafePerformIO $ do
-    r <- createMatrix RowMajor n n
-    app2 c_cholC mat x mat r "cholC"
-    return r
-  where n = rows x
-foreign import ccall "gsl-aux.h cholC" c_cholC :: TCMCM
-
-
---------------------------------------------------------
-
-{- -| efficient multiplication by the inverse of a matrix (for real matrices)
--}
-luSolveR :: Matrix Double -> Matrix Double -> Matrix Double
-luSolveR a b = luSolveR' (cmat a) (cmat b)
-
-luSolveR' a b
-    | n1==n2 && n1==r = unsafePerformIO $ do
-        s <- createMatrix RowMajor r c
-        app3 c_luSolveR mat a mat b mat s "luSolveR"
-        return s
-    | otherwise = error "luSolveR of nonsquare matrix"
-  where n1 = rows a
-        n2 = cols a
-        r  = rows b
-        c  = cols b
-foreign import ccall "gsl-aux.h luSolveR" c_luSolveR ::  TMMM
-
-{- -| efficient multiplication by the inverse of a matrix (for complex matrices). 
--}
-luSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
-luSolveC a b = luSolveC' (cmat a) (cmat b)
-
-luSolveC' a b
-    | n1==n2 && n1==r = unsafePerformIO $ do
-        s <- createMatrix RowMajor r c
-        app3 c_luSolveC mat a mat b mat s "luSolveC"
-        return s
-    | otherwise = error "luSolveC of nonsquare matrix"
-  where n1 = rows a
-        n2 = cols a
-        r  = rows b
-        c  = cols b
-foreign import ccall "gsl-aux.h luSolveC" c_luSolveC ::  TCMCMCM
-
-{- | lu decomposition of real matrix (packed as a vector including l, u, the permutation and sign)
--}
-luRaux  :: Matrix Double -> Vector Double
-luRaux = luRaux' . cmat
-
-luRaux' x = unsafePerformIO $ do
-    res <- createVector (r*r+r+1)
-    app2 c_luRaux mat x vec res "luRaux"
-    return res
-  where r = rows x
-foreign import ccall "gsl-aux.h luRaux" c_luRaux :: TMV
-
-{- | lu decomposition of complex matrix (packed as a vector including l, u, the permutation and sign)
--}
-luCaux  :: Matrix (Complex Double) -> Vector (Complex Double)
-luCaux = luCaux' . cmat
-
-luCaux' x = unsafePerformIO $ do
-    res <- createVector (r*r+r+1)
-    app2 c_luCaux mat x vec res "luCaux"
-    return res
-  where r = rows x
-foreign import ccall "gsl-aux.h luCaux" c_luCaux :: TCMCV
-
-{- | The LU decomposition of a square matrix. Is based on /gsl_linalg_LU_decomp/ and  /gsl_linalg_complex_LU_decomp/ as described in <http://www.gnu.org/software/Numeric.GSL/manual/Numeric.GSL-ref_13.html#SEC223>.
-
-@\> let m = 'fromLists' [[1,2,-3],[2+3*i,-7,0],[1,-i,2*i]]
-\> let (l,u,p,s) = luR m@
-
-L is the lower triangular:
-
-@\> l
-          1.            0.  0.
-0.154-0.231i            1.  0.
-0.154-0.231i  0.624-0.522i  1.@
-
-U is the upper triangular:
-
-@\> u
-2.+3.i           -7.            0.
-    0.  3.077-1.615i           -3.
-    0.            0.  1.873+0.433i@
-
-p is a permutation:
-
-@\> p
-[1,0,2]@
-
-L \* U obtains a permuted version of the original matrix:
-
-@\> extractRows p m
-  2.+3.i   -7.   0.
-      1.    2.  -3.
-      1.  -1.i  2.i
-\ -- CPP
-\> l \<\> u
- 2.+3.i   -7.   0.
-     1.    2.  -3.
-     1.  -1.i  2.i@
-
-s is the sign of the permutation, required to obtain sign of the determinant:
-
-@\> s * product ('toList' $ 'takeDiag' u)
-(-18.0) :+ (-16.000000000000004)
-\> 'LinearAlgebra.Algorithms.det' m
-(-18.0) :+ (-16.000000000000004)@
-
- -}
-luR :: Matrix Double -> (Matrix Double, Matrix Double, [Int], Double)
-luR m = (l,u,p, fromIntegral s') where
-    r = rows m
-    v = luRaux m
-    lu = reshape r $ subVector 0 (r*r) v
-    s':p = map round . toList . subVector (r*r) (r+1) $ v
-    u = triang r r 0 1`mul` lu
-    l = (triang r r 0 0 `mul` lu) `add` ident r
-    add = liftMatrix2 $ vectorZipR Add
-    mul = liftMatrix2 $ vectorZipR Mul
-
--- | Complex version of 'luR'.
-luC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double), [Int], Complex Double)
-luC m = (l,u,p, fromIntegral s') where
-    r = rows m
-    v = luCaux m
-    lu = reshape r $ subVector 0 (r*r) v
-    s':p = map (round.realPart) . toList . subVector (r*r) (r+1) $ v
-    u = triang r r 0 1 `mul` lu
-    l = (triang r r 0 0 `mul` lu) `add` liftMatrix comp (ident r)
-    add = liftMatrix2 $ vectorZipC Add
-    mul = liftMatrix2 $ vectorZipC Mul
-
-{- auxiliary function to get triangular matrices
--}
-triang r c h v = reshape c $ fromList [el i j | i<-[0..r-1], j<-[0..c-1]]
-    where el i j = if j-i>=h then v else 1 - v
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index bd0a6bd..052cafd 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -1,11 +1,8 @@
 #include "gsl-aux.h"
 #include <gsl/gsl_blas.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_matrix.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_fft_complex.h>
-#include <gsl/gsl_eigen.h>
 #include <gsl/gsl_integration.h>
 #include <gsl/gsl_deriv.h>
 #include <gsl/gsl_poly.h>
@@ -161,47 +158,6 @@ int mapC(int code, KCVEC(x), CVEC(r)) {
 }
 
 
-/*
-int scaleR(double* alpha, KRVEC(x), RVEC(r)) {
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("scaleR");
-    KDVVIEW(x);
-    DVVIEW(r);
-    CHECK( gsl_vector_memcpy(V(r),V(x)) , MEM);
-    int res = gsl_vector_scale(V(r),*alpha);
-    CHECK(res,res);
-    OK
-}
-
-int scaleC(gsl_complex *alpha, KCVEC(x), CVEC(r)) {
-    REQUIRES(xn == rn,BAD_SIZE);
-    DEBUGMSG("scaleC");
-    //KCVVIEW(x);
-    CVVIEW(r);
-    gsl_vector_const_view vrx = gsl_vector_const_view_array((double*)xp,xn*2);
-    gsl_vector_view vrr = gsl_vector_view_array((double*)rp,rn*2);
-    CHECK(gsl_vector_memcpy(V(vrr),V(vrx)) , MEM);
-    gsl_blas_zscal(*alpha,V(r)); // void !
-    int res = 0; 
-    CHECK(res,res);
-    OK
-}
-
-int addConstantR(double offs, KRVEC(x), RVEC(r)) { 
-    REQUIRES(xn == rn,BAD_SIZE); 
-    DEBUGMSG("addConstantR");
-    KDVVIEW(x);
-    DVVIEW(r);
-    CHECK(gsl_vector_memcpy(V(r),V(x)), MEM);
-    int res = gsl_vector_add_constant(V(r),offs);
-    CHECK(res,res);
-    OK
-}
-
-*/
-
-
-
 int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) {
     int k;
     double val = *pval;
@@ -291,248 +247,6 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) {
 
 
 
-
-int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)) {
-    REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("luSolveR");
-    KDMVIEW(a);
-    KDMVIEW(b);
-    DMVIEW(r);
-    int res;
-    gsl_matrix *LU = gsl_matrix_alloc(ar,ar);
-    CHECK(!LU,MEM);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_memcpy(LU,M(a)),MEM);
-    res = gsl_linalg_LU_decomp(LU, p, &s);
-    CHECK(res,res);
-    int c;
-
-    for (c=0; c<bc; c++) {
-        gsl_vector_const_view colb = gsl_matrix_const_column (M(b), c);
-        gsl_vector_view colr = gsl_matrix_column (M(r), c);
-        res = gsl_linalg_LU_solve (LU, p, V(colb), V(colr));
-        CHECK(res,res);
-    }
-    gsl_permutation_free(p);
-    gsl_matrix_free(LU);
-    OK
-}
-
-
-int luSolveC(KCMAT(a),KCMAT(b),CMAT(r)) {
-    REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE);
-    DEBUGMSG("luSolveC");
-    KCMVIEW(a);
-    KCMVIEW(b);
-    CMVIEW(r);
-    gsl_matrix_complex *LU = gsl_matrix_complex_alloc(ar,ar);
-    CHECK(!LU,MEM);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_complex_memcpy(LU,M(a)),MEM);
-    int res;
-    res = gsl_linalg_complex_LU_decomp(LU, p, &s);
-    CHECK(res,res);
-    int c;
-    for (c=0; c<bc; c++) {
-        gsl_vector_complex_const_view colb = gsl_matrix_complex_const_column (M(b), c);
-        gsl_vector_complex_view colr = gsl_matrix_complex_column (M(r), c);
-        res = gsl_linalg_complex_LU_solve (LU, p, V(colb), V(colr));
-        CHECK(res,res);
-    }
-    gsl_permutation_free(p);
-    gsl_matrix_complex_free(LU);
-    OK
-}
-
-
-int luRaux(KRMAT(a),RVEC(b)) {
-    REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
-    DEBUGMSG("luRaux");
-    KDMVIEW(a);
-    //DVVIEW(b);
-    gsl_matrix_view LU = gsl_matrix_view_array(bp,ar,ac);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_memcpy(M(LU),M(a)),MEM);
-    gsl_linalg_LU_decomp(M(LU), p, &s);
-    bp[ar*ar] = s;
-    int k;
-    for (k=0; k<ar; k++) {
-        bp[ar*ar+k+1] = gsl_permutation_get(p,k);
-    }
-    gsl_permutation_free(p);
-    OK
-}
-
-int luCaux(KCMAT(a),CVEC(b)) {
-    REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE);
-    DEBUGMSG("luCaux");
-    KCMVIEW(a);
-    //DVVIEW(b);
-    gsl_matrix_complex_view LU = gsl_matrix_complex_view_array((double*)bp,ar,ac);
-    int s;
-    gsl_permutation * p = gsl_permutation_alloc (ar);
-    CHECK(!p,MEM);
-    CHECK(gsl_matrix_complex_memcpy(M(LU),M(a)),MEM);
-    int res;
-    res = gsl_linalg_complex_LU_decomp(M(LU), p, &s);
-    CHECK(res,res);
-    ((double*)bp)[2*ar*ar] = s;
-    ((double*)bp)[2*ar*ar+1] = 0;
-    int k;
-    for (k=0; k<ar; k++) {
-        ((double*)bp)[2*ar*ar+2*k+2] = gsl_permutation_get(p,k);
-        ((double*)bp)[2*ar*ar+2*k+2+1] = 0;
-    }
-    gsl_permutation_free(p);
-    OK
-}
-
-int svd(KRMAT(a),RMAT(u), RVEC(s),RMAT(v)) {
-    REQUIRES(ar==ur && ac==uc && ac==sn && ac==vr && ac==vc,BAD_SIZE);
-    DEBUGMSG("svd");
-    KDMVIEW(a);
-    DMVIEW(u);
-    DVVIEW(s);
-    DMVIEW(v);
-    gsl_vector  *workv = gsl_vector_alloc(ac);
-    CHECK(!workv,MEM);
-    gsl_matrix  *workm = gsl_matrix_alloc(ac,ac);
-    CHECK(!workm,MEM);
-    CHECK(gsl_matrix_memcpy(M(u),M(a)),MEM);
-    // int res = gsl_linalg_SV_decomp_jacobi (&U.matrix, &V.matrix, &S.vector);
-    // doesn't work 
-    //int res = gsl_linalg_SV_decomp (&U.matrix, &V.matrix, &S.vector, workv);
-    int res = gsl_linalg_SV_decomp_mod (M(u), workm, M(v), V(s), workv);
-    CHECK(res,res);
-    //gsl_matrix_transpose(M(v));
-    gsl_vector_free(workv);
-    gsl_matrix_free(workm);
-    OK
-}
-
-
-// for real symmetric matrices
-int eigensystemR(KRMAT(x),RVEC(l),RMAT(v)) {
-    REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
-    DEBUGMSG("eigensystemR (gsl_eigen_symmv)");
-    KDMVIEW(x);
-    DVVIEW(l);
-    DMVIEW(v);
-    gsl_matrix *XC = gsl_matrix_alloc(xr,xr);
-    gsl_matrix_memcpy(XC,M(x)); // needed because the argument is destroyed
-                                // many thanks to Nico Mahlo for the bug report 
-    gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (xc);
-    int res = gsl_eigen_symmv (XC, V(l), M(v), w);
-    CHECK(res,res);
-    gsl_eigen_symmv_free (w);
-    gsl_matrix_free(XC);
-    gsl_eigen_symmv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
-    OK
-}
-
-// for hermitian matrices
-int eigensystemC(KCMAT(x),RVEC(l),CMAT(v)) {
-    REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE);
-    DEBUGMSG("eigensystemC");
-    KCMVIEW(x);
-    DVVIEW(l);
-    CMVIEW(v);
-    gsl_matrix_complex *XC = gsl_matrix_complex_alloc(xr,xr);
-    gsl_matrix_complex_memcpy(XC,M(x)); // again needed because the argument is destroyed
-    gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc (xc);
-    int res = gsl_eigen_hermv (XC, V(l), M(v), w);
-    CHECK(res,res);
-    gsl_eigen_hermv_free (w);
-    gsl_matrix_complex_free(XC);
-    gsl_eigen_hermv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC);
-    OK
-}
-
-int QR(KRMAT(x),RMAT(q),RMAT(r)) {
-    REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QR");
-    KDMVIEW(x);
-    DMVIEW(q);
-    DMVIEW(r);
-    gsl_matrix * a = gsl_matrix_alloc(xr,xc);
-    gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
-    gsl_matrix_memcpy(a,M(x));
-    int res = gsl_linalg_QR_decomp(a,tau);
-    CHECK(res,res);
-    gsl_linalg_QR_unpack(a,tau,M(q),M(r));
-    gsl_vector_free(tau);
-    gsl_matrix_free(a);
-    OK
-}
-
-int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau)) {
-    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QRpacked");
-    KDMVIEW(x);
-    DMVIEW(qr);
-    DVVIEW(tau);
-    //gsl_matrix * a = gsl_matrix_alloc(xr,xc);
-    //gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc));
-    gsl_matrix_memcpy(M(qr),M(x));
-    int res = gsl_linalg_QR_decomp(M(qr),V(tau));
-    CHECK(res,res);
-    OK
-}
-
-
-int QRunpack(KRMAT(xqr),KRVEC(tau),RMAT(q),RMAT(r)) {
-    //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE);
-    DEBUGMSG("QRunpack");
-    KDMVIEW(xqr);
-    KDVVIEW(tau);
-    DMVIEW(q);
-    DMVIEW(r);
-    gsl_linalg_QR_unpack(M(xqr),V(tau),M(q),M(r));
-    OK
-}
-
-
-int cholR(KRMAT(x),RMAT(l)) {
-    REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
-    DEBUGMSG("cholR");
-    KDMVIEW(x);
-    DMVIEW(l);
-    gsl_matrix_memcpy(M(l),M(x));
-    int res = gsl_linalg_cholesky_decomp(M(l));
-    CHECK(res,res);
-    int r,c;
-    for (r=0; r<xr-1; r++) {
-        for(c=r+1; c<xc; c++) {
-            lp[r*lc+c] = 0.;
-        }
-    }
-    OK
-}
-
-int cholC(KCMAT(x),CMAT(l)) {
-    REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE);
-    DEBUGMSG("cholC");
-    KCMVIEW(x);
-    CMVIEW(l);
-    gsl_matrix_complex_memcpy(M(l),M(x));
-    int res = 0; // gsl_linalg_complex_cholesky_decomp(M(l));
-    CHECK(res,res);
-    gsl_complex zero = {0.,0.};
-    int r,c;
-    for (r=0; r<xr-1; r++) {
-        for(c=r+1; c<xc; c++) {
-            lp[r*lc+c] = zero;
-        }
-    }
-    OK
-}
-
 int fft(int code, KCVEC(X), CVEC(R)) {
     REQUIRES(Xn == Rn,BAD_SIZE);
     DEBUGMSG("fft");
diff --git a/lib/Numeric/GSL/gsl-aux.h b/lib/Numeric/GSL/gsl-aux.h
index eee15e7..cd17ef0 100644
--- a/lib/Numeric/GSL/gsl-aux.h
+++ b/lib/Numeric/GSL/gsl-aux.h
@@ -26,25 +26,6 @@ int zipR(int code, KRVEC(a), KRVEC(b), RVEC(r));
 int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r));
 
 
-int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)); 
-int luSolveC(KCMAT(a),KCMAT(b),CMAT(r));
-int luRaux(KRMAT(a),RVEC(b));
-int luCaux(KCMAT(a),CVEC(b));
-
-int svd(KRMAT(x),RMAT(u), RVEC(s),RMAT(v));
-
-int eigensystemR(KRMAT(x),RVEC(l),RMAT(v));
-int eigensystemC(KCMAT(x),RVEC(l),CMAT(v));
-
-int QR(KRMAT(x),RMAT(q),RMAT(r));
-
-int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau));
-int QRunpack(KRMAT(qr),KRVEC(tau),RMAT(q),RMAT(r));
-
-int cholR(KRMAT(x),RMAT(l));
-
-int cholC(KCMAT(x),CMAT(l));
-
 int fft(int code, KCVEC(a), CVEC(b));
 
 int integrate_qng(double f(double, void*), double a, double b, double prec,
-- 
cgit v1.2.3