summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/C/lapack-aux.c
diff options
context:
space:
mode:
Diffstat (limited to 'packages/base/src/Internal/C/lapack-aux.c')
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c107
1 files changed, 107 insertions, 0 deletions
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 30689bf..177d373 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -1086,6 +1086,113 @@ int luS_l_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
1086 OK 1086 OK
1087} 1087}
1088 1088
1089
1090//////////////////// LDL factorization /////////////////////////
1091
1092int dsytrf_(char *uplo, integer *n, doublereal *a, integer *lda, integer *ipiv,
1093 doublereal *work, integer *lwork, integer *info);
1094
1095int ldl_R(DVEC(ipiv), ODMAT(r)) {
1096 integer n = rr;
1097 REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
1098 DEBUGMSG("ldl_R");
1099 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1100 integer res;
1101 integer lda = rXc;
1102 integer lwork = -1;
1103 doublereal ans;
1104 dsytrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
1105 lwork = ceil(ans);
1106 doublereal* work = (doublereal*)malloc(lwork*sizeof(doublereal));
1107 dsytrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
1108 CHECK(res,res);
1109 int k;
1110 for (k=0; k<n; k++) {
1111 ipivp[k] = auxipiv[k];
1112 }
1113 free(auxipiv);
1114 free(work);
1115 OK
1116}
1117
1118
1119int zhetrf_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *ipiv,
1120 doublecomplex *work, integer *lwork, integer *info);
1121
1122int ldl_C(DVEC(ipiv), OCMAT(r)) {
1123 integer n = rr;
1124 REQUIRES(n>=1 && rc==n && ipivn == n, BAD_SIZE);
1125 DEBUGMSG("ldl_R");
1126 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1127 integer res;
1128 integer lda = rXc;
1129 integer lwork = -1;
1130 doublecomplex ans;
1131 zhetrf_ ("L",&n,rp,&lda,auxipiv,&ans,&lwork,&res);
1132 lwork = ceil(ans.r);
1133 doublecomplex* work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
1134 zhetrf_ ("L",&n,rp,&lda,auxipiv,work,&lwork,&res);
1135 CHECK(res,res);
1136 int k;
1137 for (k=0; k<n; k++) {
1138 ipivp[k] = auxipiv[k];
1139 }
1140 free(auxipiv);
1141 free(work);
1142 OK
1143
1144}
1145
1146//////////////////// LDL solve /////////////////////////
1147
1148int dsytrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda,
1149 integer *ipiv, doublereal *b, integer *ldb, integer *info);
1150
1151int ldl_S_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) {
1152 integer m = ar;
1153 integer n = ac;
1154 integer lda = aXc;
1155 integer mrhs = br;
1156 integer nrhs = bc;
1157
1158 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1159 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1160 int k;
1161 for (k=0; k<n; k++) {
1162 auxipiv[k] = (integer)ipivp[k];
1163 }
1164 integer res;
1165 dsytrs_ ("L",&n,&nrhs,(/*no const (!?)*/ double*)ap,&lda,auxipiv,bp,&mrhs,&res);
1166 CHECK(res,res);
1167 free(auxipiv);
1168 OK
1169}
1170
1171
1172int zhetrs_(char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda,
1173 integer *ipiv, doublecomplex *b, integer *ldb, integer *info);
1174
1175int ldl_S_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) {
1176 integer m = ar;
1177 integer n = ac;
1178 integer lda = aXc;
1179 integer mrhs = br;
1180 integer nrhs = bc;
1181
1182 REQUIRES(m==n && m==mrhs && m==ipivn,BAD_SIZE);
1183 integer* auxipiv = (integer*)malloc(n*sizeof(integer));
1184 int k;
1185 for (k=0; k<n; k++) {
1186 auxipiv[k] = (integer)ipivp[k];
1187 }
1188 integer res;
1189 zhetrs_ ("L",&n,&nrhs,(doublecomplex*)ap,&lda,auxipiv,bp,&mrhs,&res);
1190 CHECK(res,res);
1191 free(auxipiv);
1192 OK
1193}
1194
1195
1089//////////////////// Matrix Product ///////////////////////// 1196//////////////////// Matrix Product /////////////////////////
1090 1197
1091void dgemm_(char *, char *, integer *, integer *, integer *, 1198void dgemm_(char *, char *, integer *, integer *, integer *,