diff options
Diffstat (limited to 'lib/Numeric/GSL/gsl-aux.c')
-rw-r--r-- | lib/Numeric/GSL/gsl-aux.c | 286 |
1 files changed, 0 insertions, 286 deletions
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c index bd0a6bd..052cafd 100644 --- a/lib/Numeric/GSL/gsl-aux.c +++ b/lib/Numeric/GSL/gsl-aux.c | |||
@@ -1,11 +1,8 @@ | |||
1 | #include "gsl-aux.h" | 1 | #include "gsl-aux.h" |
2 | #include <gsl/gsl_blas.h> | 2 | #include <gsl/gsl_blas.h> |
3 | #include <gsl/gsl_linalg.h> | ||
4 | #include <gsl/gsl_matrix.h> | ||
5 | #include <gsl/gsl_math.h> | 3 | #include <gsl/gsl_math.h> |
6 | #include <gsl/gsl_errno.h> | 4 | #include <gsl/gsl_errno.h> |
7 | #include <gsl/gsl_fft_complex.h> | 5 | #include <gsl/gsl_fft_complex.h> |
8 | #include <gsl/gsl_eigen.h> | ||
9 | #include <gsl/gsl_integration.h> | 6 | #include <gsl/gsl_integration.h> |
10 | #include <gsl/gsl_deriv.h> | 7 | #include <gsl/gsl_deriv.h> |
11 | #include <gsl/gsl_poly.h> | 8 | #include <gsl/gsl_poly.h> |
@@ -161,47 +158,6 @@ int mapC(int code, KCVEC(x), CVEC(r)) { | |||
161 | } | 158 | } |
162 | 159 | ||
163 | 160 | ||
164 | /* | ||
165 | int scaleR(double* alpha, KRVEC(x), RVEC(r)) { | ||
166 | REQUIRES(xn == rn,BAD_SIZE); | ||
167 | DEBUGMSG("scaleR"); | ||
168 | KDVVIEW(x); | ||
169 | DVVIEW(r); | ||
170 | CHECK( gsl_vector_memcpy(V(r),V(x)) , MEM); | ||
171 | int res = gsl_vector_scale(V(r),*alpha); | ||
172 | CHECK(res,res); | ||
173 | OK | ||
174 | } | ||
175 | |||
176 | int scaleC(gsl_complex *alpha, KCVEC(x), CVEC(r)) { | ||
177 | REQUIRES(xn == rn,BAD_SIZE); | ||
178 | DEBUGMSG("scaleC"); | ||
179 | //KCVVIEW(x); | ||
180 | CVVIEW(r); | ||
181 | gsl_vector_const_view vrx = gsl_vector_const_view_array((double*)xp,xn*2); | ||
182 | gsl_vector_view vrr = gsl_vector_view_array((double*)rp,rn*2); | ||
183 | CHECK(gsl_vector_memcpy(V(vrr),V(vrx)) , MEM); | ||
184 | gsl_blas_zscal(*alpha,V(r)); // void ! | ||
185 | int res = 0; | ||
186 | CHECK(res,res); | ||
187 | OK | ||
188 | } | ||
189 | |||
190 | int addConstantR(double offs, KRVEC(x), RVEC(r)) { | ||
191 | REQUIRES(xn == rn,BAD_SIZE); | ||
192 | DEBUGMSG("addConstantR"); | ||
193 | KDVVIEW(x); | ||
194 | DVVIEW(r); | ||
195 | CHECK(gsl_vector_memcpy(V(r),V(x)), MEM); | ||
196 | int res = gsl_vector_add_constant(V(r),offs); | ||
197 | CHECK(res,res); | ||
198 | OK | ||
199 | } | ||
200 | |||
201 | */ | ||
202 | |||
203 | |||
204 | |||
205 | int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { | 161 | int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { |
206 | int k; | 162 | int k; |
207 | double val = *pval; | 163 | double val = *pval; |
@@ -291,248 +247,6 @@ int zipC(int code, KCVEC(a), KCVEC(b), CVEC(r)) { | |||
291 | 247 | ||
292 | 248 | ||
293 | 249 | ||
294 | |||
295 | int luSolveR(KRMAT(a),KRMAT(b),RMAT(r)) { | ||
296 | REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
297 | DEBUGMSG("luSolveR"); | ||
298 | KDMVIEW(a); | ||
299 | KDMVIEW(b); | ||
300 | DMVIEW(r); | ||
301 | int res; | ||
302 | gsl_matrix *LU = gsl_matrix_alloc(ar,ar); | ||
303 | CHECK(!LU,MEM); | ||
304 | int s; | ||
305 | gsl_permutation * p = gsl_permutation_alloc (ar); | ||
306 | CHECK(!p,MEM); | ||
307 | CHECK(gsl_matrix_memcpy(LU,M(a)),MEM); | ||
308 | res = gsl_linalg_LU_decomp(LU, p, &s); | ||
309 | CHECK(res,res); | ||
310 | int c; | ||
311 | |||
312 | for (c=0; c<bc; c++) { | ||
313 | gsl_vector_const_view colb = gsl_matrix_const_column (M(b), c); | ||
314 | gsl_vector_view colr = gsl_matrix_column (M(r), c); | ||
315 | res = gsl_linalg_LU_solve (LU, p, V(colb), V(colr)); | ||
316 | CHECK(res,res); | ||
317 | } | ||
318 | gsl_permutation_free(p); | ||
319 | gsl_matrix_free(LU); | ||
320 | OK | ||
321 | } | ||
322 | |||
323 | |||
324 | int luSolveC(KCMAT(a),KCMAT(b),CMAT(r)) { | ||
325 | REQUIRES(ar==ac && ac==br && ar==rr && bc==rc,BAD_SIZE); | ||
326 | DEBUGMSG("luSolveC"); | ||
327 | KCMVIEW(a); | ||
328 | KCMVIEW(b); | ||
329 | CMVIEW(r); | ||
330 | gsl_matrix_complex *LU = gsl_matrix_complex_alloc(ar,ar); | ||
331 | CHECK(!LU,MEM); | ||
332 | int s; | ||
333 | gsl_permutation * p = gsl_permutation_alloc (ar); | ||
334 | CHECK(!p,MEM); | ||
335 | CHECK(gsl_matrix_complex_memcpy(LU,M(a)),MEM); | ||
336 | int res; | ||
337 | res = gsl_linalg_complex_LU_decomp(LU, p, &s); | ||
338 | CHECK(res,res); | ||
339 | int c; | ||
340 | for (c=0; c<bc; c++) { | ||
341 | gsl_vector_complex_const_view colb = gsl_matrix_complex_const_column (M(b), c); | ||
342 | gsl_vector_complex_view colr = gsl_matrix_complex_column (M(r), c); | ||
343 | res = gsl_linalg_complex_LU_solve (LU, p, V(colb), V(colr)); | ||
344 | CHECK(res,res); | ||
345 | } | ||
346 | gsl_permutation_free(p); | ||
347 | gsl_matrix_complex_free(LU); | ||
348 | OK | ||
349 | } | ||
350 | |||
351 | |||
352 | int luRaux(KRMAT(a),RVEC(b)) { | ||
353 | REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE); | ||
354 | DEBUGMSG("luRaux"); | ||
355 | KDMVIEW(a); | ||
356 | //DVVIEW(b); | ||
357 | gsl_matrix_view LU = gsl_matrix_view_array(bp,ar,ac); | ||
358 | int s; | ||
359 | gsl_permutation * p = gsl_permutation_alloc (ar); | ||
360 | CHECK(!p,MEM); | ||
361 | CHECK(gsl_matrix_memcpy(M(LU),M(a)),MEM); | ||
362 | gsl_linalg_LU_decomp(M(LU), p, &s); | ||
363 | bp[ar*ar] = s; | ||
364 | int k; | ||
365 | for (k=0; k<ar; k++) { | ||
366 | bp[ar*ar+k+1] = gsl_permutation_get(p,k); | ||
367 | } | ||
368 | gsl_permutation_free(p); | ||
369 | OK | ||
370 | } | ||
371 | |||
372 | int luCaux(KCMAT(a),CVEC(b)) { | ||
373 | REQUIRES(ar==ac && bn==ar*ar+ar+1,BAD_SIZE); | ||
374 | DEBUGMSG("luCaux"); | ||
375 | KCMVIEW(a); | ||
376 | //DVVIEW(b); | ||
377 | gsl_matrix_complex_view LU = gsl_matrix_complex_view_array((double*)bp,ar,ac); | ||
378 | int s; | ||
379 | gsl_permutation * p = gsl_permutation_alloc (ar); | ||
380 | CHECK(!p,MEM); | ||
381 | CHECK(gsl_matrix_complex_memcpy(M(LU),M(a)),MEM); | ||
382 | int res; | ||
383 | res = gsl_linalg_complex_LU_decomp(M(LU), p, &s); | ||
384 | CHECK(res,res); | ||
385 | ((double*)bp)[2*ar*ar] = s; | ||
386 | ((double*)bp)[2*ar*ar+1] = 0; | ||
387 | int k; | ||
388 | for (k=0; k<ar; k++) { | ||
389 | ((double*)bp)[2*ar*ar+2*k+2] = gsl_permutation_get(p,k); | ||
390 | ((double*)bp)[2*ar*ar+2*k+2+1] = 0; | ||
391 | } | ||
392 | gsl_permutation_free(p); | ||
393 | OK | ||
394 | } | ||
395 | |||
396 | int svd(KRMAT(a),RMAT(u), RVEC(s),RMAT(v)) { | ||
397 | REQUIRES(ar==ur && ac==uc && ac==sn && ac==vr && ac==vc,BAD_SIZE); | ||
398 | DEBUGMSG("svd"); | ||
399 | KDMVIEW(a); | ||
400 | DMVIEW(u); | ||
401 | DVVIEW(s); | ||
402 | DMVIEW(v); | ||
403 | gsl_vector *workv = gsl_vector_alloc(ac); | ||
404 | CHECK(!workv,MEM); | ||
405 | gsl_matrix *workm = gsl_matrix_alloc(ac,ac); | ||
406 | CHECK(!workm,MEM); | ||
407 | CHECK(gsl_matrix_memcpy(M(u),M(a)),MEM); | ||
408 | // int res = gsl_linalg_SV_decomp_jacobi (&U.matrix, &V.matrix, &S.vector); | ||
409 | // doesn't work | ||
410 | //int res = gsl_linalg_SV_decomp (&U.matrix, &V.matrix, &S.vector, workv); | ||
411 | int res = gsl_linalg_SV_decomp_mod (M(u), workm, M(v), V(s), workv); | ||
412 | CHECK(res,res); | ||
413 | //gsl_matrix_transpose(M(v)); | ||
414 | gsl_vector_free(workv); | ||
415 | gsl_matrix_free(workm); | ||
416 | OK | ||
417 | } | ||
418 | |||
419 | |||
420 | // for real symmetric matrices | ||
421 | int eigensystemR(KRMAT(x),RVEC(l),RMAT(v)) { | ||
422 | REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE); | ||
423 | DEBUGMSG("eigensystemR (gsl_eigen_symmv)"); | ||
424 | KDMVIEW(x); | ||
425 | DVVIEW(l); | ||
426 | DMVIEW(v); | ||
427 | gsl_matrix *XC = gsl_matrix_alloc(xr,xr); | ||
428 | gsl_matrix_memcpy(XC,M(x)); // needed because the argument is destroyed | ||
429 | // many thanks to Nico Mahlo for the bug report | ||
430 | gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (xc); | ||
431 | int res = gsl_eigen_symmv (XC, V(l), M(v), w); | ||
432 | CHECK(res,res); | ||
433 | gsl_eigen_symmv_free (w); | ||
434 | gsl_matrix_free(XC); | ||
435 | gsl_eigen_symmv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC); | ||
436 | OK | ||
437 | } | ||
438 | |||
439 | // for hermitian matrices | ||
440 | int eigensystemC(KCMAT(x),RVEC(l),CMAT(v)) { | ||
441 | REQUIRES(xr==xc && xr==ln && xr==vr && vr==vc,BAD_SIZE); | ||
442 | DEBUGMSG("eigensystemC"); | ||
443 | KCMVIEW(x); | ||
444 | DVVIEW(l); | ||
445 | CMVIEW(v); | ||
446 | gsl_matrix_complex *XC = gsl_matrix_complex_alloc(xr,xr); | ||
447 | gsl_matrix_complex_memcpy(XC,M(x)); // again needed because the argument is destroyed | ||
448 | gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc (xc); | ||
449 | int res = gsl_eigen_hermv (XC, V(l), M(v), w); | ||
450 | CHECK(res,res); | ||
451 | gsl_eigen_hermv_free (w); | ||
452 | gsl_matrix_complex_free(XC); | ||
453 | gsl_eigen_hermv_sort (V(l), M(v), GSL_EIGEN_SORT_ABS_DESC); | ||
454 | OK | ||
455 | } | ||
456 | |||
457 | int QR(KRMAT(x),RMAT(q),RMAT(r)) { | ||
458 | REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE); | ||
459 | DEBUGMSG("QR"); | ||
460 | KDMVIEW(x); | ||
461 | DMVIEW(q); | ||
462 | DMVIEW(r); | ||
463 | gsl_matrix * a = gsl_matrix_alloc(xr,xc); | ||
464 | gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc)); | ||
465 | gsl_matrix_memcpy(a,M(x)); | ||
466 | int res = gsl_linalg_QR_decomp(a,tau); | ||
467 | CHECK(res,res); | ||
468 | gsl_linalg_QR_unpack(a,tau,M(q),M(r)); | ||
469 | gsl_vector_free(tau); | ||
470 | gsl_matrix_free(a); | ||
471 | OK | ||
472 | } | ||
473 | |||
474 | int QRpacked(KRMAT(x),RMAT(qr),RVEC(tau)) { | ||
475 | //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE); | ||
476 | DEBUGMSG("QRpacked"); | ||
477 | KDMVIEW(x); | ||
478 | DMVIEW(qr); | ||
479 | DVVIEW(tau); | ||
480 | //gsl_matrix * a = gsl_matrix_alloc(xr,xc); | ||
481 | //gsl_vector * tau = gsl_vector_alloc(MIN(xr,xc)); | ||
482 | gsl_matrix_memcpy(M(qr),M(x)); | ||
483 | int res = gsl_linalg_QR_decomp(M(qr),V(tau)); | ||
484 | CHECK(res,res); | ||
485 | OK | ||
486 | } | ||
487 | |||
488 | |||
489 | int QRunpack(KRMAT(xqr),KRVEC(tau),RMAT(q),RMAT(r)) { | ||
490 | //REQUIRES(xr==rr && xc==rc && qr==qc && xr==qr,BAD_SIZE); | ||
491 | DEBUGMSG("QRunpack"); | ||
492 | KDMVIEW(xqr); | ||
493 | KDVVIEW(tau); | ||
494 | DMVIEW(q); | ||
495 | DMVIEW(r); | ||
496 | gsl_linalg_QR_unpack(M(xqr),V(tau),M(q),M(r)); | ||
497 | OK | ||
498 | } | ||
499 | |||
500 | |||
501 | int cholR(KRMAT(x),RMAT(l)) { | ||
502 | REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE); | ||
503 | DEBUGMSG("cholR"); | ||
504 | KDMVIEW(x); | ||
505 | DMVIEW(l); | ||
506 | gsl_matrix_memcpy(M(l),M(x)); | ||
507 | int res = gsl_linalg_cholesky_decomp(M(l)); | ||
508 | CHECK(res,res); | ||
509 | int r,c; | ||
510 | for (r=0; r<xr-1; r++) { | ||
511 | for(c=r+1; c<xc; c++) { | ||
512 | lp[r*lc+c] = 0.; | ||
513 | } | ||
514 | } | ||
515 | OK | ||
516 | } | ||
517 | |||
518 | int cholC(KCMAT(x),CMAT(l)) { | ||
519 | REQUIRES(xr==xc && lr==xr && lr==lc,BAD_SIZE); | ||
520 | DEBUGMSG("cholC"); | ||
521 | KCMVIEW(x); | ||
522 | CMVIEW(l); | ||
523 | gsl_matrix_complex_memcpy(M(l),M(x)); | ||
524 | int res = 0; // gsl_linalg_complex_cholesky_decomp(M(l)); | ||
525 | CHECK(res,res); | ||
526 | gsl_complex zero = {0.,0.}; | ||
527 | int r,c; | ||
528 | for (r=0; r<xr-1; r++) { | ||
529 | for(c=r+1; c<xc; c++) { | ||
530 | lp[r*lc+c] = zero; | ||
531 | } | ||
532 | } | ||
533 | OK | ||
534 | } | ||
535 | |||
536 | int fft(int code, KCVEC(X), CVEC(R)) { | 250 | int fft(int code, KCVEC(X), CVEC(R)) { |
537 | REQUIRES(Xn == Rn,BAD_SIZE); | 251 | REQUIRES(Xn == Rn,BAD_SIZE); |
538 | DEBUGMSG("fft"); | 252 | DEBUGMSG("fft"); |