summaryrefslogtreecommitdiff
path: root/packages/base
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c165
-rw-r--r--packages/base/src/Internal/LAPACK.hs74
2 files changed, 98 insertions, 141 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
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index c91cddd..13340f2 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -349,57 +349,75 @@ eigOnlyH = vrev . fst. eigSHAux (zheev 0) "eigH'"
349vrev = flatten . flipud . reshape 1 349vrev = flatten . flipud . reshape 1
350 350
351----------------------------------------------------------------------------- 351-----------------------------------------------------------------------------
352foreign import ccall unsafe "linearSolveR_l" dgesv :: TMMM R 352foreign import ccall unsafe "linearSolveR_l" dgesv :: R ::> R ::> Ok
353foreign import ccall unsafe "linearSolveC_l" zgesv :: TMMM C 353foreign import ccall unsafe "linearSolveC_l" zgesv :: C ::> C ::> Ok
354foreign import ccall unsafe "cholSolveR_l" dpotrs :: TMMM R
355foreign import ccall unsafe "cholSolveC_l" zpotrs :: TMMM C
356 354
357linearSolveSQAux g f st a b 355linearSolveSQAux g f st a b
358 | n1==n2 && n1==r = unsafePerformIO . g $ do 356 | n1==n2 && n1==r = unsafePerformIO . g $ do
359 s <- createMatrix ColumnMajor r c 357 a' <- copy ColumnMajor a
360 f # a # b # s #| st 358 s <- copy ColumnMajor b
359 f # a' # s #| st
361 return s 360 return s
362 | otherwise = error $ st ++ " of nonsquare matrix" 361 | otherwise = error $ st ++ " of nonsquare matrix"
363 where 362 where
364 n1 = rows a 363 n1 = rows a
365 n2 = cols a 364 n2 = cols a
366 r = rows b 365 r = rows b
367 c = cols b
368 366
369-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'. 367-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
370linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double 368linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
371linearSolveR a b = linearSolveSQAux id dgesv "linearSolveR" (fmat a) (fmat b) 369linearSolveR a b = linearSolveSQAux id dgesv "linearSolveR" a b
372 370
373mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double) 371mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
374mbLinearSolveR a b = linearSolveSQAux mbCatch dgesv "linearSolveR" (fmat a) (fmat b) 372mbLinearSolveR a b = linearSolveSQAux mbCatch dgesv "linearSolveR" a b
375 373
376 374
377-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'. 375-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.
378linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 376linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
379linearSolveC a b = linearSolveSQAux id zgesv "linearSolveC" (fmat a) (fmat b) 377linearSolveC a b = linearSolveSQAux id zgesv "linearSolveC" a b
380 378
381mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double)) 379mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
382mbLinearSolveC a b = linearSolveSQAux mbCatch zgesv "linearSolveC" (fmat a) (fmat b) 380mbLinearSolveC a b = linearSolveSQAux mbCatch zgesv "linearSolveC" a b
381
382--------------------------------------------------------------------------------
383foreign import ccall unsafe "cholSolveR_l" dpotrs :: R ::> R ::> Ok
384foreign import ccall unsafe "cholSolveC_l" zpotrs :: C ::> C ::> Ok
385
386
387linearSolveSQAux2 g f st a b
388 | n1==n2 && n1==r = unsafePerformIO . g $ do
389 s <- copy ColumnMajor b
390 f # a # s #| st
391 return s
392 | otherwise = error $ st ++ " of nonsquare matrix"
393 where
394 n1 = rows a
395 n2 = cols a
396 r = rows b
383 397
384-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'. 398-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
385cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double 399cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
386cholSolveR a b = linearSolveSQAux id dpotrs "cholSolveR" (fmat a) (fmat b) 400cholSolveR a b = linearSolveSQAux2 id dpotrs "cholSolveR" (fmat a) b
387 401
388-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'. 402-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.
389cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 403cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
390cholSolveC a b = linearSolveSQAux id zpotrs "cholSolveC" (fmat a) (fmat b) 404cholSolveC a b = linearSolveSQAux2 id zpotrs "cholSolveC" (fmat a) b
391 405
392----------------------------------------------------------------------------------- 406-----------------------------------------------------------------------------------
393 407
394foreign import ccall unsafe "linearSolveLSR_l" dgels :: TMMM R 408foreign import ccall unsafe "linearSolveLSR_l" dgels :: R ::> R ::> Ok
395foreign import ccall unsafe "linearSolveLSC_l" zgels :: TMMM C 409foreign import ccall unsafe "linearSolveLSC_l" zgels :: C ::> C ::> Ok
396foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> TMMM R 410foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> R ::> R ::> Ok
397foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> TMMM C 411foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> C ::> C ::> Ok
398 412
399linearSolveAux f st a b = unsafePerformIO $ do 413linearSolveAux f st a b
400 r <- createMatrix ColumnMajor (max m n) nrhs 414 | m == rows b = unsafePerformIO $ do
401 f # a # b # r #| st 415 a' <- copy ColumnMajor a
402 return r 416 r <- createMatrix ColumnMajor (max m n) nrhs
417 setRect 0 0 b r
418 f # a' # r #| st
419 return r
420 | otherwise = error $ "different number of rows in linearSolve ("++st++")"
403 where 421 where
404 m = rows a 422 m = rows a
405 n = cols a 423 n = cols a
@@ -408,12 +426,12 @@ linearSolveAux f st a b = unsafePerformIO $ do
408-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'. 426-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'.
409linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double 427linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
410linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $ 428linearSolveLSR a b = subMatrix (0,0) (cols a, cols b) $
411 linearSolveAux dgels "linearSolverLSR" (fmat a) (fmat b) 429 linearSolveAux dgels "linearSolverLSR" a b
412 430
413-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'. 431-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'.
414linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) 432linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
415linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $ 433linearSolveLSC a b = subMatrix (0,0) (cols a, cols b) $
416 linearSolveAux zgels "linearSolveLSC" (fmat a) (fmat b) 434 linearSolveAux zgels "linearSolveLSC" a b
417 435
418-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. 436-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
419linearSolveSVDR :: Maybe Double -- ^ rcond 437linearSolveSVDR :: Maybe Double -- ^ rcond
@@ -421,8 +439,8 @@ linearSolveSVDR :: Maybe Double -- ^ rcond
421 -> Matrix Double -- ^ right hand sides (as columns) 439 -> Matrix Double -- ^ right hand sides (as columns)
422 -> Matrix Double -- ^ solution vectors (as columns) 440 -> Matrix Double -- ^ solution vectors (as columns)
423linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ 441linearSolveSVDR (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
424 linearSolveAux (dgelss rcond) "linearSolveSVDR" (fmat a) (fmat b) 442 linearSolveAux (dgelss rcond) "linearSolveSVDR" a b
425linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) (fmat a) (fmat b) 443linearSolveSVDR Nothing a b = linearSolveSVDR (Just (-1)) a b
426 444
427-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used. 445-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
428linearSolveSVDC :: Maybe Double -- ^ rcond 446linearSolveSVDC :: Maybe Double -- ^ rcond
@@ -430,8 +448,8 @@ linearSolveSVDC :: Maybe Double -- ^ rcond
430 -> Matrix (Complex Double) -- ^ right hand sides (as columns) 448 -> Matrix (Complex Double) -- ^ right hand sides (as columns)
431 -> Matrix (Complex Double) -- ^ solution vectors (as columns) 449 -> Matrix (Complex Double) -- ^ solution vectors (as columns)
432linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $ 450linearSolveSVDC (Just rcond) a b = subMatrix (0,0) (cols a, cols b) $
433 linearSolveAux (zgelss rcond) "linearSolveSVDC" (fmat a) (fmat b) 451 linearSolveAux (zgelss rcond) "linearSolveSVDC" a b
434linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) (fmat a) (fmat b) 452linearSolveSVDC Nothing a b = linearSolveSVDC (Just (-1)) a b
435 453
436----------------------------------------------------------------------------------- 454-----------------------------------------------------------------------------------
437 455