diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:18:56 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-05-08 12:18:56 +0200 |
commit | 561a6c0e21bb77c21114ccbbd86d3af5ddb5a3f1 (patch) | |
tree | d49b67d75d63938229f2d5cbed5c49d06dc02bcf /packages/hmatrix | |
parent | 5992d92357cfd911c8f2e9f5faaa4fd8a323fd9a (diff) |
Conversion, LAPACK -> base
Diffstat (limited to 'packages/hmatrix')
4 files changed, 0 insertions, 2195 deletions
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 @@ | |||
1 | {-# LANGUAGE TypeFamilies #-} | ||
2 | {-# LANGUAGE FlexibleContexts #-} | ||
3 | {-# LANGUAGE FlexibleInstances #-} | ||
4 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
5 | {-# LANGUAGE FunctionalDependencies #-} | ||
6 | {-# LANGUAGE UndecidableInstances #-} | ||
7 | |||
8 | ----------------------------------------------------------------------------- | ||
9 | -- | | ||
10 | -- Module : Numeric.Conversion | ||
11 | -- Copyright : (c) Alberto Ruiz 2010 | ||
12 | -- License : GPL-style | ||
13 | -- | ||
14 | -- Maintainer : Alberto Ruiz <aruiz@um.es> | ||
15 | -- Stability : provisional | ||
16 | -- Portability : portable | ||
17 | -- | ||
18 | -- Conversion routines | ||
19 | -- | ||
20 | ----------------------------------------------------------------------------- | ||
21 | |||
22 | module Numeric.Conversion ( | ||
23 | Complexable(..), RealElement, | ||
24 | module Data.Complex | ||
25 | ) where | ||
26 | |||
27 | import Data.Packed.Internal.Vector | ||
28 | import Data.Packed.Internal.Matrix | ||
29 | import Data.Complex | ||
30 | import Control.Arrow((***)) | ||
31 | |||
32 | ------------------------------------------------------------------- | ||
33 | |||
34 | -- | Supported single-double precision type pairs | ||
35 | class (Element s, Element d) => Precision s d | s -> d, d -> s where | ||
36 | double2FloatG :: Vector d -> Vector s | ||
37 | float2DoubleG :: Vector s -> Vector d | ||
38 | |||
39 | instance Precision Float Double where | ||
40 | double2FloatG = double2FloatV | ||
41 | float2DoubleG = float2DoubleV | ||
42 | |||
43 | instance Precision (Complex Float) (Complex Double) where | ||
44 | double2FloatG = asComplex . double2FloatV . asReal | ||
45 | float2DoubleG = asComplex . float2DoubleV . asReal | ||
46 | |||
47 | -- | Supported real types | ||
48 | class (Element t, Element (Complex t), RealFloat t | ||
49 | -- , RealOf t ~ t, RealOf (Complex t) ~ t | ||
50 | ) | ||
51 | => RealElement t | ||
52 | |||
53 | instance RealElement Double | ||
54 | instance RealElement Float | ||
55 | |||
56 | |||
57 | -- | Structures that may contain complex numbers | ||
58 | class Complexable c where | ||
59 | toComplex' :: (RealElement e) => (c e, c e) -> c (Complex e) | ||
60 | fromComplex' :: (RealElement e) => c (Complex e) -> (c e, c e) | ||
61 | comp' :: (RealElement e) => c e -> c (Complex e) | ||
62 | single' :: Precision a b => c b -> c a | ||
63 | double' :: Precision a b => c a -> c b | ||
64 | |||
65 | |||
66 | instance Complexable Vector where | ||
67 | toComplex' = toComplexV | ||
68 | fromComplex' = fromComplexV | ||
69 | comp' v = toComplex' (v,constantD 0 (dim v)) | ||
70 | single' = double2FloatG | ||
71 | double' = float2DoubleG | ||
72 | |||
73 | |||
74 | -- | creates a complex vector from vectors with real and imaginary parts | ||
75 | toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) | ||
76 | toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] | ||
77 | |||
78 | -- | the inverse of 'toComplex' | ||
79 | fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) | ||
80 | fromComplexV z = (r,i) where | ||
81 | [r,i] = toColumns $ reshape 2 $ asReal z | ||
82 | |||
83 | |||
84 | instance Complexable Matrix where | ||
85 | toComplex' = uncurry $ liftMatrix2 $ curry toComplex' | ||
86 | fromComplex' z = (reshape c *** reshape c) . fromComplex' . flatten $ z | ||
87 | where c = cols z | ||
88 | comp' = liftMatrix comp' | ||
89 | single' = liftMatrix single' | ||
90 | double' = liftMatrix double' | ||
91 | |||
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 @@ | |||
1 | ----------------------------------------------------------------------------- | ||
2 | -- | | ||
3 | -- Module : Numeric.LinearAlgebra.LAPACK | ||
4 | -- Copyright : (c) Alberto Ruiz 2006-7 | ||
5 | -- License : GPL-style | ||
6 | -- | ||
7 | -- Maintainer : Alberto Ruiz (aruiz at um dot es) | ||
8 | -- Stability : provisional | ||
9 | -- Portability : portable (uses FFI) | ||
10 | -- | ||
11 | -- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>). | ||
12 | -- | ||
13 | ----------------------------------------------------------------------------- | ||
14 | {-# OPTIONS_HADDOCK hide #-} | ||
15 | |||
16 | module Numeric.LinearAlgebra.LAPACK ( | ||
17 | -- * Matrix product | ||
18 | multiplyR, multiplyC, multiplyF, multiplyQ, | ||
19 | -- * Linear systems | ||
20 | linearSolveR, linearSolveC, | ||
21 | lusR, lusC, | ||
22 | cholSolveR, cholSolveC, | ||
23 | linearSolveLSR, linearSolveLSC, | ||
24 | linearSolveSVDR, linearSolveSVDC, | ||
25 | -- * SVD | ||
26 | svR, svRd, svC, svCd, | ||
27 | svdR, svdRd, svdC, svdCd, | ||
28 | thinSVDR, thinSVDRd, thinSVDC, thinSVDCd, | ||
29 | rightSVR, rightSVC, leftSVR, leftSVC, | ||
30 | -- * Eigensystems | ||
31 | eigR, eigC, eigS, eigS', eigH, eigH', | ||
32 | eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, | ||
33 | -- * LU | ||
34 | luR, luC, | ||
35 | -- * Cholesky | ||
36 | cholS, cholH, mbCholS, mbCholH, | ||
37 | -- * QR | ||
38 | qrR, qrC, qrgrR, qrgrC, | ||
39 | -- * Hessenberg | ||
40 | hessR, hessC, | ||
41 | -- * Schur | ||
42 | schurR, schurC | ||
43 | ) where | ||
44 | |||
45 | import Data.Packed.Internal | ||
46 | import Data.Packed.Matrix | ||
47 | import Numeric.Conversion | ||
48 | import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) | ||
49 | |||
50 | import Foreign.Ptr(nullPtr) | ||
51 | import Foreign.C.Types | ||
52 | import Control.Monad(when) | ||
53 | import System.IO.Unsafe(unsafePerformIO) | ||
54 | |||
55 | ----------------------------------------------------------------------------------- | ||
56 | |||
57 | foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM | ||
58 | foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TCMCMCM | ||
59 | foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TFMFMFM | ||
60 | foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TQMQMQM | ||
61 | |||
62 | isT Matrix{order = ColumnMajor} = 0 | ||
63 | isT Matrix{order = RowMajor} = 1 | ||
64 | |||
65 | tt x@Matrix{order = ColumnMajor} = x | ||
66 | tt x@Matrix{order = RowMajor} = trans x | ||
67 | |||
68 | multiplyAux f st a b = unsafePerformIO $ do | ||
69 | when (cols a /= rows b) $ error $ "inconsistent dimensions in matrix product "++ | ||
70 | show (rows a,cols a) ++ " x " ++ show (rows b, cols b) | ||
71 | s <- createMatrix ColumnMajor (rows a) (cols b) | ||
72 | app3 (f (isT a) (isT b)) mat (tt a) mat (tt b) mat s st | ||
73 | return s | ||
74 | |||
75 | -- | Matrix product based on BLAS's /dgemm/. | ||
76 | multiplyR :: Matrix Double -> Matrix Double -> Matrix Double | ||
77 | multiplyR a b = {-# SCC "multiplyR" #-} multiplyAux dgemmc "dgemmc" a b | ||
78 | |||
79 | -- | Matrix product based on BLAS's /zgemm/. | ||
80 | multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
81 | multiplyC a b = multiplyAux zgemmc "zgemmc" a b | ||
82 | |||
83 | -- | Matrix product based on BLAS's /sgemm/. | ||
84 | multiplyF :: Matrix Float -> Matrix Float -> Matrix Float | ||
85 | multiplyF a b = multiplyAux sgemmc "sgemmc" a b | ||
86 | |||
87 | -- | Matrix product based on BLAS's /cgemm/. | ||
88 | multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float) | ||
89 | multiplyQ a b = multiplyAux cgemmc "cgemmc" a b | ||
90 | |||
91 | ----------------------------------------------------------------------------- | ||
92 | foreign import ccall unsafe "svd_l_R" dgesvd :: TMMVM | ||
93 | foreign import ccall unsafe "svd_l_C" zgesvd :: TCMCMVCM | ||
94 | foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TMMVM | ||
95 | foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TCMCMVCM | ||
96 | |||
97 | -- | Full SVD of a real matrix using LAPACK's /dgesvd/. | ||
98 | svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
99 | svdR = svdAux dgesvd "svdR" . fmat | ||
100 | |||
101 | -- | Full SVD of a real matrix using LAPACK's /dgesdd/. | ||
102 | svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
103 | svdRd = svdAux dgesdd "svdRdd" . fmat | ||
104 | |||
105 | -- | Full SVD of a complex matrix using LAPACK's /zgesvd/. | ||
106 | svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
107 | svdC = svdAux zgesvd "svdC" . fmat | ||
108 | |||
109 | -- | Full SVD of a complex matrix using LAPACK's /zgesdd/. | ||
110 | svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
111 | svdCd = svdAux zgesdd "svdCdd" . fmat | ||
112 | |||
113 | svdAux f st x = unsafePerformIO $ do | ||
114 | u <- createMatrix ColumnMajor r r | ||
115 | s <- createVector (min r c) | ||
116 | v <- createMatrix ColumnMajor c c | ||
117 | app4 f mat x mat u vec s mat v st | ||
118 | return (u,s,trans v) | ||
119 | where r = rows x | ||
120 | c = cols x | ||
121 | |||
122 | |||
123 | -- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. | ||
124 | thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
125 | thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat | ||
126 | |||
127 | -- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. | ||
128 | thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
129 | thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat | ||
130 | |||
131 | -- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. | ||
132 | thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) | ||
133 | thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat | ||
134 | |||
135 | -- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. | ||
136 | thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) | ||
137 | thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat | ||
138 | |||
139 | thinSVDAux f st x = unsafePerformIO $ do | ||
140 | u <- createMatrix ColumnMajor r q | ||
141 | s <- createVector q | ||
142 | v <- createMatrix ColumnMajor q c | ||
143 | app4 f mat x mat u vec s mat v st | ||
144 | return (u,s,trans v) | ||
145 | where r = rows x | ||
146 | c = cols x | ||
147 | q = min r c | ||
148 | |||
149 | |||
150 | -- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. | ||
151 | svR :: Matrix Double -> Vector Double | ||
152 | svR = svAux dgesvd "svR" . fmat | ||
153 | |||
154 | -- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. | ||
155 | svC :: Matrix (Complex Double) -> Vector Double | ||
156 | svC = svAux zgesvd "svC" . fmat | ||
157 | |||
158 | -- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. | ||
159 | svRd :: Matrix Double -> Vector Double | ||
160 | svRd = svAux dgesdd "svRd" . fmat | ||
161 | |||
162 | -- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. | ||
163 | svCd :: Matrix (Complex Double) -> Vector Double | ||
164 | svCd = svAux zgesdd "svCd" . fmat | ||
165 | |||
166 | svAux f st x = unsafePerformIO $ do | ||
167 | s <- createVector q | ||
168 | app2 g mat x vec s st | ||
169 | return s | ||
170 | where r = rows x | ||
171 | c = cols x | ||
172 | q = min r c | ||
173 | g ra ca pa nb pb = f ra ca pa 0 0 nullPtr nb pb 0 0 nullPtr | ||
174 | |||
175 | |||
176 | -- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. | ||
177 | rightSVR :: Matrix Double -> (Vector Double, Matrix Double) | ||
178 | rightSVR = rightSVAux dgesvd "rightSVR" . fmat | ||
179 | |||
180 | -- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. | ||
181 | rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
182 | rightSVC = rightSVAux zgesvd "rightSVC" . fmat | ||
183 | |||
184 | rightSVAux f st x = unsafePerformIO $ do | ||
185 | s <- createVector q | ||
186 | v <- createMatrix ColumnMajor c c | ||
187 | app3 g mat x vec s mat v st | ||
188 | return (s,trans v) | ||
189 | where r = rows x | ||
190 | c = cols x | ||
191 | q = min r c | ||
192 | g ra ca pa = f ra ca pa 0 0 nullPtr | ||
193 | |||
194 | |||
195 | -- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. | ||
196 | leftSVR :: Matrix Double -> (Matrix Double, Vector Double) | ||
197 | leftSVR = leftSVAux dgesvd "leftSVR" . fmat | ||
198 | |||
199 | -- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. | ||
200 | leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) | ||
201 | leftSVC = leftSVAux zgesvd "leftSVC" . fmat | ||
202 | |||
203 | leftSVAux f st x = unsafePerformIO $ do | ||
204 | u <- createMatrix ColumnMajor r r | ||
205 | s <- createVector q | ||
206 | app3 g mat x mat u vec s st | ||
207 | return (u,s) | ||
208 | where r = rows x | ||
209 | c = cols x | ||
210 | q = min r c | ||
211 | g ra ca pa ru cu pu nb pb = f ra ca pa ru cu pu nb pb 0 0 nullPtr | ||
212 | |||
213 | ----------------------------------------------------------------------------- | ||
214 | |||
215 | foreign import ccall unsafe "eig_l_R" dgeev :: TMMCVM | ||
216 | foreign import ccall unsafe "eig_l_C" zgeev :: TCMCMCVCM | ||
217 | foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> TMVM | ||
218 | foreign import ccall unsafe "eig_l_H" zheev :: CInt -> TCMVCM | ||
219 | |||
220 | eigAux f st m = unsafePerformIO $ do | ||
221 | l <- createVector r | ||
222 | v <- createMatrix ColumnMajor r r | ||
223 | app3 g mat m vec l mat v st | ||
224 | return (l,v) | ||
225 | where r = rows m | ||
226 | g ra ca pa = f ra ca pa 0 0 nullPtr | ||
227 | |||
228 | |||
229 | -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. | ||
230 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. | ||
231 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) | ||
232 | eigC = eigAux zgeev "eigC" . fmat | ||
233 | |||
234 | eigOnlyAux f st m = unsafePerformIO $ do | ||
235 | l <- createVector r | ||
236 | app2 g mat m vec l st | ||
237 | return l | ||
238 | where r = rows m | ||
239 | g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr | ||
240 | |||
241 | -- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. | ||
242 | -- The eigenvalues are not sorted. | ||
243 | eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) | ||
244 | eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat | ||
245 | |||
246 | -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. | ||
247 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. | ||
248 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) | ||
249 | eigR m = (s', v'') | ||
250 | where (s,v) = eigRaux (fmat m) | ||
251 | s' = fixeig1 s | ||
252 | v' = toRows $ trans v | ||
253 | v'' = fromColumns $ fixeig (toList s') v' | ||
254 | |||
255 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) | ||
256 | eigRaux m = unsafePerformIO $ do | ||
257 | l <- createVector r | ||
258 | v <- createMatrix ColumnMajor r r | ||
259 | app3 g mat m vec l mat v "eigR" | ||
260 | return (l,v) | ||
261 | where r = rows m | ||
262 | g ra ca pa = dgeev ra ca pa 0 0 nullPtr | ||
263 | |||
264 | fixeig1 s = toComplex' (subVector 0 r (asReal s), subVector r r (asReal s)) | ||
265 | where r = dim s | ||
266 | |||
267 | fixeig [] _ = [] | ||
268 | fixeig [_] [v] = [comp' v] | ||
269 | fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | ||
270 | | r1 == r2 && i1 == (-i2) = toComplex' (v1,v2) : toComplex' (v1,scale (-1) v2) : fixeig r vs | ||
271 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) | ||
272 | where scale = vectorMapValR Scale | ||
273 | fixeig _ _ = error "fixeig with impossible inputs" | ||
274 | |||
275 | |||
276 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | ||
277 | -- The eigenvalues are not sorted. | ||
278 | eigOnlyR :: Matrix Double -> Vector (Complex Double) | ||
279 | eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat | ||
280 | |||
281 | |||
282 | ----------------------------------------------------------------------------- | ||
283 | |||
284 | eigSHAux f st m = unsafePerformIO $ do | ||
285 | l <- createVector r | ||
286 | v <- createMatrix ColumnMajor r r | ||
287 | app3 f mat m vec l mat v st | ||
288 | return (l,v) | ||
289 | where r = rows m | ||
290 | |||
291 | -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. | ||
292 | -- The eigenvectors are the columns of v. | ||
293 | -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). | ||
294 | eigS :: Matrix Double -> (Vector Double, Matrix Double) | ||
295 | eigS m = (s', fliprl v) | ||
296 | where (s,v) = eigS' (fmat m) | ||
297 | s' = fromList . reverse . toList $ s | ||
298 | |||
299 | -- | 'eigS' in ascending order | ||
300 | eigS' :: Matrix Double -> (Vector Double, Matrix Double) | ||
301 | eigS' = eigSHAux (dsyev 1) "eigS'" . fmat | ||
302 | |||
303 | -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. | ||
304 | -- The eigenvectors are the columns of v. | ||
305 | -- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order). | ||
306 | eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
307 | eigH m = (s', fliprl v) | ||
308 | where (s,v) = eigH' (fmat m) | ||
309 | s' = fromList . reverse . toList $ s | ||
310 | |||
311 | -- | 'eigH' in ascending order | ||
312 | eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
313 | eigH' = eigSHAux (zheev 1) "eigH'" . fmat | ||
314 | |||
315 | |||
316 | -- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. | ||
317 | -- The eigenvalues are sorted in descending order. | ||
318 | eigOnlyS :: Matrix Double -> Vector Double | ||
319 | eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat | ||
320 | |||
321 | -- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. | ||
322 | -- The eigenvalues are sorted in descending order. | ||
323 | eigOnlyH :: Matrix (Complex Double) -> Vector Double | ||
324 | eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat | ||
325 | |||
326 | vrev = flatten . flipud . reshape 1 | ||
327 | |||
328 | ----------------------------------------------------------------------------- | ||
329 | foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM | ||
330 | foreign import ccall unsafe "linearSolveC_l" zgesv :: TCMCMCM | ||
331 | foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM | ||
332 | foreign import ccall unsafe "cholSolveC_l" zpotrs :: TCMCMCM | ||
333 | |||
334 | linearSolveSQAux f st a b | ||
335 | | n1==n2 && n1==r = unsafePerformIO $ do | ||
336 | s <- createMatrix ColumnMajor r c | ||
337 | app3 f mat a mat b mat s st | ||
338 | return s | ||
339 | | otherwise = error $ st ++ " of nonsquare matrix" | ||
340 | where n1 = rows a | ||
341 | n2 = cols a | ||
342 | r = rows b | ||
343 | c = cols b | ||
344 | |||
345 | -- | 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'. | ||
346 | linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
347 | linearSolveR a b = linearSolveSQAux dgesv "linearSolveR" (fmat a) (fmat b) | ||
348 | |||
349 | -- | 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'. | ||
350 | linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
351 | linearSolveC a b = linearSolveSQAux zgesv "linearSolveC" (fmat a) (fmat b) | ||
352 | |||
353 | |||
354 | -- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. | ||
355 | cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double | ||
356 | cholSolveR a b = linearSolveSQAux dpotrs "cholSolveR" (fmat a) (fmat b) | ||
357 | |||
358 | -- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. | ||
359 | cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
360 | cholSolveC a b = linearSolveSQAux zpotrs "cholSolveC" (fmat a) (fmat b) | ||
361 | |||
362 | ----------------------------------------------------------------------------------- | ||
363 | foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM | ||
364 | foreign import ccall unsafe "linearSolveLSC_l" zgels :: TCMCMCM | ||
365 | foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM | ||
366 | foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TCMCMCM | ||
367 | |||
368 | linearSolveAux f st a b = unsafePerformIO $ do | ||
369 | r <- createMatrix ColumnMajor (max m n) nrhs | ||
370 | app3 f mat a mat b mat r st | ||
371 | return r | ||
372 | where m = rows a | ||
373 | n = cols a | ||
374 | nrhs = cols b | ||
375 | |||
376 | -- | 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'. | ||
377 | linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double | ||
378 | linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ | ||
379 | linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) | ||
380 | |||
381 | -- | 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'. | ||
382 | linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
383 | linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ | ||
384 | linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) | ||
385 | |||
386 | -- | 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. | ||
387 | linearSolveSVDR :: Maybe Double -- ^ rcond | ||
388 | -> Matrix Double -- ^ coefficient matrix | ||
389 | -> Matrix Double -- ^ right hand sides (as columns) | ||
390 | -> Matrix Double -- ^ solution vectors (as columns) | ||
391 | linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ | ||
392 | linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) | ||
393 | linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) | ||
394 | |||
395 | -- | 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. | ||
396 | linearSolveSVDC :: Maybe Double -- ^ rcond | ||
397 | -> Matrix (Complex Double) -- ^ coefficient matrix | ||
398 | -> Matrix (Complex Double) -- ^ right hand sides (as columns) | ||
399 | -> Matrix (Complex Double) -- ^ solution vectors (as columns) | ||
400 | linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ | ||
401 | linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) | ||
402 | linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) | ||
403 | |||
404 | ----------------------------------------------------------------------------------- | ||
405 | foreign import ccall unsafe "chol_l_H" zpotrf :: TCMCM | ||
406 | foreign import ccall unsafe "chol_l_S" dpotrf :: TMM | ||
407 | |||
408 | cholAux f st a = do | ||
409 | r <- createMatrix ColumnMajor n n | ||
410 | app2 f mat a mat r st | ||
411 | return r | ||
412 | where n = rows a | ||
413 | |||
414 | -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/. | ||
415 | cholH :: Matrix (Complex Double) -> Matrix (Complex Double) | ||
416 | cholH = unsafePerformIO . cholAux zpotrf "cholH" . fmat | ||
417 | |||
418 | -- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/. | ||
419 | cholS :: Matrix Double -> Matrix Double | ||
420 | cholS = unsafePerformIO . cholAux dpotrf "cholS" . fmat | ||
421 | |||
422 | -- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version). | ||
423 | mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) | ||
424 | mbCholH = unsafePerformIO . mbCatch . cholAux zpotrf "cholH" . fmat | ||
425 | |||
426 | -- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/ ('Maybe' version). | ||
427 | mbCholS :: Matrix Double -> Maybe (Matrix Double) | ||
428 | mbCholS = unsafePerformIO . mbCatch . cholAux dpotrf "cholS" . fmat | ||
429 | |||
430 | ----------------------------------------------------------------------------------- | ||
431 | foreign import ccall unsafe "qr_l_R" dgeqr2 :: TMVM | ||
432 | foreign import ccall unsafe "qr_l_C" zgeqr2 :: TCMCVCM | ||
433 | |||
434 | -- | QR factorization of a real matrix, using LAPACK's /dgeqr2/. | ||
435 | qrR :: Matrix Double -> (Matrix Double, Vector Double) | ||
436 | qrR = qrAux dgeqr2 "qrR" . fmat | ||
437 | |||
438 | -- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/. | ||
439 | qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) | ||
440 | qrC = qrAux zgeqr2 "qrC" . fmat | ||
441 | |||
442 | qrAux f st a = unsafePerformIO $ do | ||
443 | r <- createMatrix ColumnMajor m n | ||
444 | tau <- createVector mn | ||
445 | app3 f mat a vec tau mat r st | ||
446 | return (r,tau) | ||
447 | where | ||
448 | m = rows a | ||
449 | n = cols a | ||
450 | mn = min m n | ||
451 | |||
452 | foreign import ccall unsafe "c_dorgqr" dorgqr :: TMVM | ||
453 | foreign import ccall unsafe "c_zungqr" zungqr :: TCMCVCM | ||
454 | |||
455 | -- | build rotation from reflectors | ||
456 | qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double | ||
457 | qrgrR = qrgrAux dorgqr "qrgrR" | ||
458 | -- | build rotation from reflectors | ||
459 | qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double) | ||
460 | qrgrC = qrgrAux zungqr "qrgrC" | ||
461 | |||
462 | qrgrAux f st n (a, tau) = unsafePerformIO $ do | ||
463 | res <- createMatrix ColumnMajor (rows a) n | ||
464 | app3 f mat (fmat a) vec (subVector 0 n tau') mat res st | ||
465 | return res | ||
466 | where | ||
467 | tau' = vjoin [tau, constantD 0 n] | ||
468 | |||
469 | ----------------------------------------------------------------------------------- | ||
470 | foreign import ccall unsafe "hess_l_R" dgehrd :: TMVM | ||
471 | foreign import ccall unsafe "hess_l_C" zgehrd :: TCMCVCM | ||
472 | |||
473 | -- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/. | ||
474 | hessR :: Matrix Double -> (Matrix Double, Vector Double) | ||
475 | hessR = hessAux dgehrd "hessR" . fmat | ||
476 | |||
477 | -- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/. | ||
478 | hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double)) | ||
479 | hessC = hessAux zgehrd "hessC" . fmat | ||
480 | |||
481 | hessAux f st a = unsafePerformIO $ do | ||
482 | r <- createMatrix ColumnMajor m n | ||
483 | tau <- createVector (mn-1) | ||
484 | app3 f mat a vec tau mat r st | ||
485 | return (r,tau) | ||
486 | where m = rows a | ||
487 | n = cols a | ||
488 | mn = min m n | ||
489 | |||
490 | ----------------------------------------------------------------------------------- | ||
491 | foreign import ccall unsafe "schur_l_R" dgees :: TMMM | ||
492 | foreign import ccall unsafe "schur_l_C" zgees :: TCMCMCM | ||
493 | |||
494 | -- | Schur factorization of a square real matrix, using LAPACK's /dgees/. | ||
495 | schurR :: Matrix Double -> (Matrix Double, Matrix Double) | ||
496 | schurR = schurAux dgees "schurR" . fmat | ||
497 | |||
498 | -- | Schur factorization of a square complex matrix, using LAPACK's /zgees/. | ||
499 | schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double)) | ||
500 | schurC = schurAux zgees "schurC" . fmat | ||
501 | |||
502 | schurAux f st a = unsafePerformIO $ do | ||
503 | u <- createMatrix ColumnMajor n n | ||
504 | s <- createMatrix ColumnMajor n n | ||
505 | app3 f mat a mat u mat s st | ||
506 | return (u,s) | ||
507 | where n = rows a | ||
508 | |||
509 | ----------------------------------------------------------------------------------- | ||
510 | foreign import ccall unsafe "lu_l_R" dgetrf :: TMVM | ||
511 | foreign import ccall unsafe "lu_l_C" zgetrf :: TCMVCM | ||
512 | |||
513 | -- | LU factorization of a general real matrix, using LAPACK's /dgetrf/. | ||
514 | luR :: Matrix Double -> (Matrix Double, [Int]) | ||
515 | luR = luAux dgetrf "luR" . fmat | ||
516 | |||
517 | -- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/. | ||
518 | luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int]) | ||
519 | luC = luAux zgetrf "luC" . fmat | ||
520 | |||
521 | luAux f st a = unsafePerformIO $ do | ||
522 | lu <- createMatrix ColumnMajor n m | ||
523 | piv <- createVector (min n m) | ||
524 | app3 f mat a vec piv mat lu st | ||
525 | return (lu, map (pred.round) (toList piv)) | ||
526 | where n = rows a | ||
527 | m = cols a | ||
528 | |||
529 | ----------------------------------------------------------------------------------- | ||
530 | type TW a = CInt -> PD -> a | ||
531 | type TQ a = CInt -> CInt -> PC -> a | ||
532 | |||
533 | foreign import ccall unsafe "luS_l_R" dgetrs :: TMVMM | ||
534 | foreign import ccall unsafe "luS_l_C" zgetrs :: TQ (TW (TQ (TQ (IO CInt)))) | ||
535 | |||
536 | -- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/. | ||
537 | lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double | ||
538 | lusR a piv b = lusAux dgetrs "lusR" (fmat a) piv (fmat b) | ||
539 | |||
540 | -- | Solve a real linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/. | ||
541 | lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double) | ||
542 | lusC a piv b = lusAux zgetrs "lusC" (fmat a) piv (fmat b) | ||
543 | |||
544 | lusAux f st a piv b | ||
545 | | n1==n2 && n2==n =unsafePerformIO $ do | ||
546 | x <- createMatrix ColumnMajor n m | ||
547 | app4 f mat a vec piv' mat b mat x st | ||
548 | return x | ||
549 | | otherwise = error $ st ++ " on LU factorization of nonsquare matrix" | ||
550 | where n1 = rows a | ||
551 | n2 = cols a | ||
552 | n = rows b | ||
553 | m = cols b | ||
554 | piv' = fromList (map (fromIntegral.succ) piv) :: Vector Double | ||
555 | |||
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 @@ | |||
1 | #include <stdio.h> | ||
2 | #include <stdlib.h> | ||
3 | #include <string.h> | ||
4 | #include <math.h> | ||
5 | #include <time.h> | ||
6 | #include "lapack-aux.h" | ||
7 | |||
8 | #define MACRO(B) do {B} while (0) | ||
9 | #define ERROR(CODE) MACRO(return CODE;) | ||
10 | #define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);}) | ||
11 | |||
12 | #define MIN(A,B) ((A)<(B)?(A):(B)) | ||
13 | #define MAX(A,B) ((A)>(B)?(A):(B)) | ||
14 | |||
15 | // #define DBGL | ||
16 | |||
17 | #ifdef DBGL | ||
18 | #define DEBUGMSG(M) printf("\nLAPACK "M"\n"); | ||
19 | #else | ||
20 | #define DEBUGMSG(M) | ||
21 | #endif | ||
22 | |||
23 | #define OK return 0; | ||
24 | |||
25 | // #ifdef DBGL | ||
26 | // #define DEBUGMSG(M) printf("LAPACK Wrapper "M"\n: "); size_t t0 = time(NULL); | ||
27 | // #define OK MACRO(printf("%ld s\n",time(0)-t0); return 0;); | ||
28 | // #else | ||
29 | // #define DEBUGMSG(M) | ||
30 | // #define OK return 0; | ||
31 | // #endif | ||
32 | |||
33 | #define TRACEMAT(M) {int q; printf(" %d x %d: ",M##r,M##c); \ | ||
34 | for(q=0;q<M##r*M##c;q++) printf("%.1f ",M##p[q]); printf("\n");} | ||
35 | |||
36 | #define CHECK(RES,CODE) MACRO(if(RES) return CODE;) | ||
37 | |||
38 | #define BAD_SIZE 2000 | ||
39 | #define BAD_CODE 2001 | ||
40 | #define MEM 2002 | ||
41 | #define BAD_FILE 2003 | ||
42 | #define SINGULAR 2004 | ||
43 | #define NOCONVER 2005 | ||
44 | #define NODEFPOS 2006 | ||
45 | #define NOSPRTD 2007 | ||
46 | |||
47 | //--------------------------------------- | ||
48 | void asm_finit() { | ||
49 | #ifdef i386 | ||
50 | |||
51 | // asm("finit"); | ||
52 | |||
53 | static unsigned char buf[108]; | ||
54 | asm("FSAVE %0":"=m" (buf)); | ||
55 | |||
56 | #if FPUDEBUG | ||
57 | if(buf[8]!=255 || buf[9]!=255) { // print warning in red | ||
58 | printf("%c[;31mWarning: FPU TAG = %x %x\%c[0m\n",0x1B,buf[8],buf[9],0x1B); | ||
59 | } | ||
60 | #endif | ||
61 | |||
62 | #if NANDEBUG | ||
63 | asm("FRSTOR %0":"=m" (buf)); | ||
64 | #endif | ||
65 | |||
66 | #endif | ||
67 | } | ||
68 | |||
69 | //--------------------------------------- | ||
70 | |||
71 | #if NANDEBUG | ||
72 | |||
73 | #define CHECKNANR(M,msg) \ | ||
74 | { int k; \ | ||
75 | for(k=0; k<(M##r * M##c); k++) { \ | ||
76 | if(M##p[k] != M##p[k]) { \ | ||
77 | printf(msg); \ | ||
78 | TRACEMAT(M) \ | ||
79 | /*exit(1);*/ \ | ||
80 | } \ | ||
81 | } \ | ||
82 | } | ||
83 | |||
84 | #define CHECKNANC(M,msg) \ | ||
85 | { int k; \ | ||
86 | for(k=0; k<(M##r * M##c); k++) { \ | ||
87 | if( M##p[k].r != M##p[k].r \ | ||
88 | || M##p[k].i != M##p[k].i) { \ | ||
89 | printf(msg); \ | ||
90 | /*exit(1);*/ \ | ||
91 | } \ | ||
92 | } \ | ||
93 | } | ||
94 | |||
95 | #else | ||
96 | #define CHECKNANC(M,msg) | ||
97 | #define CHECKNANR(M,msg) | ||
98 | #endif | ||
99 | |||
100 | //--------------------------------------- | ||
101 | |||
102 | //////////////////// real svd //////////////////////////////////// | ||
103 | |||
104 | /* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n, | ||
105 | doublereal *a, integer *lda, doublereal *s, doublereal *u, integer * | ||
106 | ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, | ||
107 | integer *info); | ||
108 | |||
109 | int svd_l_R(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { | ||
110 | integer m = ar; | ||
111 | integer n = ac; | ||
112 | integer q = MIN(m,n); | ||
113 | REQUIRES(sn==q,BAD_SIZE); | ||
114 | REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE); | ||
115 | char* jobu = "A"; | ||
116 | if (up==NULL) { | ||
117 | jobu = "N"; | ||
118 | } else { | ||
119 | if (uc==q) { | ||
120 | jobu = "S"; | ||
121 | } | ||
122 | } | ||
123 | REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE); | ||
124 | char* jobvt = "A"; | ||
125 | integer ldvt = n; | ||
126 | if (vp==NULL) { | ||
127 | jobvt = "N"; | ||
128 | } else { | ||
129 | if (vr==q) { | ||
130 | jobvt = "S"; | ||
131 | ldvt = q; | ||
132 | } | ||
133 | } | ||
134 | DEBUGMSG("svd_l_R"); | ||
135 | double *B = (double*)malloc(m*n*sizeof(double)); | ||
136 | CHECK(!B,MEM); | ||
137 | memcpy(B,ap,m*n*sizeof(double)); | ||
138 | integer lwork = -1; | ||
139 | integer res; | ||
140 | // ask for optimal lwork | ||
141 | double ans; | ||
142 | dgesvd_ (jobu,jobvt, | ||
143 | &m,&n,B,&m, | ||
144 | sp, | ||
145 | up,&m, | ||
146 | vp,&ldvt, | ||
147 | &ans, &lwork, | ||
148 | &res); | ||
149 | lwork = ceil(ans); | ||
150 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
151 | CHECK(!work,MEM); | ||
152 | dgesvd_ (jobu,jobvt, | ||
153 | &m,&n,B,&m, | ||
154 | sp, | ||
155 | up,&m, | ||
156 | vp,&ldvt, | ||
157 | work, &lwork, | ||
158 | &res); | ||
159 | CHECK(res,res); | ||
160 | free(work); | ||
161 | free(B); | ||
162 | OK | ||
163 | } | ||
164 | |||
165 | // (alternative version) | ||
166 | |||
167 | /* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal * | ||
168 | a, integer *lda, doublereal *s, doublereal *u, integer *ldu, | ||
169 | doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, | ||
170 | integer *iwork, integer *info); | ||
171 | |||
172 | int svd_l_Rdd(KDMAT(a),DMAT(u), DVEC(s),DMAT(v)) { | ||
173 | integer m = ar; | ||
174 | integer n = ac; | ||
175 | integer q = MIN(m,n); | ||
176 | REQUIRES(sn==q,BAD_SIZE); | ||
177 | REQUIRES((up == NULL && vp == NULL) | ||
178 | || (ur==m && vc==n | ||
179 | && ((uc == q && vr == q) | ||
180 | || (uc == m && vc==n))),BAD_SIZE); | ||
181 | char* jobz = "A"; | ||
182 | integer ldvt = n; | ||
183 | if (up==NULL) { | ||
184 | jobz = "N"; | ||
185 | } else { | ||
186 | if (uc==q && vr == q) { | ||
187 | jobz = "S"; | ||
188 | ldvt = q; | ||
189 | } | ||
190 | } | ||
191 | DEBUGMSG("svd_l_Rdd"); | ||
192 | double *B = (double*)malloc(m*n*sizeof(double)); | ||
193 | CHECK(!B,MEM); | ||
194 | memcpy(B,ap,m*n*sizeof(double)); | ||
195 | integer* iwk = (integer*) malloc(8*q*sizeof(integer)); | ||
196 | CHECK(!iwk,MEM); | ||
197 | integer lwk = -1; | ||
198 | integer res; | ||
199 | // ask for optimal lwk | ||
200 | double ans; | ||
201 | dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res); | ||
202 | lwk = ans; | ||
203 | double * workv = (double*)malloc(lwk*sizeof(double)); | ||
204 | CHECK(!workv,MEM); | ||
205 | dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res); | ||
206 | CHECK(res,res); | ||
207 | free(iwk); | ||
208 | free(workv); | ||
209 | free(B); | ||
210 | OK | ||
211 | } | ||
212 | |||
213 | //////////////////// complex svd //////////////////////////////////// | ||
214 | |||
215 | // not in clapack.h | ||
216 | |||
217 | int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, | ||
218 | doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, | ||
219 | integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, | ||
220 | integer *lwork, doublereal *rwork, integer *info); | ||
221 | |||
222 | int svd_l_C(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { | ||
223 | integer m = ar; | ||
224 | integer n = ac; | ||
225 | integer q = MIN(m,n); | ||
226 | REQUIRES(sn==q,BAD_SIZE); | ||
227 | REQUIRES(up==NULL || (ur==m && (uc==m || uc==q)),BAD_SIZE); | ||
228 | char* jobu = "A"; | ||
229 | if (up==NULL) { | ||
230 | jobu = "N"; | ||
231 | } else { | ||
232 | if (uc==q) { | ||
233 | jobu = "S"; | ||
234 | } | ||
235 | } | ||
236 | REQUIRES(vp==NULL || (vc==n && (vr==n || vr==q)),BAD_SIZE); | ||
237 | char* jobvt = "A"; | ||
238 | integer ldvt = n; | ||
239 | if (vp==NULL) { | ||
240 | jobvt = "N"; | ||
241 | } else { | ||
242 | if (vr==q) { | ||
243 | jobvt = "S"; | ||
244 | ldvt = q; | ||
245 | } | ||
246 | }DEBUGMSG("svd_l_C"); | ||
247 | doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); | ||
248 | CHECK(!B,MEM); | ||
249 | memcpy(B,ap,m*n*sizeof(doublecomplex)); | ||
250 | |||
251 | double *rwork = (double*) malloc(5*q*sizeof(double)); | ||
252 | CHECK(!rwork,MEM); | ||
253 | integer lwork = -1; | ||
254 | integer res; | ||
255 | // ask for optimal lwork | ||
256 | doublecomplex ans; | ||
257 | zgesvd_ (jobu,jobvt, | ||
258 | &m,&n,B,&m, | ||
259 | sp, | ||
260 | up,&m, | ||
261 | vp,&ldvt, | ||
262 | &ans, &lwork, | ||
263 | rwork, | ||
264 | &res); | ||
265 | lwork = ceil(ans.r); | ||
266 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
267 | CHECK(!work,MEM); | ||
268 | zgesvd_ (jobu,jobvt, | ||
269 | &m,&n,B,&m, | ||
270 | sp, | ||
271 | up,&m, | ||
272 | vp,&ldvt, | ||
273 | work, &lwork, | ||
274 | rwork, | ||
275 | &res); | ||
276 | CHECK(res,res); | ||
277 | free(work); | ||
278 | free(rwork); | ||
279 | free(B); | ||
280 | OK | ||
281 | } | ||
282 | |||
283 | int zgesdd_ (char *jobz, integer *m, integer *n, | ||
284 | doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, | ||
285 | integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, | ||
286 | integer *lwork, doublereal *rwork, integer* iwork, integer *info); | ||
287 | |||
288 | int svd_l_Cdd(KCMAT(a),CMAT(u), DVEC(s),CMAT(v)) { | ||
289 | //printf("entro\n"); | ||
290 | integer m = ar; | ||
291 | integer n = ac; | ||
292 | integer q = MIN(m,n); | ||
293 | REQUIRES(sn==q,BAD_SIZE); | ||
294 | REQUIRES((up == NULL && vp == NULL) | ||
295 | || (ur==m && vc==n | ||
296 | && ((uc == q && vr == q) | ||
297 | || (uc == m && vc==n))),BAD_SIZE); | ||
298 | char* jobz = "A"; | ||
299 | integer ldvt = n; | ||
300 | if (up==NULL) { | ||
301 | jobz = "N"; | ||
302 | } else { | ||
303 | if (uc==q && vr == q) { | ||
304 | jobz = "S"; | ||
305 | ldvt = q; | ||
306 | } | ||
307 | } | ||
308 | DEBUGMSG("svd_l_Cdd"); | ||
309 | doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); | ||
310 | CHECK(!B,MEM); | ||
311 | memcpy(B,ap,m*n*sizeof(doublecomplex)); | ||
312 | integer* iwk = (integer*) malloc(8*q*sizeof(integer)); | ||
313 | CHECK(!iwk,MEM); | ||
314 | int lrwk; | ||
315 | if (0 && *jobz == 'N') { | ||
316 | lrwk = 5*q; // does not work, crash at free below | ||
317 | } else { | ||
318 | lrwk = 5*q*q + 7*q; | ||
319 | } | ||
320 | double *rwk = (double*)malloc(lrwk*sizeof(double));; | ||
321 | CHECK(!rwk,MEM); | ||
322 | //printf("%s %ld %d\n",jobz,q,lrwk); | ||
323 | integer lwk = -1; | ||
324 | integer res; | ||
325 | // ask for optimal lwk | ||
326 | doublecomplex ans; | ||
327 | zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res); | ||
328 | lwk = ans.r; | ||
329 | //printf("lwk = %ld\n",lwk); | ||
330 | doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex)); | ||
331 | CHECK(!workv,MEM); | ||
332 | zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res); | ||
333 | //printf("res = %ld\n",res); | ||
334 | CHECK(res,res); | ||
335 | free(workv); // printf("freed workv\n"); | ||
336 | free(rwk); // printf("freed rwk\n"); | ||
337 | free(iwk); // printf("freed iwk\n"); | ||
338 | free(B); // printf("freed B, salgo\n"); | ||
339 | OK | ||
340 | } | ||
341 | |||
342 | //////////////////// general complex eigensystem //////////// | ||
343 | |||
344 | /* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n, | ||
345 | doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl, | ||
346 | integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work, | ||
347 | integer *lwork, doublereal *rwork, integer *info); | ||
348 | |||
349 | int eig_l_C(KCMAT(a), CMAT(u), CVEC(s),CMAT(v)) { | ||
350 | integer n = ar; | ||
351 | REQUIRES(ac==n && sn==n, BAD_SIZE); | ||
352 | REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); | ||
353 | char jobvl = up==NULL?'N':'V'; | ||
354 | REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE); | ||
355 | char jobvr = vp==NULL?'N':'V'; | ||
356 | DEBUGMSG("eig_l_C"); | ||
357 | doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); | ||
358 | CHECK(!B,MEM); | ||
359 | memcpy(B,ap,n*n*sizeof(doublecomplex)); | ||
360 | double *rwork = (double*) malloc(2*n*sizeof(double)); | ||
361 | CHECK(!rwork,MEM); | ||
362 | integer lwork = -1; | ||
363 | integer res; | ||
364 | // ask for optimal lwork | ||
365 | doublecomplex ans; | ||
366 | //printf("ask zgeev\n"); | ||
367 | zgeev_ (&jobvl,&jobvr, | ||
368 | &n,B,&n, | ||
369 | sp, | ||
370 | up,&n, | ||
371 | vp,&n, | ||
372 | &ans, &lwork, | ||
373 | rwork, | ||
374 | &res); | ||
375 | lwork = ceil(ans.r); | ||
376 | //printf("ans = %d\n",lwork); | ||
377 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
378 | CHECK(!work,MEM); | ||
379 | //printf("zgeev\n"); | ||
380 | zgeev_ (&jobvl,&jobvr, | ||
381 | &n,B,&n, | ||
382 | sp, | ||
383 | up,&n, | ||
384 | vp,&n, | ||
385 | work, &lwork, | ||
386 | rwork, | ||
387 | &res); | ||
388 | CHECK(res,res); | ||
389 | free(work); | ||
390 | free(rwork); | ||
391 | free(B); | ||
392 | OK | ||
393 | } | ||
394 | |||
395 | |||
396 | |||
397 | //////////////////// general real eigensystem //////////// | ||
398 | |||
399 | /* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal * | ||
400 | a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl, | ||
401 | integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work, | ||
402 | integer *lwork, integer *info); | ||
403 | |||
404 | int eig_l_R(KDMAT(a),DMAT(u), CVEC(s),DMAT(v)) { | ||
405 | integer n = ar; | ||
406 | REQUIRES(ac==n && sn==n, BAD_SIZE); | ||
407 | REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE); | ||
408 | char jobvl = up==NULL?'N':'V'; | ||
409 | REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE); | ||
410 | char jobvr = vp==NULL?'N':'V'; | ||
411 | DEBUGMSG("eig_l_R"); | ||
412 | double *B = (double*)malloc(n*n*sizeof(double)); | ||
413 | CHECK(!B,MEM); | ||
414 | memcpy(B,ap,n*n*sizeof(double)); | ||
415 | integer lwork = -1; | ||
416 | integer res; | ||
417 | // ask for optimal lwork | ||
418 | double ans; | ||
419 | //printf("ask dgeev\n"); | ||
420 | dgeev_ (&jobvl,&jobvr, | ||
421 | &n,B,&n, | ||
422 | (double*)sp, (double*)sp+n, | ||
423 | up,&n, | ||
424 | vp,&n, | ||
425 | &ans, &lwork, | ||
426 | &res); | ||
427 | lwork = ceil(ans); | ||
428 | //printf("ans = %d\n",lwork); | ||
429 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
430 | CHECK(!work,MEM); | ||
431 | //printf("dgeev\n"); | ||
432 | dgeev_ (&jobvl,&jobvr, | ||
433 | &n,B,&n, | ||
434 | (double*)sp, (double*)sp+n, | ||
435 | up,&n, | ||
436 | vp,&n, | ||
437 | work, &lwork, | ||
438 | &res); | ||
439 | CHECK(res,res); | ||
440 | free(work); | ||
441 | free(B); | ||
442 | OK | ||
443 | } | ||
444 | |||
445 | |||
446 | //////////////////// symmetric real eigensystem //////////// | ||
447 | |||
448 | /* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a, | ||
449 | integer *lda, doublereal *w, doublereal *work, integer *lwork, | ||
450 | integer *info); | ||
451 | |||
452 | int eig_l_S(int wantV,KDMAT(a),DVEC(s),DMAT(v)) { | ||
453 | integer n = ar; | ||
454 | REQUIRES(ac==n && sn==n, BAD_SIZE); | ||
455 | REQUIRES(vr==n && vc==n, BAD_SIZE); | ||
456 | char jobz = wantV?'V':'N'; | ||
457 | DEBUGMSG("eig_l_S"); | ||
458 | memcpy(vp,ap,n*n*sizeof(double)); | ||
459 | integer lwork = -1; | ||
460 | char uplo = 'U'; | ||
461 | integer res; | ||
462 | // ask for optimal lwork | ||
463 | double ans; | ||
464 | //printf("ask dsyev\n"); | ||
465 | dsyev_ (&jobz,&uplo, | ||
466 | &n,vp,&n, | ||
467 | sp, | ||
468 | &ans, &lwork, | ||
469 | &res); | ||
470 | lwork = ceil(ans); | ||
471 | //printf("ans = %d\n",lwork); | ||
472 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
473 | CHECK(!work,MEM); | ||
474 | dsyev_ (&jobz,&uplo, | ||
475 | &n,vp,&n, | ||
476 | sp, | ||
477 | work, &lwork, | ||
478 | &res); | ||
479 | CHECK(res,res); | ||
480 | free(work); | ||
481 | OK | ||
482 | } | ||
483 | |||
484 | //////////////////// hermitian complex eigensystem //////////// | ||
485 | |||
486 | /* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex | ||
487 | *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork, | ||
488 | doublereal *rwork, integer *info); | ||
489 | |||
490 | int eig_l_H(int wantV,KCMAT(a),DVEC(s),CMAT(v)) { | ||
491 | integer n = ar; | ||
492 | REQUIRES(ac==n && sn==n, BAD_SIZE); | ||
493 | REQUIRES(vr==n && vc==n, BAD_SIZE); | ||
494 | char jobz = wantV?'V':'N'; | ||
495 | DEBUGMSG("eig_l_H"); | ||
496 | memcpy(vp,ap,2*n*n*sizeof(double)); | ||
497 | double *rwork = (double*) malloc((3*n-2)*sizeof(double)); | ||
498 | CHECK(!rwork,MEM); | ||
499 | integer lwork = -1; | ||
500 | char uplo = 'U'; | ||
501 | integer res; | ||
502 | // ask for optimal lwork | ||
503 | doublecomplex ans; | ||
504 | //printf("ask zheev\n"); | ||
505 | zheev_ (&jobz,&uplo, | ||
506 | &n,vp,&n, | ||
507 | sp, | ||
508 | &ans, &lwork, | ||
509 | rwork, | ||
510 | &res); | ||
511 | lwork = ceil(ans.r); | ||
512 | //printf("ans = %d\n",lwork); | ||
513 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
514 | CHECK(!work,MEM); | ||
515 | zheev_ (&jobz,&uplo, | ||
516 | &n,vp,&n, | ||
517 | sp, | ||
518 | work, &lwork, | ||
519 | rwork, | ||
520 | &res); | ||
521 | CHECK(res,res); | ||
522 | free(work); | ||
523 | free(rwork); | ||
524 | OK | ||
525 | } | ||
526 | |||
527 | //////////////////// general real linear system //////////// | ||
528 | |||
529 | /* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer | ||
530 | *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info); | ||
531 | |||
532 | int linearSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) { | ||
533 | integer n = ar; | ||
534 | integer nhrs = bc; | ||
535 | REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); | ||
536 | DEBUGMSG("linearSolveR_l"); | ||
537 | double*AC = (double*)malloc(n*n*sizeof(double)); | ||
538 | memcpy(AC,ap,n*n*sizeof(double)); | ||
539 | memcpy(xp,bp,n*nhrs*sizeof(double)); | ||
540 | integer * ipiv = (integer*)malloc(n*sizeof(integer)); | ||
541 | integer res; | ||
542 | dgesv_ (&n,&nhrs, | ||
543 | AC, &n, | ||
544 | ipiv, | ||
545 | xp, &n, | ||
546 | &res); | ||
547 | if(res>0) { | ||
548 | return SINGULAR; | ||
549 | } | ||
550 | CHECK(res,res); | ||
551 | free(ipiv); | ||
552 | free(AC); | ||
553 | OK | ||
554 | } | ||
555 | |||
556 | //////////////////// general complex linear system //////////// | ||
557 | |||
558 | /* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a, | ||
559 | integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer * | ||
560 | info); | ||
561 | |||
562 | int linearSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { | ||
563 | integer n = ar; | ||
564 | integer nhrs = bc; | ||
565 | REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); | ||
566 | DEBUGMSG("linearSolveC_l"); | ||
567 | doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); | ||
568 | memcpy(AC,ap,n*n*sizeof(doublecomplex)); | ||
569 | memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); | ||
570 | integer * ipiv = (integer*)malloc(n*sizeof(integer)); | ||
571 | integer res; | ||
572 | zgesv_ (&n,&nhrs, | ||
573 | AC, &n, | ||
574 | ipiv, | ||
575 | xp, &n, | ||
576 | &res); | ||
577 | if(res>0) { | ||
578 | return SINGULAR; | ||
579 | } | ||
580 | CHECK(res,res); | ||
581 | free(ipiv); | ||
582 | free(AC); | ||
583 | OK | ||
584 | } | ||
585 | |||
586 | //////// symmetric positive definite real linear system using Cholesky //////////// | ||
587 | |||
588 | /* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs, | ||
589 | doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * | ||
590 | info); | ||
591 | |||
592 | int cholSolveR_l(KDMAT(a),KDMAT(b),DMAT(x)) { | ||
593 | integer n = ar; | ||
594 | integer nhrs = bc; | ||
595 | REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); | ||
596 | DEBUGMSG("cholSolveR_l"); | ||
597 | memcpy(xp,bp,n*nhrs*sizeof(double)); | ||
598 | integer res; | ||
599 | dpotrs_ ("U", | ||
600 | &n,&nhrs, | ||
601 | (double*)ap, &n, | ||
602 | xp, &n, | ||
603 | &res); | ||
604 | CHECK(res,res); | ||
605 | OK | ||
606 | } | ||
607 | |||
608 | //////// Hermitian positive definite real linear system using Cholesky //////////// | ||
609 | |||
610 | /* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, | ||
611 | doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, | ||
612 | integer *info); | ||
613 | |||
614 | int cholSolveC_l(KCMAT(a),KCMAT(b),CMAT(x)) { | ||
615 | integer n = ar; | ||
616 | integer nhrs = bc; | ||
617 | REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); | ||
618 | DEBUGMSG("cholSolveC_l"); | ||
619 | memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); | ||
620 | integer res; | ||
621 | zpotrs_ ("U", | ||
622 | &n,&nhrs, | ||
623 | (doublecomplex*)ap, &n, | ||
624 | xp, &n, | ||
625 | &res); | ||
626 | CHECK(res,res); | ||
627 | OK | ||
628 | } | ||
629 | |||
630 | //////////////////// least squares real linear system //////////// | ||
631 | |||
632 | /* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer * | ||
633 | nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, | ||
634 | doublereal *work, integer *lwork, integer *info); | ||
635 | |||
636 | int linearSolveLSR_l(KDMAT(a),KDMAT(b),DMAT(x)) { | ||
637 | integer m = ar; | ||
638 | integer n = ac; | ||
639 | integer nrhs = bc; | ||
640 | integer ldb = xr; | ||
641 | REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); | ||
642 | DEBUGMSG("linearSolveLSR_l"); | ||
643 | double*AC = (double*)malloc(m*n*sizeof(double)); | ||
644 | memcpy(AC,ap,m*n*sizeof(double)); | ||
645 | if (m>=n) { | ||
646 | memcpy(xp,bp,m*nrhs*sizeof(double)); | ||
647 | } else { | ||
648 | int k; | ||
649 | for(k = 0; k<nrhs; k++) { | ||
650 | memcpy(xp+ldb*k,bp+m*k,m*sizeof(double)); | ||
651 | } | ||
652 | } | ||
653 | integer res; | ||
654 | integer lwork = -1; | ||
655 | double ans; | ||
656 | dgels_ ("N",&m,&n,&nrhs, | ||
657 | AC,&m, | ||
658 | xp,&ldb, | ||
659 | &ans,&lwork, | ||
660 | &res); | ||
661 | lwork = ceil(ans); | ||
662 | //printf("ans = %d\n",lwork); | ||
663 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
664 | dgels_ ("N",&m,&n,&nrhs, | ||
665 | AC,&m, | ||
666 | xp,&ldb, | ||
667 | work,&lwork, | ||
668 | &res); | ||
669 | if(res>0) { | ||
670 | return SINGULAR; | ||
671 | } | ||
672 | CHECK(res,res); | ||
673 | free(work); | ||
674 | free(AC); | ||
675 | OK | ||
676 | } | ||
677 | |||
678 | //////////////////// least squares complex linear system //////////// | ||
679 | |||
680 | /* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer * | ||
681 | nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, | ||
682 | doublecomplex *work, integer *lwork, integer *info); | ||
683 | |||
684 | int linearSolveLSC_l(KCMAT(a),KCMAT(b),CMAT(x)) { | ||
685 | integer m = ar; | ||
686 | integer n = ac; | ||
687 | integer nrhs = bc; | ||
688 | integer ldb = xr; | ||
689 | REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); | ||
690 | DEBUGMSG("linearSolveLSC_l"); | ||
691 | doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); | ||
692 | memcpy(AC,ap,m*n*sizeof(doublecomplex)); | ||
693 | if (m>=n) { | ||
694 | memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); | ||
695 | } else { | ||
696 | int k; | ||
697 | for(k = 0; k<nrhs; k++) { | ||
698 | memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex)); | ||
699 | } | ||
700 | } | ||
701 | integer res; | ||
702 | integer lwork = -1; | ||
703 | doublecomplex ans; | ||
704 | zgels_ ("N",&m,&n,&nrhs, | ||
705 | AC,&m, | ||
706 | xp,&ldb, | ||
707 | &ans,&lwork, | ||
708 | &res); | ||
709 | lwork = ceil(ans.r); | ||
710 | //printf("ans = %d\n",lwork); | ||
711 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
712 | zgels_ ("N",&m,&n,&nrhs, | ||
713 | AC,&m, | ||
714 | xp,&ldb, | ||
715 | work,&lwork, | ||
716 | &res); | ||
717 | if(res>0) { | ||
718 | return SINGULAR; | ||
719 | } | ||
720 | CHECK(res,res); | ||
721 | free(work); | ||
722 | free(AC); | ||
723 | OK | ||
724 | } | ||
725 | |||
726 | //////////////////// least squares real linear system using SVD //////////// | ||
727 | |||
728 | /* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, | ||
729 | doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * | ||
730 | s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, | ||
731 | integer *info); | ||
732 | |||
733 | int linearSolveSVDR_l(double rcond,KDMAT(a),KDMAT(b),DMAT(x)) { | ||
734 | integer m = ar; | ||
735 | integer n = ac; | ||
736 | integer nrhs = bc; | ||
737 | integer ldb = xr; | ||
738 | REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); | ||
739 | DEBUGMSG("linearSolveSVDR_l"); | ||
740 | double*AC = (double*)malloc(m*n*sizeof(double)); | ||
741 | double*S = (double*)malloc(MIN(m,n)*sizeof(double)); | ||
742 | memcpy(AC,ap,m*n*sizeof(double)); | ||
743 | if (m>=n) { | ||
744 | memcpy(xp,bp,m*nrhs*sizeof(double)); | ||
745 | } else { | ||
746 | int k; | ||
747 | for(k = 0; k<nrhs; k++) { | ||
748 | memcpy(xp+ldb*k,bp+m*k,m*sizeof(double)); | ||
749 | } | ||
750 | } | ||
751 | integer res; | ||
752 | integer lwork = -1; | ||
753 | integer rank; | ||
754 | double ans; | ||
755 | dgelss_ (&m,&n,&nrhs, | ||
756 | AC,&m, | ||
757 | xp,&ldb, | ||
758 | S, | ||
759 | &rcond,&rank, | ||
760 | &ans,&lwork, | ||
761 | &res); | ||
762 | lwork = ceil(ans); | ||
763 | //printf("ans = %d\n",lwork); | ||
764 | double * work = (double*)malloc(lwork*sizeof(double)); | ||
765 | dgelss_ (&m,&n,&nrhs, | ||
766 | AC,&m, | ||
767 | xp,&ldb, | ||
768 | S, | ||
769 | &rcond,&rank, | ||
770 | work,&lwork, | ||
771 | &res); | ||
772 | if(res>0) { | ||
773 | return NOCONVER; | ||
774 | } | ||
775 | CHECK(res,res); | ||
776 | free(work); | ||
777 | free(S); | ||
778 | free(AC); | ||
779 | OK | ||
780 | } | ||
781 | |||
782 | //////////////////// least squares complex linear system using SVD //////////// | ||
783 | |||
784 | // not in clapack.h | ||
785 | |||
786 | int zgelss_(integer *m, integer *n, integer *nhrs, | ||
787 | doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s, | ||
788 | doublereal *rcond, integer* rank, | ||
789 | doublecomplex *work, integer* lwork, doublereal* rwork, | ||
790 | integer *info); | ||
791 | |||
792 | int linearSolveSVDC_l(double rcond, KCMAT(a),KCMAT(b),CMAT(x)) { | ||
793 | integer m = ar; | ||
794 | integer n = ac; | ||
795 | integer nrhs = bc; | ||
796 | integer ldb = xr; | ||
797 | REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); | ||
798 | DEBUGMSG("linearSolveSVDC_l"); | ||
799 | doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); | ||
800 | double*S = (double*)malloc(MIN(m,n)*sizeof(double)); | ||
801 | double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); | ||
802 | memcpy(AC,ap,m*n*sizeof(doublecomplex)); | ||
803 | if (m>=n) { | ||
804 | memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); | ||
805 | } else { | ||
806 | int k; | ||
807 | for(k = 0; k<nrhs; k++) { | ||
808 | memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex)); | ||
809 | } | ||
810 | } | ||
811 | integer res; | ||
812 | integer lwork = -1; | ||
813 | integer rank; | ||
814 | doublecomplex ans; | ||
815 | zgelss_ (&m,&n,&nrhs, | ||
816 | AC,&m, | ||
817 | xp,&ldb, | ||
818 | S, | ||
819 | &rcond,&rank, | ||
820 | &ans,&lwork, | ||
821 | RWORK, | ||
822 | &res); | ||
823 | lwork = ceil(ans.r); | ||
824 | //printf("ans = %d\n",lwork); | ||
825 | doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
826 | zgelss_ (&m,&n,&nrhs, | ||
827 | AC,&m, | ||
828 | xp,&ldb, | ||
829 | S, | ||
830 | &rcond,&rank, | ||
831 | work,&lwork, | ||
832 | RWORK, | ||
833 | &res); | ||
834 | if(res>0) { | ||
835 | return NOCONVER; | ||
836 | } | ||
837 | CHECK(res,res); | ||
838 | free(work); | ||
839 | free(RWORK); | ||
840 | free(S); | ||
841 | free(AC); | ||
842 | OK | ||
843 | } | ||
844 | |||
845 | //////////////////// Cholesky factorization ///////////////////////// | ||
846 | |||
847 | /* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, | ||
848 | integer *lda, integer *info); | ||
849 | |||
850 | int chol_l_H(KCMAT(a),CMAT(l)) { | ||
851 | integer n = ar; | ||
852 | REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE); | ||
853 | DEBUGMSG("chol_l_H"); | ||
854 | memcpy(lp,ap,n*n*sizeof(doublecomplex)); | ||
855 | char uplo = 'U'; | ||
856 | integer res; | ||
857 | zpotrf_ (&uplo,&n,lp,&n,&res); | ||
858 | CHECK(res>0,NODEFPOS); | ||
859 | CHECK(res,res); | ||
860 | doublecomplex zero = {0.,0.}; | ||
861 | int r,c; | ||
862 | for (r=0; r<lr-1; r++) { | ||
863 | for(c=r+1; c<lc; c++) { | ||
864 | lp[r*lc+c] = zero; | ||
865 | } | ||
866 | } | ||
867 | OK | ||
868 | } | ||
869 | |||
870 | |||
871 | /* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer * | ||
872 | lda, integer *info); | ||
873 | |||
874 | int chol_l_S(KDMAT(a),DMAT(l)) { | ||
875 | integer n = ar; | ||
876 | REQUIRES(n>=1 && ac == n && lr==n && lc==n,BAD_SIZE); | ||
877 | DEBUGMSG("chol_l_S"); | ||
878 | memcpy(lp,ap,n*n*sizeof(double)); | ||
879 | char uplo = 'U'; | ||
880 | integer res; | ||
881 | dpotrf_ (&uplo,&n,lp,&n,&res); | ||
882 | CHECK(res>0,NODEFPOS); | ||
883 | CHECK(res,res); | ||
884 | int r,c; | ||
885 | for (r=0; r<lr-1; r++) { | ||
886 | for(c=r+1; c<lc; c++) { | ||
887 | lp[r*lc+c] = 0.; | ||
888 | } | ||
889 | } | ||
890 | OK | ||
891 | } | ||
892 | |||
893 | //////////////////// QR factorization ///////////////////////// | ||
894 | |||
895 | /* Subroutine */ int dgeqr2_(integer *m, integer *n, doublereal *a, integer * | ||
896 | lda, doublereal *tau, doublereal *work, integer *info); | ||
897 | |||
898 | int qr_l_R(KDMAT(a), DVEC(tau), DMAT(r)) { | ||
899 | integer m = ar; | ||
900 | integer n = ac; | ||
901 | integer mn = MIN(m,n); | ||
902 | REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); | ||
903 | DEBUGMSG("qr_l_R"); | ||
904 | double *WORK = (double*)malloc(n*sizeof(double)); | ||
905 | CHECK(!WORK,MEM); | ||
906 | memcpy(rp,ap,m*n*sizeof(double)); | ||
907 | integer res; | ||
908 | dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); | ||
909 | CHECK(res,res); | ||
910 | free(WORK); | ||
911 | OK | ||
912 | } | ||
913 | |||
914 | /* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a, | ||
915 | integer *lda, doublecomplex *tau, doublecomplex *work, integer *info); | ||
916 | |||
917 | int qr_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { | ||
918 | integer m = ar; | ||
919 | integer n = ac; | ||
920 | integer mn = MIN(m,n); | ||
921 | REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); | ||
922 | DEBUGMSG("qr_l_C"); | ||
923 | doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex)); | ||
924 | CHECK(!WORK,MEM); | ||
925 | memcpy(rp,ap,m*n*sizeof(doublecomplex)); | ||
926 | integer res; | ||
927 | zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); | ||
928 | CHECK(res,res); | ||
929 | free(WORK); | ||
930 | OK | ||
931 | } | ||
932 | |||
933 | /* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal * | ||
934 | a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, | ||
935 | integer *info); | ||
936 | |||
937 | int c_dorgqr(KDMAT(a), KDVEC(tau), DMAT(r)) { | ||
938 | integer m = ar; | ||
939 | integer n = MIN(ac,ar); | ||
940 | integer k = taun; | ||
941 | DEBUGMSG("c_dorgqr"); | ||
942 | integer lwork = 8*n; // FIXME | ||
943 | double *WORK = (double*)malloc(lwork*sizeof(double)); | ||
944 | CHECK(!WORK,MEM); | ||
945 | memcpy(rp,ap,m*k*sizeof(double)); | ||
946 | integer res; | ||
947 | dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res); | ||
948 | CHECK(res,res); | ||
949 | free(WORK); | ||
950 | OK | ||
951 | } | ||
952 | |||
953 | /* Subroutine */ int zungqr_(integer *m, integer *n, integer *k, | ||
954 | doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * | ||
955 | work, integer *lwork, integer *info); | ||
956 | |||
957 | int c_zungqr(KCMAT(a), KCVEC(tau), CMAT(r)) { | ||
958 | integer m = ar; | ||
959 | integer n = MIN(ac,ar); | ||
960 | integer k = taun; | ||
961 | DEBUGMSG("z_ungqr"); | ||
962 | integer lwork = 8*n; // FIXME | ||
963 | doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
964 | CHECK(!WORK,MEM); | ||
965 | memcpy(rp,ap,m*k*sizeof(doublecomplex)); | ||
966 | integer res; | ||
967 | zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res); | ||
968 | CHECK(res,res); | ||
969 | free(WORK); | ||
970 | OK | ||
971 | } | ||
972 | |||
973 | |||
974 | //////////////////// Hessenberg factorization ///////////////////////// | ||
975 | |||
976 | /* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, | ||
977 | doublereal *a, integer *lda, doublereal *tau, doublereal *work, | ||
978 | integer *lwork, integer *info); | ||
979 | |||
980 | int hess_l_R(KDMAT(a), DVEC(tau), DMAT(r)) { | ||
981 | integer m = ar; | ||
982 | integer n = ac; | ||
983 | integer mn = MIN(m,n); | ||
984 | REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); | ||
985 | DEBUGMSG("hess_l_R"); | ||
986 | integer lwork = 5*n; // fixme | ||
987 | double *WORK = (double*)malloc(lwork*sizeof(double)); | ||
988 | CHECK(!WORK,MEM); | ||
989 | memcpy(rp,ap,m*n*sizeof(double)); | ||
990 | integer res; | ||
991 | integer one = 1; | ||
992 | dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); | ||
993 | CHECK(res,res); | ||
994 | free(WORK); | ||
995 | OK | ||
996 | } | ||
997 | |||
998 | |||
999 | /* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi, | ||
1000 | doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * | ||
1001 | work, integer *lwork, integer *info); | ||
1002 | |||
1003 | int hess_l_C(KCMAT(a), CVEC(tau), CMAT(r)) { | ||
1004 | integer m = ar; | ||
1005 | integer n = ac; | ||
1006 | integer mn = MIN(m,n); | ||
1007 | REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); | ||
1008 | DEBUGMSG("hess_l_C"); | ||
1009 | integer lwork = 5*n; // fixme | ||
1010 | doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
1011 | CHECK(!WORK,MEM); | ||
1012 | memcpy(rp,ap,m*n*sizeof(doublecomplex)); | ||
1013 | integer res; | ||
1014 | integer one = 1; | ||
1015 | zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); | ||
1016 | CHECK(res,res); | ||
1017 | free(WORK); | ||
1018 | OK | ||
1019 | } | ||
1020 | |||
1021 | //////////////////// Schur factorization ///////////////////////// | ||
1022 | |||
1023 | /* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n, | ||
1024 | doublereal *a, integer *lda, integer *sdim, doublereal *wr, | ||
1025 | doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, | ||
1026 | integer *lwork, logical *bwork, integer *info); | ||
1027 | |||
1028 | int schur_l_R(KDMAT(a), DMAT(u), DMAT(s)) { | ||
1029 | integer m = ar; | ||
1030 | integer n = ac; | ||
1031 | REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); | ||
1032 | DEBUGMSG("schur_l_R"); | ||
1033 | //int k; | ||
1034 | //printf("---------------------------\n"); | ||
1035 | //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n"); | ||
1036 | //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n"); | ||
1037 | //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n"); | ||
1038 | memcpy(sp,ap,n*n*sizeof(double)); | ||
1039 | integer lwork = 6*n; // fixme | ||
1040 | double *WORK = (double*)malloc(lwork*sizeof(double)); | ||
1041 | double *WR = (double*)malloc(n*sizeof(double)); | ||
1042 | double *WI = (double*)malloc(n*sizeof(double)); | ||
1043 | // WR and WI not really required in this call | ||
1044 | logical *BWORK = (logical*)malloc(n*sizeof(logical)); | ||
1045 | integer res; | ||
1046 | integer sdim; | ||
1047 | dgees_ ("V","N",NULL,&n,sp,&n,&sdim,WR,WI,up,&n,WORK,&lwork,BWORK,&res); | ||
1048 | //printf("%p: ",ap); for(k=0;k<n*n;k++) printf("%f ",ap[k]); printf("\n"); | ||
1049 | //printf("%p: ",up); for(k=0;k<n*n;k++) printf("%f ",up[k]); printf("\n"); | ||
1050 | //printf("%p: ",sp); for(k=0;k<n*n;k++) printf("%f ",sp[k]); printf("\n"); | ||
1051 | if(res>0) { | ||
1052 | return NOCONVER; | ||
1053 | } | ||
1054 | CHECK(res,res); | ||
1055 | free(WR); | ||
1056 | free(WI); | ||
1057 | free(BWORK); | ||
1058 | free(WORK); | ||
1059 | OK | ||
1060 | } | ||
1061 | |||
1062 | |||
1063 | /* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n, | ||
1064 | doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w, | ||
1065 | doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork, | ||
1066 | doublereal *rwork, logical *bwork, integer *info); | ||
1067 | |||
1068 | int schur_l_C(KCMAT(a), CMAT(u), CMAT(s)) { | ||
1069 | integer m = ar; | ||
1070 | integer n = ac; | ||
1071 | REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); | ||
1072 | DEBUGMSG("schur_l_C"); | ||
1073 | memcpy(sp,ap,n*n*sizeof(doublecomplex)); | ||
1074 | integer lwork = 6*n; // fixme | ||
1075 | doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); | ||
1076 | doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); | ||
1077 | // W not really required in this call | ||
1078 | logical *BWORK = (logical*)malloc(n*sizeof(logical)); | ||
1079 | double *RWORK = (double*)malloc(n*sizeof(double)); | ||
1080 | integer res; | ||
1081 | integer sdim; | ||
1082 | zgees_ ("V","N",NULL,&n,sp,&n,&sdim,W, | ||
1083 | up,&n, | ||
1084 | WORK,&lwork,RWORK,BWORK,&res); | ||
1085 | if(res>0) { | ||
1086 | return NOCONVER; | ||
1087 | } | ||
1088 | CHECK(res,res); | ||
1089 | free(W); | ||
1090 | free(BWORK); | ||
1091 | free(WORK); | ||
1092 | OK | ||
1093 | } | ||
1094 | |||
1095 | //////////////////// LU factorization ///////////////////////// | ||
1096 | |||
1097 | /* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer * | ||
1098 | lda, integer *ipiv, integer *info); | ||
1099 | |||
1100 | int lu_l_R(KDMAT(a), DVEC(ipiv), DMAT(r)) { | ||
1101 | integer m = ar; | ||
1102 | integer n = ac; | ||
1103 | integer mn = MIN(m,n); | ||
1104 | REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); | ||
1105 | DEBUGMSG("lu_l_R"); | ||
1106 | integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); | ||
1107 | memcpy(rp,ap,m*n*sizeof(double)); | ||
1108 | integer res; | ||
1109 | dgetrf_ (&m,&n,rp,&m,auxipiv,&res); | ||
1110 | if(res>0) { | ||
1111 | res = 0; // fixme | ||
1112 | } | ||
1113 | CHECK(res,res); | ||
1114 | int k; | ||
1115 | for (k=0; k<mn; k++) { | ||
1116 | ipivp[k] = auxipiv[k]; | ||
1117 | } | ||
1118 | free(auxipiv); | ||
1119 | OK | ||
1120 | } | ||
1121 | |||
1122 | |||
1123 | /* Subroutine */ int zgetrf_(integer *m, integer *n, doublecomplex *a, | ||
1124 | integer *lda, integer *ipiv, integer *info); | ||
1125 | |||
1126 | int lu_l_C(KCMAT(a), DVEC(ipiv), CMAT(r)) { | ||
1127 | integer m = ar; | ||
1128 | integer n = ac; | ||
1129 | integer mn = MIN(m,n); | ||
1130 | REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); | ||
1131 | DEBUGMSG("lu_l_C"); | ||
1132 | integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); | ||
1133 | memcpy(rp,ap,m*n*sizeof(doublecomplex)); | ||
1134 | integer res; | ||
1135 | zgetrf_ (&m,&n,rp,&m,auxipiv,&res); | ||
1136 | if(res>0) { | ||
1137 | res = 0; // fixme | ||
1138 | } | ||
1139 | CHECK(res,res); | ||
1140 | int k; | ||
1141 | for (k=0; k<mn; k++) { | ||
1142 | ipivp[k] = auxipiv[k]; | ||
1143 | } | ||
1144 | free(auxipiv); | ||
1145 | OK | ||
1146 | } | ||
1147 | |||
1148 | |||
1149 | //////////////////// LU substitution ///////////////////////// | ||
1150 | |||
1151 | /* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs, | ||
1152 | doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer * | ||
1153 | ldb, integer *info); | ||
1154 | |||
1155 | int luS_l_R(KDMAT(a), KDVEC(ipiv), KDMAT(b), DMAT(x)) { | ||
1156 | integer m = ar; | ||
1157 | integer n = ac; | ||
1158 | integer mrhs = br; | ||
1159 | integer nrhs = bc; | ||
1160 | |||
1161 | REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE); | ||
1162 | integer* auxipiv = (integer*)malloc(n*sizeof(integer)); | ||
1163 | int k; | ||
1164 | for (k=0; k<n; k++) { | ||
1165 | auxipiv[k] = (integer)ipivp[k]; | ||
1166 | } | ||
1167 | integer res; | ||
1168 | memcpy(xp,bp,mrhs*nrhs*sizeof(double)); | ||
1169 | dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res); | ||
1170 | CHECK(res,res); | ||
1171 | free(auxipiv); | ||
1172 | OK | ||
1173 | } | ||
1174 | |||
1175 | |||
1176 | /* Subroutine */ int zgetrs_(char *trans, integer *n, integer *nrhs, | ||
1177 | doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, | ||
1178 | integer *ldb, integer *info); | ||
1179 | |||
1180 | int luS_l_C(KCMAT(a), KDVEC(ipiv), KCMAT(b), CMAT(x)) { | ||
1181 | integer m = ar; | ||
1182 | integer n = ac; | ||
1183 | integer mrhs = br; | ||
1184 | integer nrhs = bc; | ||
1185 | |||
1186 | REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE); | ||
1187 | integer* auxipiv = (integer*)malloc(n*sizeof(integer)); | ||
1188 | int k; | ||
1189 | for (k=0; k<n; k++) { | ||
1190 | auxipiv[k] = (integer)ipivp[k]; | ||
1191 | } | ||
1192 | integer res; | ||
1193 | memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex)); | ||
1194 | zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res); | ||
1195 | CHECK(res,res); | ||
1196 | free(auxipiv); | ||
1197 | OK | ||
1198 | } | ||
1199 | |||
1200 | //////////////////// Matrix Product ///////////////////////// | ||
1201 | |||
1202 | void dgemm_(char *, char *, integer *, integer *, integer *, | ||
1203 | double *, const double *, integer *, const double *, | ||
1204 | integer *, double *, double *, integer *); | ||
1205 | |||
1206 | int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { | ||
1207 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
1208 | DEBUGMSG("dgemm_"); | ||
1209 | CHECKNANR(a,"NaN multR Input\n") | ||
1210 | CHECKNANR(b,"NaN multR Input\n") | ||
1211 | integer m = ta?ac:ar; | ||
1212 | integer n = tb?br:bc; | ||
1213 | integer k = ta?ar:ac; | ||
1214 | integer lda = ar; | ||
1215 | integer ldb = br; | ||
1216 | integer ldc = rr; | ||
1217 | double alpha = 1; | ||
1218 | double beta = 0; | ||
1219 | dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); | ||
1220 | CHECKNANR(r,"NaN multR Output\n") | ||
1221 | OK | ||
1222 | } | ||
1223 | |||
1224 | void zgemm_(char *, char *, integer *, integer *, integer *, | ||
1225 | doublecomplex *, const doublecomplex *, integer *, const doublecomplex *, | ||
1226 | integer *, doublecomplex *, doublecomplex *, integer *); | ||
1227 | |||
1228 | int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { | ||
1229 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
1230 | DEBUGMSG("zgemm_"); | ||
1231 | CHECKNANC(a,"NaN multC Input\n") | ||
1232 | CHECKNANC(b,"NaN multC Input\n") | ||
1233 | integer m = ta?ac:ar; | ||
1234 | integer n = tb?br:bc; | ||
1235 | integer k = ta?ar:ac; | ||
1236 | integer lda = ar; | ||
1237 | integer ldb = br; | ||
1238 | integer ldc = rr; | ||
1239 | doublecomplex alpha = {1,0}; | ||
1240 | doublecomplex beta = {0,0}; | ||
1241 | zgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha, | ||
1242 | ap,&lda, | ||
1243 | bp,&ldb,&beta, | ||
1244 | rp,&ldc); | ||
1245 | CHECKNANC(r,"NaN multC Output\n") | ||
1246 | OK | ||
1247 | } | ||
1248 | |||
1249 | void sgemm_(char *, char *, integer *, integer *, integer *, | ||
1250 | float *, const float *, integer *, const float *, | ||
1251 | integer *, float *, float *, integer *); | ||
1252 | |||
1253 | int multiplyF(int ta, int tb, KFMAT(a),KFMAT(b),FMAT(r)) { | ||
1254 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
1255 | DEBUGMSG("sgemm_"); | ||
1256 | integer m = ta?ac:ar; | ||
1257 | integer n = tb?br:bc; | ||
1258 | integer k = ta?ar:ac; | ||
1259 | integer lda = ar; | ||
1260 | integer ldb = br; | ||
1261 | integer ldc = rr; | ||
1262 | float alpha = 1; | ||
1263 | float beta = 0; | ||
1264 | sgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); | ||
1265 | OK | ||
1266 | } | ||
1267 | |||
1268 | void cgemm_(char *, char *, integer *, integer *, integer *, | ||
1269 | complex *, const complex *, integer *, const complex *, | ||
1270 | integer *, complex *, complex *, integer *); | ||
1271 | |||
1272 | int multiplyQ(int ta, int tb, KQMAT(a),KQMAT(b),QMAT(r)) { | ||
1273 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
1274 | DEBUGMSG("cgemm_"); | ||
1275 | integer m = ta?ac:ar; | ||
1276 | integer n = tb?br:bc; | ||
1277 | integer k = ta?ar:ac; | ||
1278 | integer lda = ar; | ||
1279 | integer ldb = br; | ||
1280 | integer ldc = rr; | ||
1281 | complex alpha = {1,0}; | ||
1282 | complex beta = {0,0}; | ||
1283 | cgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha, | ||
1284 | ap,&lda, | ||
1285 | bp,&ldb,&beta, | ||
1286 | rp,&ldc); | ||
1287 | OK | ||
1288 | } | ||
1289 | |||
1290 | //////////////////// transpose ///////////////////////// | ||
1291 | |||
1292 | int transF(KFMAT(x),FMAT(t)) { | ||
1293 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1294 | DEBUGMSG("transF"); | ||
1295 | int i,j; | ||
1296 | for (i=0; i<tr; i++) { | ||
1297 | for (j=0; j<tc; j++) { | ||
1298 | tp[i*tc+j] = xp[j*xc+i]; | ||
1299 | } | ||
1300 | } | ||
1301 | OK | ||
1302 | } | ||
1303 | |||
1304 | int transR(KDMAT(x),DMAT(t)) { | ||
1305 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1306 | DEBUGMSG("transR"); | ||
1307 | int i,j; | ||
1308 | for (i=0; i<tr; i++) { | ||
1309 | for (j=0; j<tc; j++) { | ||
1310 | tp[i*tc+j] = xp[j*xc+i]; | ||
1311 | } | ||
1312 | } | ||
1313 | OK | ||
1314 | } | ||
1315 | |||
1316 | int transQ(KQMAT(x),QMAT(t)) { | ||
1317 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1318 | DEBUGMSG("transQ"); | ||
1319 | int i,j; | ||
1320 | for (i=0; i<tr; i++) { | ||
1321 | for (j=0; j<tc; j++) { | ||
1322 | tp[i*tc+j] = xp[j*xc+i]; | ||
1323 | } | ||
1324 | } | ||
1325 | OK | ||
1326 | } | ||
1327 | |||
1328 | int transC(KCMAT(x),CMAT(t)) { | ||
1329 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1330 | DEBUGMSG("transC"); | ||
1331 | int i,j; | ||
1332 | for (i=0; i<tr; i++) { | ||
1333 | for (j=0; j<tc; j++) { | ||
1334 | tp[i*tc+j] = xp[j*xc+i]; | ||
1335 | } | ||
1336 | } | ||
1337 | OK | ||
1338 | } | ||
1339 | |||
1340 | int transP(KPMAT(x), PMAT(t)) { | ||
1341 | REQUIRES(xr==tc && xc==tr,BAD_SIZE); | ||
1342 | REQUIRES(xs==ts,NOCONVER); | ||
1343 | DEBUGMSG("transP"); | ||
1344 | int i,j; | ||
1345 | for (i=0; i<tr; i++) { | ||
1346 | for (j=0; j<tc; j++) { | ||
1347 | memcpy(tp+(i*tc+j)*xs,xp +(j*xc+i)*xs,xs); | ||
1348 | } | ||
1349 | } | ||
1350 | OK | ||
1351 | } | ||
1352 | |||
1353 | //////////////////// constant ///////////////////////// | ||
1354 | |||
1355 | int constantF(float * pval, FVEC(r)) { | ||
1356 | DEBUGMSG("constantF") | ||
1357 | int k; | ||
1358 | double val = *pval; | ||
1359 | for(k=0;k<rn;k++) { | ||
1360 | rp[k]=val; | ||
1361 | } | ||
1362 | OK | ||
1363 | } | ||
1364 | |||
1365 | int constantR(double * pval, DVEC(r)) { | ||
1366 | DEBUGMSG("constantR") | ||
1367 | int k; | ||
1368 | double val = *pval; | ||
1369 | for(k=0;k<rn;k++) { | ||
1370 | rp[k]=val; | ||
1371 | } | ||
1372 | OK | ||
1373 | } | ||
1374 | |||
1375 | int constantQ(complex* pval, QVEC(r)) { | ||
1376 | DEBUGMSG("constantQ") | ||
1377 | int k; | ||
1378 | complex val = *pval; | ||
1379 | for(k=0;k<rn;k++) { | ||
1380 | rp[k]=val; | ||
1381 | } | ||
1382 | OK | ||
1383 | } | ||
1384 | |||
1385 | int constantC(doublecomplex* pval, CVEC(r)) { | ||
1386 | DEBUGMSG("constantC") | ||
1387 | int k; | ||
1388 | doublecomplex val = *pval; | ||
1389 | for(k=0;k<rn;k++) { | ||
1390 | rp[k]=val; | ||
1391 | } | ||
1392 | OK | ||
1393 | } | ||
1394 | |||
1395 | int constantP(void* pval, PVEC(r)) { | ||
1396 | DEBUGMSG("constantP") | ||
1397 | int k; | ||
1398 | for(k=0;k<rn;k++) { | ||
1399 | memcpy(rp+k*rs,pval,rs); | ||
1400 | } | ||
1401 | OK | ||
1402 | } | ||
1403 | |||
1404 | //////////////////// float-double conversion ///////////////////////// | ||
1405 | |||
1406 | int float2double(FVEC(x),DVEC(y)) { | ||
1407 | DEBUGMSG("float2double") | ||
1408 | int k; | ||
1409 | for(k=0;k<xn;k++) { | ||
1410 | yp[k]=xp[k]; | ||
1411 | } | ||
1412 | OK | ||
1413 | } | ||
1414 | |||
1415 | int double2float(DVEC(x),FVEC(y)) { | ||
1416 | DEBUGMSG("double2float") | ||
1417 | int k; | ||
1418 | for(k=0;k<xn;k++) { | ||
1419 | yp[k]=xp[k]; | ||
1420 | } | ||
1421 | OK | ||
1422 | } | ||
1423 | |||
1424 | //////////////////// conjugate ///////////////////////// | ||
1425 | |||
1426 | int conjugateQ(KQVEC(x),QVEC(t)) { | ||
1427 | REQUIRES(xn==tn,BAD_SIZE); | ||
1428 | DEBUGMSG("conjugateQ"); | ||
1429 | int k; | ||
1430 | for(k=0;k<xn;k++) { | ||
1431 | tp[k].r = xp[k].r; | ||
1432 | tp[k].i = -xp[k].i; | ||
1433 | } | ||
1434 | OK | ||
1435 | } | ||
1436 | |||
1437 | int conjugateC(KCVEC(x),CVEC(t)) { | ||
1438 | REQUIRES(xn==tn,BAD_SIZE); | ||
1439 | DEBUGMSG("conjugateC"); | ||
1440 | int k; | ||
1441 | for(k=0;k<xn;k++) { | ||
1442 | tp[k].r = xp[k].r; | ||
1443 | tp[k].i = -xp[k].i; | ||
1444 | } | ||
1445 | OK | ||
1446 | } | ||
1447 | |||
1448 | //////////////////// step ///////////////////////// | ||
1449 | |||
1450 | int stepF(FVEC(x),FVEC(y)) { | ||
1451 | DEBUGMSG("stepF") | ||
1452 | int k; | ||
1453 | for(k=0;k<xn;k++) { | ||
1454 | yp[k]=xp[k]>0; | ||
1455 | } | ||
1456 | OK | ||
1457 | } | ||
1458 | |||
1459 | int stepD(DVEC(x),DVEC(y)) { | ||
1460 | DEBUGMSG("stepD") | ||
1461 | int k; | ||
1462 | for(k=0;k<xn;k++) { | ||
1463 | yp[k]=xp[k]>0; | ||
1464 | } | ||
1465 | OK | ||
1466 | } | ||
1467 | |||
1468 | //////////////////// cond ///////////////////////// | ||
1469 | |||
1470 | int condF(FVEC(x),FVEC(y),FVEC(lt),FVEC(eq),FVEC(gt),FVEC(r)) { | ||
1471 | REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); | ||
1472 | DEBUGMSG("condF") | ||
1473 | int k; | ||
1474 | for(k=0;k<xn;k++) { | ||
1475 | rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]); | ||
1476 | } | ||
1477 | OK | ||
1478 | } | ||
1479 | |||
1480 | int condD(DVEC(x),DVEC(y),DVEC(lt),DVEC(eq),DVEC(gt),DVEC(r)) { | ||
1481 | REQUIRES(xn==yn && xn==ltn && xn==eqn && xn==gtn && xn==rn ,BAD_SIZE); | ||
1482 | DEBUGMSG("condD") | ||
1483 | int k; | ||
1484 | for(k=0;k<xn;k++) { | ||
1485 | rp[k] = xp[k]<yp[k]?ltp[k]:(xp[k]>yp[k]?gtp[k]:eqp[k]); | ||
1486 | } | ||
1487 | OK | ||
1488 | } | ||
1489 | |||
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 @@ | |||
1 | /* | ||
2 | * We have copied the definitions in f2c.h required | ||
3 | * to compile clapack.h, modified to support both | ||
4 | * 32 and 64 bit | ||
5 | |||
6 | http://opengrok.creo.hu/dragonfly/xref/src/contrib/gcc-3.4/libf2c/readme.netlib | ||
7 | http://www.ibm.com/developerworks/library/l-port64.html | ||
8 | */ | ||
9 | |||
10 | #ifdef _LP64 | ||
11 | typedef int integer; | ||
12 | typedef unsigned int uinteger; | ||
13 | typedef int logical; | ||
14 | typedef long longint; /* system-dependent */ | ||
15 | typedef unsigned long ulongint; /* system-dependent */ | ||
16 | #else | ||
17 | typedef long int integer; | ||
18 | typedef unsigned long int uinteger; | ||
19 | typedef long int logical; | ||
20 | typedef long long longint; /* system-dependent */ | ||
21 | typedef unsigned long long ulongint; /* system-dependent */ | ||
22 | #endif | ||
23 | |||
24 | typedef char *address; | ||
25 | typedef short int shortint; | ||
26 | typedef float real; | ||
27 | typedef double doublereal; | ||
28 | typedef struct { real r, i; } complex; | ||
29 | typedef struct { doublereal r, i; } doublecomplex; | ||
30 | typedef short int shortlogical; | ||
31 | typedef char logical1; | ||
32 | typedef char integer1; | ||
33 | |||
34 | typedef logical (*L_fp)(); | ||
35 | typedef short ftnlen; | ||
36 | |||
37 | /********************************************************/ | ||
38 | |||
39 | #define FVEC(A) int A##n, float*A##p | ||
40 | #define DVEC(A) int A##n, double*A##p | ||
41 | #define QVEC(A) int A##n, complex*A##p | ||
42 | #define CVEC(A) int A##n, doublecomplex*A##p | ||
43 | #define PVEC(A) int A##n, void* A##p, int A##s | ||
44 | #define FMAT(A) int A##r, int A##c, float* A##p | ||
45 | #define DMAT(A) int A##r, int A##c, double* A##p | ||
46 | #define QMAT(A) int A##r, int A##c, complex* A##p | ||
47 | #define CMAT(A) int A##r, int A##c, doublecomplex* A##p | ||
48 | #define PMAT(A) int A##r, int A##c, void* A##p, int A##s | ||
49 | |||
50 | #define KFVEC(A) int A##n, const float*A##p | ||
51 | #define KDVEC(A) int A##n, const double*A##p | ||
52 | #define KQVEC(A) int A##n, const complex*A##p | ||
53 | #define KCVEC(A) int A##n, const doublecomplex*A##p | ||
54 | #define KPVEC(A) int A##n, const void* A##p, int A##s | ||
55 | #define KFMAT(A) int A##r, int A##c, const float* A##p | ||
56 | #define KDMAT(A) int A##r, int A##c, const double* A##p | ||
57 | #define KQMAT(A) int A##r, int A##c, const complex* A##p | ||
58 | #define KCMAT(A) int A##r, int A##c, const doublecomplex* A##p | ||
59 | #define KPMAT(A) int A##r, int A##c, const void* A##p, int A##s | ||
60 | |||