summaryrefslogtreecommitdiff
path: root/lib/Numeric/GSL/gsl-aux.c
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/GSL/gsl-aux.c')
-rw-r--r--lib/Numeric/GSL/gsl-aux.c286
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/*
165int 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
176int 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
190int 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
205int mapValR(int code, double* pval, KRVEC(x), RVEC(r)) { 161int 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
295int 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
324int 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
352int 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
372int 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
396int 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
421int 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
440int 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
457int 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
474int 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
489int 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
501int 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
518int 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
536int fft(int code, KCVEC(X), CVEC(R)) { 250int 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");