From 561a6c0e21bb77c21114ccbbd86d3af5ddb5a3f1 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 12:18:56 +0200 Subject: Conversion, LAPACK -> base --- packages/base/hmatrix-base.cabal | 4 +- packages/base/src/Data/Packed/Development.hs | 2 +- packages/base/src/Data/Packed/Vector.hs | 13 +- packages/base/src/Numeric/Conversion.hs | 91 ++ packages/base/src/Numeric/LAPACK.hs | 552 ++++++++ packages/hmatrix/src/Numeric/Conversion.hs | 91 -- .../hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs | 555 -------- .../src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 1489 -------------------- .../src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h | 60 - 9 files changed, 659 insertions(+), 2198 deletions(-) create mode 100644 packages/base/src/Numeric/Conversion.hs create mode 100644 packages/base/src/Numeric/LAPACK.hs delete mode 100644 packages/hmatrix/src/Numeric/Conversion.hs delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c delete mode 100644 packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 96f90e2..0413d4a 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal @@ -33,7 +33,9 @@ library Data.Packed.Matrix, Data.Packed.Foreign, Data.Packed.ST, - Data.Packed.Development + Data.Packed.Development, + Numeric.Conversion + Numeric.LAPACK other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs index 777b6c5..9350acb 100644 --- a/packages/base/src/Data/Packed/Development.hs +++ b/packages/base/src/Data/Packed/Development.hs @@ -24,7 +24,7 @@ module Data.Packed.Development ( unsafeFromForeignPtr, unsafeToForeignPtr, check, (//), - at', atM' + at', atM', fi, table ) where import Data.Packed.Internal diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs index b5a4318..653a257 100644 --- a/packages/base/src/Data/Packed/Vector.hs +++ b/packages/base/src/Data/Packed/Vector.hs @@ -18,7 +18,7 @@ module Data.Packed.Vector ( Vector, - fromList, (|>), toList, buildVector, + fromList, (|>), toList, buildVector, constant, dim, (@>), subVector, takesV, vjoin, join, mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, @@ -27,6 +27,7 @@ module Data.Packed.Vector ( ) where import Data.Packed.Internal.Vector +import Data.Packed.Internal.Matrix import Foreign.Storable ------------------------------------------------------------------- @@ -94,3 +95,13 @@ unzipVector = unzipVectorWith id join :: Storable t => [Vector t] -> Vector t join = vjoin +{- | creates a vector with a given number of equal components: + +@> constant 2 7 +7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ +-} +constant :: Element a => a -> Int -> Vector a +-- constant x n = runSTVector (newVector x n) +constant = constantD-- about 2x faster + + diff --git a/packages/base/src/Numeric/Conversion.hs b/packages/base/src/Numeric/Conversion.hs new file mode 100644 index 0000000..a1f9385 --- /dev/null +++ b/packages/base/src/Numeric/Conversion.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.Conversion +-- Copyright : (c) Alberto Ruiz 2010 +-- License : BSD3 +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Conversion routines +-- +----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} + + +module Numeric.Conversion ( + Complexable(..), RealElement, + module Data.Complex +) where + +import Data.Packed.Internal.Vector +import Data.Packed.Internal.Matrix +import Data.Complex +import Control.Arrow((***)) + +------------------------------------------------------------------- + +-- | Supported single-double precision type pairs +class (Element s, Element d) => Precision s d | s -> d, d -> s where + double2FloatG :: Vector d -> Vector s + float2DoubleG :: Vector s -> Vector d + +instance Precision Float Double where + double2FloatG = double2FloatV + float2DoubleG = float2DoubleV + +instance Precision (Complex Float) (Complex Double) where + double2FloatG = asComplex . double2FloatV . asReal + float2DoubleG = asComplex . float2DoubleV . asReal + +-- | Supported real types +class (Element t, Element (Complex t), RealFloat t +-- , RealOf t ~ t, RealOf (Complex t) ~ t + ) + => RealElement t + +instance RealElement Double +instance RealElement Float + + +-- | Structures that may contain complex numbers +class Complexable c where + toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) + fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) + comp' :: (RealElement e) => c e -> c (Complex e) + single' :: Precision a b => c b -> c a + double' :: Precision a b => c a -> c b + + +instance Complexable Vector where + toComplex' = toComplexV + fromComplex' = fromComplexV + comp' v = toComplex' (v,constantD 0 (dim v)) + single' = double2FloatG + double' = float2DoubleG + + +-- | creates a complex vector from vectors with real and imaginary parts +toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) +toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] + +-- | the inverse of 'toComplex' +fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) +fromComplexV z = (r,i) where + [r,i] = toColumns $ reshape 2 $ asReal z + + +instance Complexable Matrix where + toComplex' = uncurry $ liftMatrix2 $ curry toComplex' + fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z + where c = cols z + comp' = liftMatrix comp' + single' = liftMatrix single' + double' = liftMatrix double' + diff --git a/packages/base/src/Numeric/LAPACK.hs b/packages/base/src/Numeric/LAPACK.hs new file mode 100644 index 0000000..08cd759 --- /dev/null +++ b/packages/base/src/Numeric/LAPACK.hs @@ -0,0 +1,552 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.LAPACK +-- Copyright : (c) Alberto Ruiz 2006-14 +-- License : BSD3 +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Functional interface to selected LAPACK functions (). +-- +----------------------------------------------------------------------------- + +module Numeric.LAPACK ( + -- * Matrix product + multiplyR, multiplyC, multiplyF, multiplyQ, + -- * Linear systems + linearSolveR, linearSolveC, + lusR, lusC, + cholSolveR, cholSolveC, + linearSolveLSR, linearSolveLSC, + linearSolveSVDR, linearSolveSVDC, + -- * SVD + svR, svRd, svC, svCd, + svdR, svdRd, svdC, svdCd, + thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, + rightSVR, rightSVC, leftSVR, leftSVC, + -- * Eigensystems + eigR, eigC, eigS, eigS', eigH, eigH', + eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, + -- * LU + luR, luC, + -- * Cholesky + cholS, cholH, mbCholS, mbCholH, + -- * QR + qrR, qrC, qrgrR, qrgrC, + -- * Hessenberg + hessR, hessC, + -- * Schur + schurR, schurC +) where + +import Data.Packed.Development +import Data.Packed +import Data.Packed.Internal +import Numeric.Conversion + +import Foreign.Ptr(nullPtr) +import Foreign.C.Types +import Control.Monad(when) +import System.IO.Unsafe(unsafePerformIO) + +----------------------------------------------------------------------------------- + +foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM +foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM +foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM +foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM + +isT Matrix{order = ColumnMajor} = 0 +isT Matrix{order = RowMajor} = 1 + +tt x@Matrix{order = ColumnMajor} = x +tt x@Matrix{order = RowMajor} = trans x + +multiplyAux f st a b = unsafePerformIO $ do + when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ + show (rows a,cols a) ++ " x " ++ show (rows b, cols b) + s <- createMatrix ColumnMajor (rows a) (cols b) + app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st + return s + +-- | Matrix product based on BLAS's /dgemm/. +multiplyR :: Matrix Double -> Matrix Double -> Matrix Double +multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b + +-- | Matrix product based on BLAS's /zgemm/. +multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +multiplyC a b = multiplyAux zgemmc "zgemmc" a b + +-- | Matrix product based on BLAS's /sgemm/. +multiplyF :: Matrix Float -> Matrix Float -> Matrix Float +multiplyF a b = multiplyAux sgemmc "sgemmc" a b + +-- | Matrix product based on BLAS's /cgemm/. +multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) +multiplyQ a b = multiplyAux cgemmc "cgemmc" a b + +----------------------------------------------------------------------------- +foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM +foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM +foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM +foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM + +-- | Full SVD of a real matrix using LAPACK's /dgesvd/. +svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdR = svdAux dgesvd "svdR" . fmat + +-- | Full SVD of a real matrix using LAPACK's /dgesdd/. +svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +svdRd = svdAux dgesdd "svdRdd" . fmat + +-- | Full SVD of a complex matrix using LAPACK's /zgesvd/. +svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +svdC = svdAux zgesvd "svdC" . fmat + +-- | Full SVD of a complex matrix using LAPACK's /zgesdd/. +svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +svdCd = svdAux zgesdd "svdCdd" . fmat + +svdAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r r + s <- createVector (min r c) + v <- createMatrix ColumnMajor c c + app4 f mat x mat u vec s mat v st + return (u,s,trans v) + where r = rows x + c = cols x + + +-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. +thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. +thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat + +-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. +thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) +thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat + +-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. +thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) +thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat + +thinSVDAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r q + s <- createVector q + v <- createMatrix ColumnMajor q c + app4 f mat x mat u vec s mat v st + return (u,s,trans v) + where r = rows x + c = cols x + q = min r c + + +-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. +svR :: Matrix Double -> Vector Double +svR = svAux dgesvd "svR" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. +svC :: Matrix (Complex Double) -> Vector Double +svC = svAux zgesvd "svC" . fmat + +-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. +svRd :: Matrix Double -> Vector Double +svRd = svAux dgesdd "svRd" . fmat + +-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. +svCd :: Matrix (Complex Double) -> Vector Double +svCd = svAux zgesdd "svCd" . fmat + +svAux f st x = unsafePerformIO $ do + s <- createVector q + app2 g mat x vec s st + return s + where r = rows x + c = cols x + q = min r c + g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr + + +-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVR :: Matrix Double -> (Vector Double, Matrix Double) +rightSVR = rightSVAux dgesvd "rightSVR" . fmat + +-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. +rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +rightSVC = rightSVAux zgesvd "rightSVC" . fmat + +rightSVAux f st x = unsafePerformIO $ do + s <- createVector q + v <- createMatrix ColumnMajor c c + app3 g mat x vec s mat v st + return (s,trans v) + where r = rows x + c = cols x + q = min r c + g ra ca pa = f ra ca pa 0 0 nullPtr + + +-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVR :: Matrix Double -> (Matrix Double, Vector Double) +leftSVR = leftSVAux dgesvd "leftSVR" . fmat + +-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. +leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) +leftSVC = leftSVAux zgesvd "leftSVC" . fmat + +leftSVAux f st x = unsafePerformIO $ do + u <- createMatrix ColumnMajor r r + s <- createVector q + app3 g mat x mat u vec s st + return (u,s) + where r = rows x + c = cols x + q = min r c + g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr + +----------------------------------------------------------------------------- + +foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM +foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM +foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM +foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM + +eigAux f st m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 g mat m vec l mat v st + return (l,v) + where r = rows m + g ra ca pa = f ra ca pa 0 0 nullPtr + + +-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. +eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) +eigC = eigAux zgeev "eigC" . fmat + +eigOnlyAux f st m = unsafePerformIO $ do + l <- createVector r + app2 g mat m vec l st + return l + where r = rows m + g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr + +-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. +-- The eigenvalues are not sorted. +eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) +eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat + +-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. +-- The eigenvectors are the columns of v. The eigenvalues are not sorted. +eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) +eigR m = (s', v'') + where (s,v) = eigRaux (fmat m) + s' = fixeig1 s + v' = toRows $ trans v + v'' = fromColumns $ fixeig (toList s') v' + +eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) +eigRaux m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 g mat m vec l mat v "eigR" + return (l,v) + where r = rows m + g ra ca pa = dgeev ra ca pa 0 0 nullPtr + +fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) + where r = dim s + +fixeig [] _ = [] +fixeig [_] [v] = [comp' v] +fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) + | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1, mapVector negate v2) : fixeig r vs + | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) +fixeig _ _ = error "fixeig with impossible inputs" + + +-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. +-- The eigenvalues are not sorted. +eigOnlyR :: Matrix Double -> Vector (Complex Double) +eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat + + +----------------------------------------------------------------------------- + +eigSHAux f st m = unsafePerformIO $ do + l <- createVector r + v <- createMatrix ColumnMajor r r + app3 f mat m vec l mat v st + return (l,v) + where r = rows m + +-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. +-- The eigenvectors are the columns of v. +-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). +eigS :: Matrix Double -> (Vector Double, Matrix Double) +eigS m = (s', fliprl v) + where (s,v) = eigS' (fmat m) + s' = fromList . reverse . toList $ s + +-- | 'eigS' in ascending order +eigS' :: Matrix Double -> (Vector Double, Matrix Double) +eigS' = eigSHAux (dsyev 1) "eigS'" . fmat + +-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. +-- The eigenvectors are the columns of v. +-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). +eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +eigH m = (s', fliprl v) + where (s,v) = eigH' (fmat m) + s' = fromList . reverse . toList $ s + +-- | 'eigH' in ascending order +eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) +eigH' = eigSHAux (zheev 1) "eigH'" . fmat + + +-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. +-- The eigenvalues are sorted in descending order. +eigOnlyS :: Matrix Double -> Vector Double +eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat + +-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. +-- The eigenvalues are sorted in descending order. +eigOnlyH :: Matrix (Complex Double) -> Vector Double +eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat + +vrev = flatten . flipud . reshape 1 + +----------------------------------------------------------------------------- +foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM +foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM +foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM +foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM + +linearSolveSQAux f st a b + | n1==n2 && n1==r = unsafePerformIO $ do + s <- createMatrix ColumnMajor r c + app3 f mat a mat b mat s st + return s + | otherwise = error $ st ++ " of nonsquare matrix" + where n1 = rows a + n2 = cols a + r = rows b + c = cols b + +-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. +linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double +linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) + +-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. +linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) + + +-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. +cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double +cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b) + +-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. +cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM +foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM +foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM +foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM + +linearSolveAux f st a b = unsafePerformIO $ do + r <- createMatrix ColumnMajor (max m n) nrhs + app3 f mat a mat b mat r st + return r + where m = rows a + n = cols a + nrhs = cols b + +-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. +linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double +linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) + +-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. +linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) +linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) + +-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +linearSolveSVDR :: Maybe Double -- ^ rcond + -> Matrix Double -- ^ coefficient matrix + -> Matrix Double -- ^ right hand sides (as columns) + -> Matrix Double -- ^ solution vectors (as columns) +linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) +linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) + +-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. +linearSolveSVDC :: Maybe Double -- ^ rcond + -> Matrix (Complex Double) -- ^ coefficient matrix + -> Matrix (Complex Double) -- ^ right hand sides (as columns) + -> Matrix (Complex Double) -- ^ solution vectors (as columns) +linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ + linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) +linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM +foreign import ccall unsafe "chol_l_S" dpotrf :: TMM + +cholAux f st a = do + r <- createMatrix ColumnMajor n n + app2 f mat a mat r st + return r + where n = rows a + +-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. +cholH :: Matrix (Complex Double) -> Matrix (Complex Double) +cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat + +-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. +cholS :: Matrix Double -> Matrix Double +cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat + +-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version). +mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) +mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat + +-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version). +mbCholS :: Matrix Double -> Maybe (Matrix Double) +mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM +foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM + +-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. +qrR :: Matrix Double -> (Matrix Double, Vector Double) +qrR = qrAux dgeqr2 "qrR" . fmat + +-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. +qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) +qrC = qrAux zgeqr2 "qrC" . fmat + +qrAux f st a = unsafePerformIO $ do + r <- createMatrix ColumnMajor m n + tau <- createVector mn + app3 f mat a vec tau mat r st + return (r,tau) + where + m = rows a + n = cols a + mn = min m n + +foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM +foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM + +-- | build rotation from reflectors +qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double +qrgrR = qrgrAux dorgqr "qrgrR" +-- | build rotation from reflectors +qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) +qrgrC = qrgrAux zungqr "qrgrC" + +qrgrAux f st n (a, tau) = unsafePerformIO $ do + res <- createMatrix ColumnMajor (rows a) n + app3 f mat (fmat a) vec (subVector 0 n tau') mat res st + return res + where + tau' = vjoin [tau, constantD 0 n] + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM +foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM + +-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. +hessR :: Matrix Double -> (Matrix Double, Vector Double) +hessR = hessAux dgehrd "hessR" . fmat + +-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. +hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) +hessC = hessAux zgehrd "hessC" . fmat + +hessAux f st a = unsafePerformIO $ do + r <- createMatrix ColumnMajor m n + tau <- createVector (mn-1) + app3 f mat a vec tau mat r st + return (r,tau) + where m = rows a + n = cols a + mn = min m n + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "schur_l_R" dgees :: TMMM +foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM + +-- | Schur factorization of a square real matrix, using LAPACK's /dgees/. +schurR :: Matrix Double -> (Matrix Double, Matrix Double) +schurR = schurAux dgees "schurR" . fmat + +-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. +schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) +schurC = schurAux zgees "schurC" . fmat + +schurAux f st a = unsafePerformIO $ do + u <- createMatrix ColumnMajor n n + s <- createMatrix ColumnMajor n n + app3 f mat a mat u mat s st + return (u,s) + where n = rows a + +----------------------------------------------------------------------------------- +foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM +foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM + +-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. +luR :: Matrix Double -> (Matrix Double, [Int]) +luR = luAux dgetrf "luR" . fmat + +-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. +luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) +luC = luAux zgetrf "luC" . fmat + +luAux f st a = unsafePerformIO $ do + lu <- createMatrix ColumnMajor n m + piv <- createVector (min n m) + app3 f mat a vec piv mat lu st + return (lu, map (pred.round) (toList piv)) + where n = rows a + m = cols a + +----------------------------------------------------------------------------------- +type TW a = CInt -> PD -> a +type TQ a = CInt -> CInt -> PC -> a + +foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM +foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) + +-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. +lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double +lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) + +-- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. +lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) +lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) + +lusAux f st a piv b + | n1==n2 && n2==n =unsafePerformIO $ do + x <- createMatrix ColumnMajor n m + app4 f mat a vec piv' mat b mat x st + return x + | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" + where n1 = rows a + n2 = cols a + n = rows b + m = cols b + piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double + diff --git a/packages/hmatrix/src/Numeric/Conversion.hs b/packages/hmatrix/src/Numeric/Conversion.hs deleted file mode 100644 index 8941451..0000000 --- a/packages/hmatrix/src/Numeric/Conversion.hs +++ /dev/null @@ -1,91 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE UndecidableInstances #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Conversion --- Copyright : (c) Alberto Ruiz 2010 --- License : GPL-style --- --- Maintainer : Alberto Ruiz --- Stability : provisional --- Portability : portable --- --- Conversion routines --- ------------------------------------------------------------------------------ - -module Numeric.Conversion ( - Complexable(..), RealElement, - module Data.Complex -) where - -import Data.Packed.Internal.Vector -import Data.Packed.Internal.Matrix -import Data.Complex -import Control.Arrow((***)) - -------------------------------------------------------------------- - --- | Supported single-double precision type pairs -class (Element s, Element d) => Precision s d | s -> d, d -> s where - double2FloatG :: Vector d -> Vector s - float2DoubleG :: Vector s -> Vector d - -instance Precision Float Double where - double2FloatG = double2FloatV - float2DoubleG = float2DoubleV - -instance Precision (Complex Float) (Complex Double) where - double2FloatG = asComplex . double2FloatV . asReal - float2DoubleG = asComplex . float2DoubleV . asReal - --- | Supported real types -class (Element t, Element (Complex t), RealFloat t --- , RealOf t ~ t, RealOf (Complex t) ~ t - ) - => RealElement t - -instance RealElement Double -instance RealElement Float - - --- | Structures that may contain complex numbers -class Complexable c where - toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) - fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) - comp' :: (RealElement e) => c e -> c (Complex e) - single' :: Precision a b => c b -> c a - double' :: Precision a b => c a -> c b - - -instance Complexable Vector where - toComplex' = toComplexV - fromComplex' = fromComplexV - comp' v = toComplex' (v,constantD 0 (dim v)) - single' = double2FloatG - double' = float2DoubleG - - --- | creates a complex vector from vectors with real and imaginary parts -toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) -toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] - --- | the inverse of 'toComplex' -fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) -fromComplexV z = (r,i) where - [r,i] = toColumns $ reshape 2 $ asReal z - - -instance Complexable Matrix where - toComplex' = uncurry $ liftMatrix2 $ curry toComplex' - fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z - where c = cols z - comp' = liftMatrix comp' - single' = liftMatrix single' - double' = liftMatrix double' - diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs deleted file mode 100644 index 11394a6..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK.hs +++ /dev/null @@ -1,555 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Numeric.LinearAlgebra.LAPACK --- Copyright : (c) Alberto Ruiz 2006-7 --- License : GPL-style --- --- Maintainer : Alberto Ruiz (aruiz at um dot es) --- Stability : provisional --- Portability : portable (uses FFI) --- --- Functional interface to selected LAPACK functions (). --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Numeric.LinearAlgebra.LAPACK ( - -- * Matrix product - multiplyR, multiplyC, multiplyF, multiplyQ, - -- * Linear systems - linearSolveR, linearSolveC, - lusR, lusC, - cholSolveR, cholSolveC, - linearSolveLSR, linearSolveLSC, - linearSolveSVDR, linearSolveSVDC, - -- * SVD - svR, svRd, svC, svCd, - svdR, svdRd, svdC, svdCd, - thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, - rightSVR, rightSVC, leftSVR, leftSVC, - -- * Eigensystems - eigR, eigC, eigS, eigS', eigH, eigH', - eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, - -- * LU - luR, luC, - -- * Cholesky - cholS, cholH, mbCholS, mbCholH, - -- * QR - qrR, qrC, qrgrR, qrgrC, - -- * Hessenberg - hessR, hessC, - -- * Schur - schurR, schurC -) where - -import Data.Packed.Internal -import Data.Packed.Matrix -import Numeric.Conversion -import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) - -import Foreign.Ptr(nullPtr) -import Foreign.C.Types -import Control.Monad(when) -import System.IO.Unsafe(unsafePerformIO) - ------------------------------------------------------------------------------------ - -foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM -foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM -foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM -foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM - -isT Matrix{order = ColumnMajor} = 0 -isT Matrix{order = RowMajor} = 1 - -tt x@Matrix{order = ColumnMajor} = x -tt x@Matrix{order = RowMajor} = trans x - -multiplyAux f st a b = unsafePerformIO $ do - when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ - show (rows a,cols a) ++ " x " ++ show (rows b, cols b) - s <- createMatrix ColumnMajor (rows a) (cols b) - app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st - return s - --- | Matrix product based on BLAS's /dgemm/. -multiplyR :: Matrix Double -> Matrix Double -> Matrix Double -multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b - --- | Matrix product based on BLAS's /zgemm/. -multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) -multiplyC a b = multiplyAux zgemmc "zgemmc" a b - --- | Matrix product based on BLAS's /sgemm/. -multiplyF :: Matrix Float -> Matrix Float -> Matrix Float -multiplyF a b = multiplyAux sgemmc "sgemmc" a b - --- | Matrix product based on BLAS's /cgemm/. -multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) -multiplyQ a b = multiplyAux cgemmc "cgemmc" a b - ------------------------------------------------------------------------------ -foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM -foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM -foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM -foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM - --- | Full SVD of a real matrix using LAPACK's /dgesvd/. -svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) -svdR = svdAux dgesvd "svdR" . fmat - --- | Full SVD of a real matrix using LAPACK's /dgesdd/. -svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) -svdRd = svdAux dgesdd "svdRdd" . fmat - --- | Full SVD of a complex matrix using LAPACK's /zgesvd/. -svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) -svdC = svdAux zgesvd "svdC" . fmat - --- | Full SVD of a complex matrix using LAPACK's /zgesdd/. -svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) -svdCd = svdAux zgesdd "svdCdd" . fmat - -svdAux f st x = unsafePerformIO $ do - u <- createMatrix ColumnMajor r r - s <- createVector (min r c) - v <- createMatrix ColumnMajor c c - app4 f mat x mat u vec s mat v st - return (u,s,trans v) - where r = rows x - c = cols x - - --- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. -thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) -thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat - --- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. -thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) -thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat - --- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. -thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) -thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat - --- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. -thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) -thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat - -thinSVDAux f st x = unsafePerformIO $ do - u <- createMatrix ColumnMajor r q - s <- createVector q - v <- createMatrix ColumnMajor q c - app4 f mat x mat u vec s mat v st - return (u,s,trans v) - where r = rows x - c = cols x - q = min r c - - --- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. -svR :: Matrix Double -> Vector Double -svR = svAux dgesvd "svR" . fmat - --- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. -svC :: Matrix (Complex Double) -> Vector Double -svC = svAux zgesvd "svC" . fmat - --- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. -svRd :: Matrix Double -> Vector Double -svRd = svAux dgesdd "svRd" . fmat - --- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. -svCd :: Matrix (Complex Double) -> Vector Double -svCd = svAux zgesdd "svCd" . fmat - -svAux f st x = unsafePerformIO $ do - s <- createVector q - app2 g mat x vec s st - return s - where r = rows x - c = cols x - q = min r c - g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr - - --- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. -rightSVR :: Matrix Double -> (Vector Double, Matrix Double) -rightSVR = rightSVAux dgesvd "rightSVR" . fmat - --- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. -rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) -rightSVC = rightSVAux zgesvd "rightSVC" . fmat - -rightSVAux f st x = unsafePerformIO $ do - s <- createVector q - v <- createMatrix ColumnMajor c c - app3 g mat x vec s mat v st - return (s,trans v) - where r = rows x - c = cols x - q = min r c - g ra ca pa = f ra ca pa 0 0 nullPtr - - --- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. -leftSVR :: Matrix Double -> (Matrix Double, Vector Double) -leftSVR = leftSVAux dgesvd "leftSVR" . fmat - --- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. -leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) -leftSVC = leftSVAux zgesvd "leftSVC" . fmat - -leftSVAux f st x = unsafePerformIO $ do - u <- createMatrix ColumnMajor r r - s <- createVector q - app3 g mat x mat u vec s st - return (u,s) - where r = rows x - c = cols x - q = min r c - g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr - ------------------------------------------------------------------------------ - -foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM -foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM -foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM -foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM - -eigAux f st m = unsafePerformIO $ do - l <- createVector r - v <- createMatrix ColumnMajor r r - app3 g mat m vec l mat v st - return (l,v) - where r = rows m - g ra ca pa = f ra ca pa 0 0 nullPtr - - --- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. --- The eigenvectors are the columns of v. The eigenvalues are not sorted. -eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) -eigC = eigAux zgeev "eigC" . fmat - -eigOnlyAux f st m = unsafePerformIO $ do - l <- createVector r - app2 g mat m vec l st - return l - where r = rows m - g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr - --- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. --- The eigenvalues are not sorted. -eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) -eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat - --- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. --- The eigenvectors are the columns of v. The eigenvalues are not sorted. -eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) -eigR m = (s', v'') - where (s,v) = eigRaux (fmat m) - s' = fixeig1 s - v' = toRows $ trans v - v'' = fromColumns $ fixeig (toList s') v' - -eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) -eigRaux m = unsafePerformIO $ do - l <- createVector r - v <- createMatrix ColumnMajor r r - app3 g mat m vec l mat v "eigR" - return (l,v) - where r = rows m - g ra ca pa = dgeev ra ca pa 0 0 nullPtr - -fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) - where r = dim s - -fixeig [] _ = [] -fixeig [_] [v] = [comp' v] -fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) - | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs - | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) - where scale = vectorMapValR Scale -fixeig _ _ = error "fixeig with impossible inputs" - - --- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. --- The eigenvalues are not sorted. -eigOnlyR :: Matrix Double -> Vector (Complex Double) -eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat - - ------------------------------------------------------------------------------ - -eigSHAux f st m = unsafePerformIO $ do - l <- createVector r - v <- createMatrix ColumnMajor r r - app3 f mat m vec l mat v st - return (l,v) - where r = rows m - --- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. --- The eigenvectors are the columns of v. --- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). -eigS :: Matrix Double -> (Vector Double, Matrix Double) -eigS m = (s', fliprl v) - where (s,v) = eigS' (fmat m) - s' = fromList . reverse . toList $ s - --- | 'eigS' in ascending order -eigS' :: Matrix Double -> (Vector Double, Matrix Double) -eigS' = eigSHAux (dsyev 1) "eigS'" . fmat - --- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. --- The eigenvectors are the columns of v. --- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). -eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) -eigH m = (s', fliprl v) - where (s,v) = eigH' (fmat m) - s' = fromList . reverse . toList $ s - --- | 'eigH' in ascending order -eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) -eigH' = eigSHAux (zheev 1) "eigH'" . fmat - - --- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. --- The eigenvalues are sorted in descending order. -eigOnlyS :: Matrix Double -> Vector Double -eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat - --- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. --- The eigenvalues are sorted in descending order. -eigOnlyH :: Matrix (Complex Double) -> Vector Double -eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat - -vrev = flatten . flipud . reshape 1 - ------------------------------------------------------------------------------ -foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM -foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM -foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM -foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM - -linearSolveSQAux f st a b - | n1==n2 && n1==r = unsafePerformIO $ do - s <- createMatrix ColumnMajor r c - app3 f mat a mat b mat s st - return s - | otherwise = error $ st ++ " of nonsquare matrix" - where n1 = rows a - n2 = cols a - r = rows b - c = cols b - --- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. -linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double -linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) - --- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. -linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) -linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) - - --- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. -cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double -cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b) - --- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. -cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) -cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM -foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM -foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM -foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM - -linearSolveAux f st a b = unsafePerformIO $ do - r <- createMatrix ColumnMajor (max m n) nrhs - app3 f mat a mat b mat r st - return r - where m = rows a - n = cols a - nrhs = cols b - --- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. -linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double -linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ - linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) - --- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. -linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) -linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ - linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) - --- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. -linearSolveSVDR :: Maybe Double -- ^ rcond - -> Matrix Double -- ^ coefficient matrix - -> Matrix Double -- ^ right hand sides (as columns) - -> Matrix Double -- ^ solution vectors (as columns) -linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ - linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) -linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) - --- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. -linearSolveSVDC :: Maybe Double -- ^ rcond - -> Matrix (Complex Double) -- ^ coefficient matrix - -> Matrix (Complex Double) -- ^ right hand sides (as columns) - -> Matrix (Complex Double) -- ^ solution vectors (as columns) -linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ - linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) -linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM -foreign import ccall unsafe "chol_l_S" dpotrf :: TMM - -cholAux f st a = do - r <- createMatrix ColumnMajor n n - app2 f mat a mat r st - return r - where n = rows a - --- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. -cholH :: Matrix (Complex Double) -> Matrix (Complex Double) -cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat - --- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. -cholS :: Matrix Double -> Matrix Double -cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat - --- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version). -mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) -mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat - --- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version). -mbCholS :: Matrix Double -> Maybe (Matrix Double) -mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM -foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM - --- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. -qrR :: Matrix Double -> (Matrix Double, Vector Double) -qrR = qrAux dgeqr2 "qrR" . fmat - --- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. -qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) -qrC = qrAux zgeqr2 "qrC" . fmat - -qrAux f st a = unsafePerformIO $ do - r <- createMatrix ColumnMajor m n - tau <- createVector mn - app3 f mat a vec tau mat r st - return (r,tau) - where - m = rows a - n = cols a - mn = min m n - -foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM -foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM - --- | build rotation from reflectors -qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double -qrgrR = qrgrAux dorgqr "qrgrR" --- | build rotation from reflectors -qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) -qrgrC = qrgrAux zungqr "qrgrC" - -qrgrAux f st n (a, tau) = unsafePerformIO $ do - res <- createMatrix ColumnMajor (rows a) n - app3 f mat (fmat a) vec (subVector 0 n tau') mat res st - return res - where - tau' = vjoin [tau, constantD 0 n] - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM -foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM - --- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. -hessR :: Matrix Double -> (Matrix Double, Vector Double) -hessR = hessAux dgehrd "hessR" . fmat - --- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. -hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) -hessC = hessAux zgehrd "hessC" . fmat - -hessAux f st a = unsafePerformIO $ do - r <- createMatrix ColumnMajor m n - tau <- createVector (mn-1) - app3 f mat a vec tau mat r st - return (r,tau) - where m = rows a - n = cols a - mn = min m n - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "schur_l_R" dgees :: TMMM -foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM - --- | Schur factorization of a square real matrix, using LAPACK's /dgees/. -schurR :: Matrix Double -> (Matrix Double, Matrix Double) -schurR = schurAux dgees "schurR" . fmat - --- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. -schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) -schurC = schurAux zgees "schurC" . fmat - -schurAux f st a = unsafePerformIO $ do - u <- createMatrix ColumnMajor n n - s <- createMatrix ColumnMajor n n - app3 f mat a mat u mat s st - return (u,s) - where n = rows a - ------------------------------------------------------------------------------------ -foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM -foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM - --- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. -luR :: Matrix Double -> (Matrix Double, [Int]) -luR = luAux dgetrf "luR" . fmat - --- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. -luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) -luC = luAux zgetrf "luC" . fmat - -luAux f st a = unsafePerformIO $ do - lu <- createMatrix ColumnMajor n m - piv <- createVector (min n m) - app3 f mat a vec piv mat lu st - return (lu, map (pred.round) (toList piv)) - where n = rows a - m = cols a - ------------------------------------------------------------------------------------ -type TW a = CInt -> PD -> a -type TQ a = CInt -> CInt -> PC -> a - -foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM -foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) - --- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. -lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double -lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) - --- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. -lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) -lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) - -lusAux f st a piv b - | n1==n2 && n2==n =unsafePerformIO $ do - x <- createMatrix ColumnMajor n m - app4 f mat a vec piv' mat b mat x st - return x - | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" - where n1 = rows a - n2 = cols a - n = rows b - m = cols b - piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double - diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c deleted file mode 100644 index e5e45ef..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ /dev/null @@ -1,1489 +0,0 @@ -#include -#include -#include -#include -#include -#include "lapack-aux.h" - -#define MACRO(B) do {B} while (0) -#define ERROR(CODE) MACRO(return CODE;) -#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) - -#define MIN(A,B) ((A)<(B)?(A):(B)) -#define MAX(A,B) ((A)>(B)?(A):(B)) - -// #define DBGL - -#ifdef DBGL -#define DEBUGMSG(M) printf("\nLAPACK "M"\n"); -#else -#define DEBUGMSG(M) -#endif - -#define OK return 0; - -// #ifdef DBGL -// #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL); -// #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;); -// #else -// #define DEBUGMSG(M) -// #define OK return 0; -// #endif - -#define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \ - for(q=0;q=1 && ar==ac && ar==br,BAD_SIZE); - DEBUGMSG("linearSolveR_l"); - double*AC = (double*)malloc(n*n*sizeof(double)); - memcpy(AC,ap,n*n*sizeof(double)); - memcpy(xp,bp,n*nhrs*sizeof(double)); - integer * ipiv = (integer*)malloc(n*sizeof(integer)); - integer res; - dgesv_ (&n,&nhrs, - AC, &n, - ipiv, - xp, &n, - &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); - free(ipiv); - free(AC); - OK -} - -//////////////////// general complex linear system //////////// - -/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a, - integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer * - info); - -int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { - integer n = ar; - integer nhrs = bc; - REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); - DEBUGMSG("linearSolveC_l"); - doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); - memcpy(AC,ap,n*n*sizeof(doublecomplex)); - memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); - integer * ipiv = (integer*)malloc(n*sizeof(integer)); - integer res; - zgesv_ (&n,&nhrs, - AC, &n, - ipiv, - xp, &n, - &res); - if(res>0) { - return SINGULAR; - } - CHECK(res,res); - free(ipiv); - free(AC); - OK -} - -//////// symmetric positive definite real linear system using Cholesky //////////// - -/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs, - doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * - info); - -int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) { - integer n = ar; - integer nhrs = bc; - REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); - DEBUGMSG("cholSolveR_l"); - memcpy(xp,bp,n*nhrs*sizeof(double)); - integer res; - dpotrs_ ("U", - &n,&nhrs, - (double*)ap, &n, - xp, &n, - &res); - CHECK(res,res); - OK -} - -//////// Hermitian positive definite real linear system using Cholesky //////////// - -/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, - doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, - integer *info); - -int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { - integer n = ar; - integer nhrs = bc; - REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); - DEBUGMSG("cholSolveC_l"); - memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); - integer res; - zpotrs_ ("U", - &n,&nhrs, - (doublecomplex*)ap, &n, - xp, &n, - &res); - CHECK(res,res); - OK -} - -//////////////////// least squares real linear system //////////// - -/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer * - nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, - doublereal *work, integer *lwork, integer *info); - -int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) { - integer m = ar; - integer n = ac; - integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); - DEBUGMSG("linearSolveLSR_l"); - double*AC = (double*)malloc(m*n*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(double)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(double)); - } else { - int k; - for(k = 0; k0) { - return SINGULAR; - } - CHECK(res,res); - free(work); - free(AC); - OK -} - -//////////////////// least squares complex linear system //////////// - -/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer * - nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, - doublecomplex *work, integer *lwork, integer *info); - -int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) { - integer m = ar; - integer n = ac; - integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); - DEBUGMSG("linearSolveLSC_l"); - doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); - memcpy(AC,ap,m*n*sizeof(doublecomplex)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); - } else { - int k; - for(k = 0; k0) { - return SINGULAR; - } - CHECK(res,res); - free(work); - free(AC); - OK -} - -//////////////////// least squares real linear system using SVD //////////// - -/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, - doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * - s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, - integer *info); - -int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) { - integer m = ar; - integer n = ac; - integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); - DEBUGMSG("linearSolveSVDR_l"); - double*AC = (double*)malloc(m*n*sizeof(double)); - double*S = (double*)malloc(MIN(m,n)*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(double)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(double)); - } else { - int k; - for(k = 0; k0) { - return NOCONVER; - } - CHECK(res,res); - free(work); - free(S); - free(AC); - OK -} - -//////////////////// least squares complex linear system using SVD //////////// - -// not in clapack.h - -int zgelss_(integer *m, integer *n, integer *nhrs, - doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s, - doublereal *rcond, integer* rank, - doublecomplex *work, integer* lwork, doublereal* rwork, - integer *info); - -int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) { - integer m = ar; - integer n = ac; - integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); - DEBUGMSG("linearSolveSVDC_l"); - doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); - double*S = (double*)malloc(MIN(m,n)*sizeof(double)); - double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(doublecomplex)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); - } else { - int k; - for(k = 0; k0) { - return NOCONVER; - } - CHECK(res,res); - free(work); - free(RWORK); - free(S); - free(AC); - OK -} - -//////////////////// Cholesky factorization ///////////////////////// - -/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, - integer *lda, integer *info); - -int chol_l_H(KCMAT(a),CMAT(l)) { - integer n = ar; - REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE); - DEBUGMSG("chol_l_H"); - memcpy(lp,ap,n*n*sizeof(doublecomplex)); - char uplo = 'U'; - integer res; - zpotrf_ (&uplo,&n,lp,&n,&res); - CHECK(res>0,NODEFPOS); - CHECK(res,res); - doublecomplex zero = {0.,0.}; - int r,c; - for (r=0; r=1 && ac == n && lr==n && lc==n,BAD_SIZE); - DEBUGMSG("chol_l_S"); - memcpy(lp,ap,n*n*sizeof(double)); - char uplo = 'U'; - integer res; - dpotrf_ (&uplo,&n,lp,&n,&res); - CHECK(res>0,NODEFPOS); - CHECK(res,res); - int r,c; - for (r=0; r=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); - DEBUGMSG("qr_l_R"); - double *WORK = (double*)malloc(n*sizeof(double)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(double)); - integer res; - dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); - CHECK(res,res); - free(WORK); - OK -} - -/* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a, - integer *lda, doublecomplex *tau, doublecomplex *work, integer *info); - -int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { - integer m = ar; - integer n = ac; - integer mn = MIN(m,n); - REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); - DEBUGMSG("qr_l_C"); - doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); - integer res; - zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); - CHECK(res,res); - free(WORK); - OK -} - -/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal * - a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, - integer *info); - -int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) { - integer m = ar; - integer n = MIN(ac,ar); - integer k = taun; - DEBUGMSG("c_dorgqr"); - integer lwork = 8*n; // FIXME - double *WORK = (double*)malloc(lwork*sizeof(double)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*k*sizeof(double)); - integer res; - dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res); - CHECK(res,res); - free(WORK); - OK -} - -/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k, - doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * - work, integer *lwork, integer *info); - -int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) { - integer m = ar; - integer n = MIN(ac,ar); - integer k = taun; - DEBUGMSG("z_ungqr"); - integer lwork = 8*n; // FIXME - doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*k*sizeof(doublecomplex)); - integer res; - zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res); - CHECK(res,res); - free(WORK); - OK -} - - -//////////////////// Hessenberg factorization ///////////////////////// - -/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, - doublereal *a, integer *lda, doublereal *tau, doublereal *work, - integer *lwork, integer *info); - -int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) { - integer m = ar; - integer n = ac; - integer mn = MIN(m,n); - REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); - DEBUGMSG("hess_l_R"); - integer lwork = 5*n; // fixme - double *WORK = (double*)malloc(lwork*sizeof(double)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(double)); - integer res; - integer one = 1; - dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); - CHECK(res,res); - free(WORK); - OK -} - - -/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi, - doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * - work, integer *lwork, integer *info); - -int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { - integer m = ar; - integer n = ac; - integer mn = MIN(m,n); - REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); - DEBUGMSG("hess_l_C"); - integer lwork = 5*n; // fixme - doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); - integer res; - integer one = 1; - zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); - CHECK(res,res); - free(WORK); - OK -} - -//////////////////// Schur factorization ///////////////////////// - -/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n, - doublereal *a, integer *lda, integer *sdim, doublereal *wr, - doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, - integer *lwork, logical *bwork, integer *info); - -int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) { - integer m = ar; - integer n = ac; - REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); - DEBUGMSG("schur_l_R"); - //int k; - //printf("---------------------------\n"); - //printf("%p: ",ap); for(k=0;k0) { - return NOCONVER; - } - CHECK(res,res); - free(WR); - free(WI); - free(BWORK); - free(WORK); - OK -} - - -/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n, - doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w, - doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork, - doublereal *rwork, logical *bwork, integer *info); - -int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) { - integer m = ar; - integer n = ac; - REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); - DEBUGMSG("schur_l_C"); - memcpy(sp,ap,n*n*sizeof(doublecomplex)); - integer lwork = 6*n; // fixme - doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); - doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); - // W not really required in this call - logical *BWORK = (logical*)malloc(n*sizeof(logical)); - double *RWORK = (double*)malloc(n*sizeof(double)); - integer res; - integer sdim; - zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W, - up,&n, - WORK,&lwork,RWORK,BWORK,&res); - if(res>0) { - return NOCONVER; - } - CHECK(res,res); - free(W); - free(BWORK); - free(WORK); - OK -} - -//////////////////// LU factorization ///////////////////////// - -/* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer * - lda, integer *ipiv, integer *info); - -int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) { - integer m = ar; - integer n = ac; - integer mn = MIN(m,n); - REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); - DEBUGMSG("lu_l_R"); - integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); - memcpy(rp,ap,m*n*sizeof(double)); - integer res; - dgetrf_ (&m,&n,rp,&m,auxipiv,&res); - if(res>0) { - res = 0; // fixme - } - CHECK(res,res); - int k; - for (k=0; k=1 && n >=1 && ipivn == mn, BAD_SIZE); - DEBUGMSG("lu_l_C"); - integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); - integer res; - zgetrf_ (&m,&n,rp,&m,auxipiv,&res); - if(res>0) { - res = 0; // fixme - } - CHECK(res,res); - int k; - for (k=0; k0; - } - OK -} - -int stepD(DVEC(x),DVEC(y)) { - DEBUGMSG("stepD") - int k; - for(k=0;k0; - } - OK -} - -//////////////////// cond ///////////////////////// - -int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) { - REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); - DEBUGMSG("condF") - int k; - for(k=0;kyp[k]?gtp[k]:eqp[k]); - } - OK -} - -int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) { - REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); - DEBUGMSG("condD") - int k; - for(k=0;kyp[k]?gtp[k]:eqp[k]); - } - OK -} - diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h b/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h deleted file mode 100644 index a3f1899..0000000 --- a/packages/hmatrix/src/Numeric/LinearAlgebra/LAPACK/lapack-aux.h +++ /dev/null @@ -1,60 +0,0 @@ -/* - * We have copied the definitions in f2c.h required - * to compile clapack.h, modified to support both - * 32 and 64 bit - - http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib - http://www.ibm.com/developerworks/library/l-port64.html - */ - -#ifdef _LP64 -typedef int integer; -typedef unsigned int uinteger; -typedef int logical; -typedef long longint; /* system-dependent */ -typedef unsigned long ulongint; /* system-dependent */ -#else -typedef long int integer; -typedef unsigned long int uinteger; -typedef long int logical; -typedef long long longint; /* system-dependent */ -typedef unsigned long long ulongint; /* system-dependent */ -#endif - -typedef char *address; -typedef short int shortint; -typedef float real; -typedef double doublereal; -typedef struct { real r, i; } complex; -typedef struct { doublereal r, i; } doublecomplex; -typedef short int shortlogical; -typedef char logical1; -typedef char integer1; - -typedef logical (*L_fp)(); -typedef short ftnlen; - -/********************************************************/ - -#define FVEC(A) int A##n, float*A##p -#define DVEC(A) int A##n, double*A##p -#define QVEC(A) int A##n, complex*A##p -#define CVEC(A) int A##n, doublecomplex*A##p -#define PVEC(A) int A##n, void* A##p, int A##s -#define FMAT(A) int A##r, int A##c, float* A##p -#define DMAT(A) int A##r, int A##c, double* A##p -#define QMAT(A) int A##r, int A##c, complex* A##p -#define CMAT(A) int A##r, int A##c, doublecomplex* A##p -#define PMAT(A) int A##r, int A##c, void* A##p, int A##s - -#define KFVEC(A) int A##n, const float*A##p -#define KDVEC(A) int A##n, const double*A##p -#define KQVEC(A) int A##n, const complex*A##p -#define KCVEC(A) int A##n, const doublecomplex*A##p -#define KPVEC(A) int A##n, const void* A##p, int A##s -#define KFMAT(A) int A##r, int A##c, const float* A##p -#define KDMAT(A) int A##r, int A##c, const double* A##p -#define KQMAT(A) int A##r, int A##c, const complex* A##p -#define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p -#define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s - -- cgit v1.2.3