summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c49
-rw-r--r--packages/base/src/Internal/LAPACK.hs47
2 files changed, 41 insertions, 55 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 80e5720..1018126 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -114,12 +114,12 @@ for(k=0; k<(M##r * M##c); k++) { \
114//////////////////////////////////////////////////////////////////////////////// 114////////////////////////////////////////////////////////////////////////////////
115//////////////////// real svd /////////////////////////////////////////////////// 115//////////////////// real svd ///////////////////////////////////////////////////
116 116
117/* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n, 117int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
118 doublereal *a, integer *lda, doublereal *s, doublereal *u, integer * 118 doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
119 ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, 119 ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
120 integer *info); 120 integer *info);
121 121
122int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { 122int svd_l_R(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
123 integer m = ar; 123 integer m = ar;
124 integer n = ac; 124 integer n = ac;
125 integer q = MIN(m,n); 125 integer q = MIN(m,n);
@@ -145,15 +145,12 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
145 } 145 }
146 } 146 }
147 DEBUGMSG("svd_l_R"); 147 DEBUGMSG("svd_l_R");
148 double *B = (double*)malloc(m*n*sizeof(double));
149 CHECK(!B,MEM);
150 memcpy(B,ap,m*n*sizeof(double));
151 integer lwork = -1; 148 integer lwork = -1;
152 integer res; 149 integer res;
153 // ask for optimal lwork 150 // ask for optimal lwork
154 double ans; 151 double ans;
155 dgesvd_ (jobu,jobvt, 152 dgesvd_ (jobu,jobvt,
156 &m,&n,B,&m, 153 &m,&n,ap,&m,
157 sp, 154 sp,
158 up,&m, 155 up,&m,
159 vp,&ldvt, 156 vp,&ldvt,
@@ -163,7 +160,7 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
163 double * work = (double*)malloc(lwork*sizeof(double)); 160 double * work = (double*)malloc(lwork*sizeof(double));
164 CHECK(!work,MEM); 161 CHECK(!work,MEM);
165 dgesvd_ (jobu,jobvt, 162 dgesvd_ (jobu,jobvt,
166 &m,&n,B,&m, 163 &m,&n,ap,&m,
167 sp, 164 sp,
168 up,&m, 165 up,&m,
169 vp,&ldvt, 166 vp,&ldvt,
@@ -171,18 +168,17 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
171 &res); 168 &res);
172 CHECK(res,res); 169 CHECK(res,res);
173 free(work); 170 free(work);
174 free(B);
175 OK 171 OK
176} 172}
177 173
178// (alternative version) 174// (alternative version)
179 175
180/* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal * 176int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
181 a, integer *lda, doublereal *s, doublereal *u, integer *ldu, 177 a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
182 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork, 178 doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
183 integer *iwork, integer *info); 179 integer *iwork, integer *info);
184 180
185int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { 181int svd_l_Rdd(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
186 integer m = ar; 182 integer m = ar;
187 integer n = ac; 183 integer n = ac;
188 integer q = MIN(m,n); 184 integer q = MIN(m,n);
@@ -202,37 +198,31 @@ int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
202 } 198 }
203 } 199 }
204 DEBUGMSG("svd_l_Rdd"); 200 DEBUGMSG("svd_l_Rdd");
205 double *B = (double*)malloc(m*n*sizeof(double));
206 CHECK(!B,MEM);
207 memcpy(B,ap,m*n*sizeof(double));
208 integer* iwk = (integer*) malloc(8*q*sizeof(integer)); 201 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
209 CHECK(!iwk,MEM); 202 CHECK(!iwk,MEM);
210 integer lwk = -1; 203 integer lwk = -1;
211 integer res; 204 integer res;
212 // ask for optimal lwk 205 // ask for optimal lwk
213 double ans; 206 double ans;
214 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res); 207 dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
215 lwk = ans; 208 lwk = ans;
216 double * workv = (double*)malloc(lwk*sizeof(double)); 209 double * workv = (double*)malloc(lwk*sizeof(double));
217 CHECK(!workv,MEM); 210 CHECK(!workv,MEM);
218 dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res); 211 dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
219 CHECK(res,res); 212 CHECK(res,res);
220 free(iwk); 213 free(iwk);
221 free(workv); 214 free(workv);
222 free(B);
223 OK 215 OK
224} 216}
225 217
226//////////////////// complex svd //////////////////////////////////// 218//////////////////// complex svd ////////////////////////////////////
227 219
228// not in clapack.h
229
230int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, 220int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
231 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u, 221 doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
232 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, 222 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
233 integer *lwork, doublereal *rwork, integer *info); 223 integer *lwork, doublereal *rwork, integer *info);
234 224
235int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { 225int svd_l_C(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
236 integer m = ar; 226 integer m = ar;
237 integer n = ac; 227 integer n = ac;
238 integer q = MIN(m,n); 228 integer q = MIN(m,n);
@@ -257,10 +247,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
257 ldvt = q; 247 ldvt = q;
258 } 248 }
259 }DEBUGMSG("svd_l_C"); 249 }DEBUGMSG("svd_l_C");
260 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); 250
261 CHECK(!B,MEM);
262 memcpy(B,ap,m*n*sizeof(doublecomplex));
263
264 double *rwork = (double*) malloc(5*q*sizeof(double)); 251 double *rwork = (double*) malloc(5*q*sizeof(double));
265 CHECK(!rwork,MEM); 252 CHECK(!rwork,MEM);
266 integer lwork = -1; 253 integer lwork = -1;
@@ -268,7 +255,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
268 // ask for optimal lwork 255 // ask for optimal lwork
269 doublecomplex ans; 256 doublecomplex ans;
270 zgesvd_ (jobu,jobvt, 257 zgesvd_ (jobu,jobvt,
271 &m,&n,B,&m, 258 &m,&n,ap,&m,
272 sp, 259 sp,
273 up,&m, 260 up,&m,
274 vp,&ldvt, 261 vp,&ldvt,
@@ -279,7 +266,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
279 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); 266 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
280 CHECK(!work,MEM); 267 CHECK(!work,MEM);
281 zgesvd_ (jobu,jobvt, 268 zgesvd_ (jobu,jobvt,
282 &m,&n,B,&m, 269 &m,&n,ap,&m,
283 sp, 270 sp,
284 up,&m, 271 up,&m,
285 vp,&ldvt, 272 vp,&ldvt,
@@ -289,7 +276,6 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
289 CHECK(res,res); 276 CHECK(res,res);
290 free(work); 277 free(work);
291 free(rwork); 278 free(rwork);
292 free(B);
293 OK 279 OK
294} 280}
295 281
@@ -298,8 +284,7 @@ int zgesdd_ (char *jobz, integer *m, integer *n,
298 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work, 284 integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
299 integer *lwork, doublereal *rwork, integer* iwork, integer *info); 285 integer *lwork, doublereal *rwork, integer* iwork, integer *info);
300 286
301int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { 287int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
302 //printf("entro\n");
303 integer m = ar; 288 integer m = ar;
304 integer n = ac; 289 integer n = ac;
305 integer q = MIN(m,n); 290 integer q = MIN(m,n);
@@ -319,9 +304,6 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
319 } 304 }
320 } 305 }
321 DEBUGMSG("svd_l_Cdd"); 306 DEBUGMSG("svd_l_Cdd");
322 doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
323 CHECK(!B,MEM);
324 memcpy(B,ap,m*n*sizeof(doublecomplex));
325 integer* iwk = (integer*) malloc(8*q*sizeof(integer)); 307 integer* iwk = (integer*) malloc(8*q*sizeof(integer));
326 CHECK(!iwk,MEM); 308 CHECK(!iwk,MEM);
327 int lrwk; 309 int lrwk;
@@ -337,18 +319,17 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
337 integer res; 319 integer res;
338 // ask for optimal lwk 320 // ask for optimal lwk
339 doublecomplex ans; 321 doublecomplex ans;
340 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res); 322 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
341 lwk = ans.r; 323 lwk = ans.r;
342 //printf("lwk = %ld\n",lwk); 324 //printf("lwk = %ld\n",lwk);
343 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex)); 325 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
344 CHECK(!workv,MEM); 326 CHECK(!workv,MEM);
345 zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res); 327 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
346 //printf("res = %ld\n",res); 328 //printf("res = %ld\n",res);
347 CHECK(res,res); 329 CHECK(res,res);
348 free(workv); // printf("freed workv\n"); 330 free(workv); // printf("freed workv\n");
349 free(rwk); // printf("freed rwk\n"); 331 free(rwk); // printf("freed rwk\n");
350 free(iwk); // printf("freed iwk\n"); 332 free(iwk); // printf("freed iwk\n");
351 free(B); // printf("freed B, salgo\n");
352 OK 333 OK
353} 334}
354 335
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 124e353..c91cddd 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -102,25 +102,26 @@ foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TSVD C
102 102
103-- | Full SVD of a real matrix using LAPACK's /dgesvd/. 103-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
104svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) 104svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
105svdR = svdAux dgesvd "svdR" . fmat 105svdR = svdAux dgesvd "svdR"
106 106
107-- | Full SVD of a real matrix using LAPACK's /dgesdd/. 107-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
108svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) 108svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
109svdRd = svdAux dgesdd "svdRdd" . fmat 109svdRd = svdAux dgesdd "svdRdd"
110 110
111-- | Full SVD of a complex matrix using LAPACK's /zgesvd/. 111-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
112svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) 112svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
113svdC = svdAux zgesvd "svdC" . fmat 113svdC = svdAux zgesvd "svdC"
114 114
115-- | Full SVD of a complex matrix using LAPACK's /zgesdd/. 115-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
116svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) 116svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
117svdCd = svdAux zgesdd "svdCdd" . fmat 117svdCd = svdAux zgesdd "svdCdd"
118 118
119svdAux f st x = unsafePerformIO $ do 119svdAux f st x = unsafePerformIO $ do
120 a <- copy ColumnMajor x
120 u <- createMatrix ColumnMajor r r 121 u <- createMatrix ColumnMajor r r
121 s <- createVector (min r c) 122 s <- createVector (min r c)
122 v <- createMatrix ColumnMajor c c 123 v <- createMatrix ColumnMajor c c
123 f # x # u # s # v #| st 124 f # a # u # s # v #| st
124 return (u,s,v) 125 return (u,s,v)
125 where 126 where
126 r = rows x 127 r = rows x
@@ -129,25 +130,26 @@ svdAux f st x = unsafePerformIO $ do
129 130
130-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'. 131-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
131thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) 132thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
132thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat 133thinSVDR = thinSVDAux dgesvd "thinSVDR"
133 134
134-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'. 135-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
135thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) 136thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
136thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat 137thinSVDC = thinSVDAux zgesvd "thinSVDC"
137 138
138-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'. 139-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
139thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double) 140thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
140thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat 141thinSVDRd = thinSVDAux dgesdd "thinSVDRdd"
141 142
142-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'. 143-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
143thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double)) 144thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
144thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat 145thinSVDCd = thinSVDAux zgesdd "thinSVDCdd"
145 146
146thinSVDAux f st x = unsafePerformIO $ do 147thinSVDAux f st x = unsafePerformIO $ do
148 a <- copy ColumnMajor x
147 u <- createMatrix ColumnMajor r q 149 u <- createMatrix ColumnMajor r q
148 s <- createVector q 150 s <- createVector q
149 v <- createMatrix ColumnMajor q c 151 v <- createMatrix ColumnMajor q c
150 f # x # u # s # v #| st 152 f # a # u # s # v #| st
151 return (u,s,v) 153 return (u,s,v)
152 where 154 where
153 r = rows x 155 r = rows x
@@ -157,23 +159,24 @@ thinSVDAux f st x = unsafePerformIO $ do
157 159
158-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'. 160-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
159svR :: Matrix Double -> Vector Double 161svR :: Matrix Double -> Vector Double
160svR = svAux dgesvd "svR" . fmat 162svR = svAux dgesvd "svR"
161 163
162-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'. 164-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
163svC :: Matrix (Complex Double) -> Vector Double 165svC :: Matrix (Complex Double) -> Vector Double
164svC = svAux zgesvd "svC" . fmat 166svC = svAux zgesvd "svC"
165 167
166-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'. 168-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
167svRd :: Matrix Double -> Vector Double 169svRd :: Matrix Double -> Vector Double
168svRd = svAux dgesdd "svRd" . fmat 170svRd = svAux dgesdd "svRd"
169 171
170-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'. 172-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
171svCd :: Matrix (Complex Double) -> Vector Double 173svCd :: Matrix (Complex Double) -> Vector Double
172svCd = svAux zgesdd "svCd" . fmat 174svCd = svAux zgesdd "svCd"
173 175
174svAux f st x = unsafePerformIO $ do 176svAux f st x = unsafePerformIO $ do
177 a <- copy ColumnMajor x
175 s <- createVector q 178 s <- createVector q
176 g # x # s #| st 179 g # a # s #| st
177 return s 180 return s
178 where 181 where
179 r = rows x 182 r = rows x
@@ -184,16 +187,17 @@ svAux f st x = unsafePerformIO $ do
184 187
185-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'. 188-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
186rightSVR :: Matrix Double -> (Vector Double, Matrix Double) 189rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
187rightSVR = rightSVAux dgesvd "rightSVR" . fmat 190rightSVR = rightSVAux dgesvd "rightSVR"
188 191
189-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'. 192-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
190rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) 193rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
191rightSVC = rightSVAux zgesvd "rightSVC" . fmat 194rightSVC = rightSVAux zgesvd "rightSVC"
192 195
193rightSVAux f st x = unsafePerformIO $ do 196rightSVAux f st x = unsafePerformIO $ do
197 a <- copy ColumnMajor x
194 s <- createVector q 198 s <- createVector q
195 v <- createMatrix ColumnMajor c c 199 v <- createMatrix ColumnMajor c c
196 g # x # s # v #| st 200 g # a # s # v #| st
197 return (s,v) 201 return (s,v)
198 where 202 where
199 r = rows x 203 r = rows x
@@ -204,16 +208,17 @@ rightSVAux f st x = unsafePerformIO $ do
204 208
205-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'. 209-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'A\' and jobvt == \'N\'.
206leftSVR :: Matrix Double -> (Matrix Double, Vector Double) 210leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
207leftSVR = leftSVAux dgesvd "leftSVR" . fmat 211leftSVR = leftSVAux dgesvd "leftSVR"
208 212
209-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'. 213-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
210leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double) 214leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
211leftSVC = leftSVAux zgesvd "leftSVC" . fmat 215leftSVC = leftSVAux zgesvd "leftSVC"
212 216
213leftSVAux f st x = unsafePerformIO $ do 217leftSVAux f st x = unsafePerformIO $ do
218 a <- copy ColumnMajor x
214 u <- createMatrix ColumnMajor r r 219 u <- createMatrix ColumnMajor r r
215 s <- createVector q 220 s <- createVector q
216 g # x # u # s #| st 221 g # a # u # s #| st
217 return (u,s) 222 return (u,s)
218 where 223 where
219 r = rows x 224 r = rows x