diff options
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r-- | packages/base/src/Internal/C/lapack-aux.c | 49 |
1 files changed, 15 insertions, 34 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, | 117 | int 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 | ||
122 | int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { | 122 | int 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 * | 176 | int 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 | ||
185 | int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) { | 181 | int 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 | |||
230 | int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n, | 220 | int 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 | ||
235 | int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { | 225 | int 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 | ||
301 | int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { | 287 | int 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 | ||