summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/C
diff options
context:
space:
mode:
authorDominic Steinitz <dominic@steinitz.org>2017-03-21 17:35:43 +0000
committerDominic Steinitz <dominic@steinitz.org>2017-03-21 17:35:43 +0000
commit49d718705d205d62aea2762445f95735a671f305 (patch)
tree589b5c4396647f48b941d313432647ecb53ef606 /packages/base/src/Internal/C
parentfa1642dcf26f1da0a6f4c1324bcd1e8baf9fd478 (diff)
Add tridiagonal solver and tests for it and triagonal solver.
Diffstat (limited to 'packages/base/src/Internal/C')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c72
1 files changed, 72 insertions, 0 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 4a8129c..5018a45 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -668,6 +668,78 @@ int triSolveC_l_l(KOCMAT(a),OCMAT(b)) {
668 OK 668 OK
669} 669}
670 670
671//////// tridiagonal real linear system ////////////
672
673int dgttrf_(integer *n,
674 doublereal *dl, doublereal *d, doublereal *du, doublereal *du2,
675 integer *ipiv,
676 integer *info);
677
678int dgttrs_(char *trans, integer *n, integer *nrhs,
679 doublereal *dl, doublereal *d, doublereal *du, doublereal *du2,
680 integer *ipiv, doublereal *b, integer *ldb,
681 integer *info);
682
683int triDiagSolveR_l(DVEC(dl), DVEC(d), DVEC(du), ODMAT(b)) {
684 integer n = dn;
685 integer nhrs = bc;
686 REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE);
687 DEBUGMSG("triDiagSolveR_l");
688 integer res;
689 integer* ipiv = (integer*)malloc(n*sizeof(integer));
690 double* du2 = (double*)malloc((n - 2)*sizeof(double));
691 dgttrf_ (&n,
692 dlp, dp, dup, du2,
693 ipiv,
694 &res);
695 CHECK(res,res);
696 dgttrs_ ("N",
697 &n,&nhrs,
698 dlp, dp, dup, du2,
699 ipiv, bp, &n,
700 &res);
701 CHECK(res,res);
702 free(ipiv);
703 free(du2);
704 OK
705}
706
707//////// tridiagonal complex linear system ////////////
708
709int zgttrf_(integer *n,
710 doublecomplex *dl, doublecomplex *d, doublecomplex *du, doublecomplex *du2,
711 integer *ipiv,
712 integer *info);
713
714int zgttrs_(char *trans, integer *n, integer *nrhs,
715 doublecomplex *dl, doublecomplex *d, doublecomplex *du, doublecomplex *du2,
716 integer *ipiv, doublecomplex *b, integer *ldb,
717 integer *info);
718
719int triDiagSolveC_l(CVEC(dl), CVEC(d), CVEC(du), OCMAT(b)) {
720 integer n = dn;
721 integer nhrs = bc;
722 REQUIRES(n >= 1 && dln == dn - 1 && dun == dn - 1 && br == n, BAD_SIZE);
723 DEBUGMSG("triDiagSolveC_l");
724 integer res;
725 integer* ipiv = (integer*)malloc(n*sizeof(integer));
726 doublecomplex* du2 = (doublecomplex*)malloc((n - 2)*sizeof(doublecomplex));
727 zgttrf_ (&n,
728 dlp, dp, dup, du2,
729 ipiv,
730 &res);
731 CHECK(res,res);
732 zgttrs_ ("N",
733 &n,&nhrs,
734 dlp, dp, dup, du2,
735 ipiv, bp, &n,
736 &res);
737 CHECK(res,res);
738 free(ipiv);
739 free(du2);
740 OK
741}
742
671//////////////////// least squares real linear system //////////// 743//////////////////// least squares real linear system ////////////
672 744
673int dgels_(char *trans, integer *m, integer *n, integer * 745int dgels_(char *trans, integer *m, integer *n, integer *