From 5d967adae38d7fe80443b57b40ae89cd15db5e18 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 5 Jun 2015 16:45:23 +0200 Subject: move lapack --- packages/base/src/Numeric/LinearAlgebra/LAPACK.hs | 569 ---------------------- 1 file changed, 569 deletions(-) delete mode 100644 packages/base/src/Numeric/LinearAlgebra/LAPACK.hs (limited to 'packages/base/src/Numeric') diff --git a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs b/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs deleted file mode 100644 index 6fb2b13..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/LAPACK.hs +++ /dev/null @@ -1,569 +0,0 @@ ------------------------------------------------------------------------------ --- | --- Module : Numeric.LinearAlgebra.LAPACK --- Copyright : (c) Alberto Ruiz 2006-14 --- License : BSD3 --- Maintainer : Alberto Ruiz --- Stability : provisional --- --- Functional interface to selected LAPACK functions (). --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - - -module Numeric.LinearAlgebra.LAPACK ( - -- * Matrix product - multiplyR, multiplyC, multiplyF, multiplyQ, multiplyI, - -- * Linear systems - linearSolveR, linearSolveC, - mbLinearSolveR, mbLinearSolveC, - 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 -foreign import ccall unsafe "multiplyI" c_multiplyI :: OM CInt (OM CInt (OM CInt (IO CInt))) - -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 - -multiplyI :: Matrix CInt -> Matrix CInt -> Matrix CInt -multiplyI a b = unsafePerformIO $ do - when (cols a /= rows b) $ error $ - "inconsistent dimensions in matrix product "++ shSize a ++ " x " ++ shSize b - s <- createMatrix ColumnMajor (rows a) (cols b) - app3 c_multiplyI omat a omat b omat s "c_multiplyI" - return s - ------------------------------------------------------------------------------ -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,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,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,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 0) "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 g f st a b - | n1==n2 && n1==r = unsafePerformIO . g $ 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 id dgesv "linearSolveR" (fmat a) (fmat b) - -mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double) -mbLinearSolveR a b = linearSolveSQAux mbCatch 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 id zgesv "linearSolveC" (fmat a) (fmat b) - -mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) -mbLinearSolveC a b = linearSolveSQAux mbCatch 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 id 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 id 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