From 551cf7498c33bc0948bb4cb8444ae6f8af7278ea Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 13:43:07 +0200 Subject: separation ok --- packages/base/hmatrix-base.cabal | 6 +- packages/base/src/Data/Packed/Development.hs | 3 +- packages/base/src/Data/Packed/Vector.hs | 3 +- packages/base/src/Numeric/LAPACK.hs | 552 ---------------------- packages/base/src/Numeric/LinearAlgebra/LAPACK.hs | 552 ++++++++++++++++++++++ 5 files changed, 559 insertions(+), 557 deletions(-) delete mode 100644 packages/base/src/Numeric/LAPACK.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra/LAPACK.hs (limited to 'packages/base') diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal index 0413d4a..bd356df 100644 --- a/packages/base/hmatrix-base.cabal +++ b/packages/base/hmatrix-base.cabal @@ -7,7 +7,8 @@ Maintainer: Alberto Ruiz Stability: provisional Homepage: https://github.com/albertoruiz/hmatrix Synopsis: Numeric Linear Algebra -Description: Numeric Linear Algebra +Description: Basic support for dense linear algebra + internally implemented using BLAS and LAPACK. Category: Math tested-with: GHC==7.8 @@ -35,7 +36,7 @@ library Data.Packed.ST, Data.Packed.Development, Numeric.Conversion - Numeric.LAPACK + Numeric.LinearAlgebra.LAPACK other-modules: Data.Packed.Internal, Data.Packed.Internal.Common, @@ -52,7 +53,6 @@ library ghc-options: -Wall -fno-warn-missing-signatures -fno-warn-orphans --- -fno-warn-unused-binds cc-options: -O4 -msse2 -Wall diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs index 9350acb..6dfd551 100644 --- a/packages/base/src/Data/Packed/Development.hs +++ b/packages/base/src/Data/Packed/Development.hs @@ -24,7 +24,8 @@ module Data.Packed.Development ( unsafeFromForeignPtr, unsafeToForeignPtr, check, (//), - at', atM', fi, table + at', atM', fi, table, + conformMs, conformVs, shSize, splitEvery ) where import Data.Packed.Internal diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs index 653a257..a8f6bb1 100644 --- a/packages/base/src/Data/Packed/Vector.hs +++ b/packages/base/src/Data/Packed/Vector.hs @@ -23,7 +23,8 @@ module Data.Packed.Vector ( subVector, takesV, vjoin, join, mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, - foldLoop, foldVector, foldVectorG, foldVectorWithIndex + foldLoop, foldVector, foldVectorG, foldVectorWithIndex, + stepD, stepF, condD, condF, conjugateC, conjugateQ ) where import Data.Packed.Internal.Vector diff --git a/packages/base/src/Numeric/LAPACK.hs b/packages/base/src/Numeric/LAPACK.hs deleted file mode 100644 index 08cd759..0000000 --- a/packages/base/src/Numeric/LAPACK.hs +++ /dev/null @@ -1,552 +0,0 @@ ------------------------------------------------------------------------------ --- | --- 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/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs new file mode 100644 index 0000000..af0c744 --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs @@ -0,0 +1,552 @@ +----------------------------------------------------------------------------- +-- | +-- Module : Numeric.LinearAlgebra.LAPACK +-- Copyright : (c) Alberto Ruiz 2006-14 +-- License : BSD3 +-- +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +-- Functional interface to selected LAPACK functions (). +-- +----------------------------------------------------------------------------- + +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.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 + -- cgit v1.2.3