summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/C
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c165
1 files changed, 52 insertions, 113 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 1018126..ca60846 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -314,22 +314,19 @@ int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
314 } 314 }
315 double *rwk = (double*)malloc(lrwk*sizeof(double));; 315 double *rwk = (double*)malloc(lrwk*sizeof(double));;
316 CHECK(!rwk,MEM); 316 CHECK(!rwk,MEM);
317 //printf("%s %ld %d\n",jobz,q,lrwk);
318 integer lwk = -1; 317 integer lwk = -1;
319 integer res; 318 integer res;
320 // ask for optimal lwk 319 // ask for optimal lwk
321 doublecomplex ans; 320 doublecomplex ans;
322 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res); 321 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
323 lwk = ans.r; 322 lwk = ans.r;
324 //printf("lwk = %ld\n",lwk);
325 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex)); 323 doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
326 CHECK(!workv,MEM); 324 CHECK(!workv,MEM);
327 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res); 325 zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
328 //printf("res = %ld\n",res);
329 CHECK(res,res); 326 CHECK(res,res);
330 free(workv); // printf("freed workv\n"); 327 free(workv);
331 free(rwk); // printf("freed rwk\n"); 328 free(rwk);
332 free(iwk); // printf("freed iwk\n"); 329 free(iwk);
333 OK 330 OK
334} 331}
335 332
@@ -498,80 +495,72 @@ int eig_l_H(int wantV,DVEC(s),OCMAT(v)) {
498 495
499//////////////////// general real linear system //////////// 496//////////////////// general real linear system ////////////
500 497
501/* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer 498int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
502 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info); 499 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
503 500
504int linearSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { 501int linearSolveR_l(ODMAT(a),ODMAT(b)) {
505 integer n = ar; 502 integer n = ar;
506 integer nhrs = bc; 503 integer nhrs = bc;
507 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); 504 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
508 DEBUGMSG("linearSolveR_l"); 505 DEBUGMSG("linearSolveR_l");
509 double*AC = (double*)malloc(n*n*sizeof(double));
510 memcpy(AC,ap,n*n*sizeof(double));
511 memcpy(xp,bp,n*nhrs*sizeof(double));
512 integer * ipiv = (integer*)malloc(n*sizeof(integer)); 506 integer * ipiv = (integer*)malloc(n*sizeof(integer));
513 integer res; 507 integer res;
514 dgesv_ (&n,&nhrs, 508 dgesv_ (&n,&nhrs,
515 AC, &n, 509 ap, &n,
516 ipiv, 510 ipiv,
517 xp, &n, 511 bp, &n,
518 &res); 512 &res);
519 if(res>0) { 513 if(res>0) {
520 return SINGULAR; 514 return SINGULAR;
521 } 515 }
522 CHECK(res,res); 516 CHECK(res,res);
523 free(ipiv); 517 free(ipiv);
524 free(AC);
525 OK 518 OK
526} 519}
527 520
528//////////////////// general complex linear system //////////// 521//////////////////// general complex linear system ////////////
529 522
530/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a, 523int zgesv_(integer *n, integer *nrhs, doublecomplex *a,
531 integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer * 524 integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *
532 info); 525 info);
533 526
534int linearSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { 527int linearSolveC_l(OCMAT(a),OCMAT(b)) {
535 integer n = ar; 528 integer n = ar;
536 integer nhrs = bc; 529 integer nhrs = bc;
537 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); 530 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
538 DEBUGMSG("linearSolveC_l"); 531 DEBUGMSG("linearSolveC_l");
539 doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
540 memcpy(AC,ap,n*n*sizeof(doublecomplex));
541 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
542 integer * ipiv = (integer*)malloc(n*sizeof(integer)); 532 integer * ipiv = (integer*)malloc(n*sizeof(integer));
543 integer res; 533 integer res;
544 zgesv_ (&n,&nhrs, 534 zgesv_ (&n,&nhrs,
545 AC, &n, 535 ap, &n,
546 ipiv, 536 ipiv,
547 xp, &n, 537 bp, &n,
548 &res); 538 &res);
549 if(res>0) { 539 if(res>0) {
550 return SINGULAR; 540 return SINGULAR;
551 } 541 }
552 CHECK(res,res); 542 CHECK(res,res);
553 free(ipiv); 543 free(ipiv);
554 free(AC);
555 OK 544 OK
556} 545}
557 546
558//////// symmetric positive definite real linear system using Cholesky //////////// 547//////// symmetric positive definite real linear system using Cholesky ////////////
559 548
560/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs, 549int dpotrs_(char *uplo, integer *n, integer *nrhs,
561 doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * 550 doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
562 info); 551 info);
563 552
564int cholSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { 553int cholSolveR_l(KODMAT(a),ODMAT(b)) {
565 integer n = ar; 554 integer n = ar;
555 integer lda = aXc;
566 integer nhrs = bc; 556 integer nhrs = bc;
567 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); 557 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
568 DEBUGMSG("cholSolveR_l"); 558 DEBUGMSG("cholSolveR_l");
569 memcpy(xp,bp,n*nhrs*sizeof(double));
570 integer res; 559 integer res;
571 dpotrs_ ("U", 560 dpotrs_ ("U",
572 &n,&nhrs, 561 &n,&nhrs,
573 (double*)ap, &n, 562 (double*)ap, &lda,
574 xp, &n, 563 bp, &n,
575 &res); 564 &res);
576 CHECK(res,res); 565 CHECK(res,res);
577 OK 566 OK
@@ -579,21 +568,21 @@ int cholSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) {
579 568
580//////// Hermitian positive definite real linear system using Cholesky //////////// 569//////// Hermitian positive definite real linear system using Cholesky ////////////
581 570
582/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, 571int zpotrs_(char *uplo, integer *n, integer *nrhs,
583 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, 572 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
584 integer *info); 573 integer *info);
585 574
586int cholSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { 575int cholSolveC_l(KOCMAT(a),OCMAT(b)) {
587 integer n = ar; 576 integer n = ar;
577 integer lda = aXc;
588 integer nhrs = bc; 578 integer nhrs = bc;
589 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); 579 REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE);
590 DEBUGMSG("cholSolveC_l"); 580 DEBUGMSG("cholSolveC_l");
591 memcpy(xp,bp,n*nhrs*sizeof(doublecomplex));
592 integer res; 581 integer res;
593 zpotrs_ ("U", 582 zpotrs_ ("U",
594 &n,&nhrs, 583 &n,&nhrs,
595 (doublecomplex*)ap, &n, 584 (doublecomplex*)ap, &lda,
596 xp, &n, 585 bp, &n,
597 &res); 586 &res);
598 CHECK(res,res); 587 CHECK(res,res);
599 OK 588 OK
@@ -601,41 +590,30 @@ int cholSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) {
601 590
602//////////////////// least squares real linear system //////////// 591//////////////////// least squares real linear system ////////////
603 592
604/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer * 593int dgels_(char *trans, integer *m, integer *n, integer *
605 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, 594 nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb,
606 doublereal *work, integer *lwork, integer *info); 595 doublereal *work, integer *lwork, integer *info);
607 596
608int linearSolveLSR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { 597int linearSolveLSR_l(ODMAT(a),ODMAT(b)) {
609 integer m = ar; 598 integer m = ar;
610 integer n = ac; 599 integer n = ac;
611 integer nrhs = bc; 600 integer nrhs = bc;
612 integer ldb = xr; 601 integer ldb = bXc;
613 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); 602 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
614 DEBUGMSG("linearSolveLSR_l"); 603 DEBUGMSG("linearSolveLSR_l");
615 double*AC = (double*)malloc(m*n*sizeof(double));
616 memcpy(AC,ap,m*n*sizeof(double));
617 if (m>=n) {
618 memcpy(xp,bp,m*nrhs*sizeof(double));
619 } else {
620 int k;
621 for(k = 0; k<nrhs; k++) {
622 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
623 }
624 }
625 integer res; 604 integer res;
626 integer lwork = -1; 605 integer lwork = -1;
627 double ans; 606 double ans;
628 dgels_ ("N",&m,&n,&nrhs, 607 dgels_ ("N",&m,&n,&nrhs,
629 AC,&m, 608 ap,&m,
630 xp,&ldb, 609 bp,&ldb,
631 &ans,&lwork, 610 &ans,&lwork,
632 &res); 611 &res);
633 lwork = ceil(ans); 612 lwork = ceil(ans);
634 //printf("ans = %d\n",lwork);
635 double * work = (double*)malloc(lwork*sizeof(double)); 613 double * work = (double*)malloc(lwork*sizeof(double));
636 dgels_ ("N",&m,&n,&nrhs, 614 dgels_ ("N",&m,&n,&nrhs,
637 AC,&m, 615 ap,&m,
638 xp,&ldb, 616 bp,&ldb,
639 work,&lwork, 617 work,&lwork,
640 &res); 618 &res);
641 if(res>0) { 619 if(res>0) {
@@ -643,47 +621,35 @@ int linearSolveLSR_l(KODMAT(a),KODMAT(b),ODMAT(x)) {
643 } 621 }
644 CHECK(res,res); 622 CHECK(res,res);
645 free(work); 623 free(work);
646 free(AC);
647 OK 624 OK
648} 625}
649 626
650//////////////////// least squares complex linear system //////////// 627//////////////////// least squares complex linear system ////////////
651 628
652/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer * 629int zgels_(char *trans, integer *m, integer *n, integer *
653 nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, 630 nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
654 doublecomplex *work, integer *lwork, integer *info); 631 doublecomplex *work, integer *lwork, integer *info);
655 632
656int linearSolveLSC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { 633int linearSolveLSC_l(OCMAT(a),OCMAT(b)) {
657 integer m = ar; 634 integer m = ar;
658 integer n = ac; 635 integer n = ac;
659 integer nrhs = bc; 636 integer nrhs = bc;
660 integer ldb = xr; 637 integer ldb = bXc;
661 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); 638 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
662 DEBUGMSG("linearSolveLSC_l"); 639 DEBUGMSG("linearSolveLSC_l");
663 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
664 memcpy(AC,ap,m*n*sizeof(doublecomplex));
665 if (m>=n) {
666 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
667 } else {
668 int k;
669 for(k = 0; k<nrhs; k++) {
670 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
671 }
672 }
673 integer res; 640 integer res;
674 integer lwork = -1; 641 integer lwork = -1;
675 doublecomplex ans; 642 doublecomplex ans;
676 zgels_ ("N",&m,&n,&nrhs, 643 zgels_ ("N",&m,&n,&nrhs,
677 AC,&m, 644 ap,&m,
678 xp,&ldb, 645 bp,&ldb,
679 &ans,&lwork, 646 &ans,&lwork,
680 &res); 647 &res);
681 lwork = ceil(ans.r); 648 lwork = ceil(ans.r);
682 //printf("ans = %d\n",lwork);
683 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); 649 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
684 zgels_ ("N",&m,&n,&nrhs, 650 zgels_ ("N",&m,&n,&nrhs,
685 AC,&m, 651 ap,&m,
686 xp,&ldb, 652 bp,&ldb,
687 work,&lwork, 653 work,&lwork,
688 &res); 654 &res);
689 if(res>0) { 655 if(res>0) {
@@ -691,52 +657,40 @@ int linearSolveLSC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) {
691 } 657 }
692 CHECK(res,res); 658 CHECK(res,res);
693 free(work); 659 free(work);
694 free(AC);
695 OK 660 OK
696} 661}
697 662
698//////////////////// least squares real linear system using SVD //////////// 663//////////////////// least squares real linear system using SVD ////////////
699 664
700/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, 665int dgelss_(integer *m, integer *n, integer *nrhs,
701 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * 666 doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *
702 s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, 667 s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork,
703 integer *info); 668 integer *info);
704 669
705int linearSolveSVDR_l(double rcond,KODMAT(a),KODMAT(b),ODMAT(x)) { 670int linearSolveSVDR_l(double rcond,ODMAT(a),ODMAT(b)) {
706 integer m = ar; 671 integer m = ar;
707 integer n = ac; 672 integer n = ac;
708 integer nrhs = bc; 673 integer nrhs = bc;
709 integer ldb = xr; 674 integer ldb = bXc;
710 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); 675 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
711 DEBUGMSG("linearSolveSVDR_l"); 676 DEBUGMSG("linearSolveSVDR_l");
712 double*AC = (double*)malloc(m*n*sizeof(double));
713 double*S = (double*)malloc(MIN(m,n)*sizeof(double)); 677 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
714 memcpy(AC,ap,m*n*sizeof(double));
715 if (m>=n) {
716 memcpy(xp,bp,m*nrhs*sizeof(double));
717 } else {
718 int k;
719 for(k = 0; k<nrhs; k++) {
720 memcpy(xp+ldb*k,bp+m*k,m*sizeof(double));
721 }
722 }
723 integer res; 678 integer res;
724 integer lwork = -1; 679 integer lwork = -1;
725 integer rank; 680 integer rank;
726 double ans; 681 double ans;
727 dgelss_ (&m,&n,&nrhs, 682 dgelss_ (&m,&n,&nrhs,
728 AC,&m, 683 ap,&m,
729 xp,&ldb, 684 bp,&ldb,
730 S, 685 S,
731 &rcond,&rank, 686 &rcond,&rank,
732 &ans,&lwork, 687 &ans,&lwork,
733 &res); 688 &res);
734 lwork = ceil(ans); 689 lwork = ceil(ans);
735 //printf("ans = %d\n",lwork);
736 double * work = (double*)malloc(lwork*sizeof(double)); 690 double * work = (double*)malloc(lwork*sizeof(double));
737 dgelss_ (&m,&n,&nrhs, 691 dgelss_ (&m,&n,&nrhs,
738 AC,&m, 692 ap,&m,
739 xp,&ldb, 693 bp,&ldb,
740 S, 694 S,
741 &rcond,&rank, 695 &rcond,&rank,
742 work,&lwork, 696 work,&lwork,
@@ -747,57 +701,43 @@ int linearSolveSVDR_l(double rcond,KODMAT(a),KODMAT(b),ODMAT(x)) {
747 CHECK(res,res); 701 CHECK(res,res);
748 free(work); 702 free(work);
749 free(S); 703 free(S);
750 free(AC);
751 OK 704 OK
752} 705}
753 706
754//////////////////// least squares complex linear system using SVD //////////// 707//////////////////// least squares complex linear system using SVD ////////////
755 708
756// not in clapack.h
757
758int zgelss_(integer *m, integer *n, integer *nhrs, 709int zgelss_(integer *m, integer *n, integer *nhrs,
759 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s, 710 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublereal *s,
760 doublereal *rcond, integer* rank, 711 doublereal *rcond, integer* rank,
761 doublecomplex *work, integer* lwork, doublereal* rwork, 712 doublecomplex *work, integer* lwork, doublereal* rwork,
762 integer *info); 713 integer *info);
763 714
764int linearSolveSVDC_l(double rcond, KOCMAT(a),KOCMAT(b),OCMAT(x)) { 715int linearSolveSVDC_l(double rcond, OCMAT(a),OCMAT(b)) {
765 integer m = ar; 716 integer m = ar;
766 integer n = ac; 717 integer n = ac;
767 integer nrhs = bc; 718 integer nrhs = bc;
768 integer ldb = xr; 719 integer ldb = bXc;
769 REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); 720 REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE);
770 DEBUGMSG("linearSolveSVDC_l"); 721 DEBUGMSG("linearSolveSVDC_l");
771 doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
772 double*S = (double*)malloc(MIN(m,n)*sizeof(double)); 722 double*S = (double*)malloc(MIN(m,n)*sizeof(double));
773 double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); 723 double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double));
774 memcpy(AC,ap,m*n*sizeof(doublecomplex));
775 if (m>=n) {
776 memcpy(xp,bp,m*nrhs*sizeof(doublecomplex));
777 } else {
778 int k;
779 for(k = 0; k<nrhs; k++) {
780 memcpy(xp+ldb*k,bp+m*k,m*sizeof(doublecomplex));
781 }
782 }
783 integer res; 724 integer res;
784 integer lwork = -1; 725 integer lwork = -1;
785 integer rank; 726 integer rank;
786 doublecomplex ans; 727 doublecomplex ans;
787 zgelss_ (&m,&n,&nrhs, 728 zgelss_ (&m,&n,&nrhs,
788 AC,&m, 729 ap,&m,
789 xp,&ldb, 730 bp,&ldb,
790 S, 731 S,
791 &rcond,&rank, 732 &rcond,&rank,
792 &ans,&lwork, 733 &ans,&lwork,
793 RWORK, 734 RWORK,
794 &res); 735 &res);
795 lwork = ceil(ans.r); 736 lwork = ceil(ans.r);
796 //printf("ans = %d\n",lwork);
797 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); 737 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
798 zgelss_ (&m,&n,&nrhs, 738 zgelss_ (&m,&n,&nrhs,
799 AC,&m, 739 ap,&m,
800 xp,&ldb, 740 bp,&ldb,
801 S, 741 S,
802 &rcond,&rank, 742 &rcond,&rank,
803 work,&lwork, 743 work,&lwork,
@@ -810,7 +750,6 @@ int linearSolveSVDC_l(double rcond, KOCMAT(a),KOCMAT(b),OCMAT(x)) {
810 free(work); 750 free(work);
811 free(RWORK); 751 free(RWORK);
812 free(S); 752 free(S);
813 free(AC);
814 OK 753 OK
815} 754}
816 755