6#define max(a,b) ((a) >= (b) ? (a) : (b))
7#define dmax(a,b) (double)max(a,b)
9static int lsame2(
const char c1[],
const char c2[]){
10 if (c1[0]==c2[0])
return 1;
14 int dtrsm2(
char *side,
char *uplo,
char *transa,
char *diag,
15 integer *m, integer *n,
double *alpha,
double *a, integer *
16 lda,
double *b, integer *ldb)
19 integer a_dim1, a_offset, b_dim1, b_offset, i1, i2, i3;
22 static double temp,alpha2,aref;
23 static integer i, j, k;
27 static logical nounit;
28#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
29#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
32 a_offset = 1 + a_dim1 * 1;
35 b_offset = 1 + b_dim1 * 1;
38 lside = lsame2(side,
"L");
44 nounit = lsame2(diag,
"N");
45 upper = lsame2(uplo,
"U");
47 if (! lside && ! lsame2(side,
"R")) {
49 }
else if (! upper && ! lsame2(uplo,
"L")) {
51 }
else if (! lsame2(transa,
"N") && ! lsame2(transa,
52 "T") && ! lsame2(transa,
"C")) {
54 }
else if (! lsame2(diag,
"U") && ! lsame2(diag,
61 }
else if (*lda < max(1,nrowa)) {
63 }
else if (*ldb < max(1,*m)) {
76 for (j = 1; j <= i1; ++j) {
78 for (i = 1; i <= i2; ++i) {
88 if (lsame2(transa,
"N")) {
92 for (j = 1; j <= i1; ++j) {
96 for (i = 1; i <= i2; ++i) {
97 b_ref(i, j) *= alpha2;
101 for (k = *m; k >= 1; --k) {
102 if (b_ref(k, j) != 0.) {
104 b_ref(k, j) /= a_ref(k, k);
108 for (i = 1; i <= i2; ++i) {
109 b_ref(i, j) += aref * a_ref(i, k);
119 for (j = 1; j <= i1; ++j) {
123 for (i = 1; i <= i2; ++i) {
124 b_ref(i, j) *= alpha2;
129 for (k = 1; k <= i2; ++k) {
130 if (b_ref(k, j) != 0.) {
132 b_ref(k, j) /= a_ref(k, k);
136 for (i = k + 1; i <= i3; ++i) {
137 b_ref(i, j) += aref * a_ref(i, k);
150 for (j = 1; j <= i1; ++j) {
153 for (i = 1; i <= i2; ++i) {
154 temp = alpha2 * b_ref(i, j);
156 for (k = 1; k <= i3; ++k) {
157 temp -= a_ref(k, i) * b_ref(k, j);
170 for (j = 1; j <= i1; ++j) {
171 for (i = *m; i >= 1; --i) {
172 temp = alpha2 * b_ref(i, j);
174 for (k = i + 1; k <= i2; ++k) {
175 temp -= a_ref(k, i) * b_ref(k, j);
189 if (lsame2(transa,
"N")) {
193 for (j = 1; j <= i1; ++j) {
196 for (i = 1; i <= i2; ++i) {
197 b_ref(i, j) *= alpha2;
202 for (k = 1; k <= i2; ++k) {
203 if (a_ref(k, j) != 0.) {
206 for (i = 1; i <= i3; ++i) {
207 b_ref(i, j) += aref * b_ref(i, k);
214 temp = 1. / a_ref(j, j);
216 for (i = 1; i <= i2; ++i) {
224 for (j = *n; j >= 1; --j) {
227 for (i = 1; i <= i1; ++i) {
228 b_ref(i, j) *= alpha2;
233 for (k = j + 1; k <= i1; ++k) {
234 if (a_ref(k, j) != 0.) {
237 for (i = 1; i <= i2; ++i) {
238 b_ref(i, j) += aref * b_ref(i, k);
245 temp = 1. / a_ref(j, j);
247 for (i = 1; i <= i1; ++i) {
258 for (k = *n; k >= 1; --k) {
260 temp = 1. / a_ref(k, k);
262 for (i = 1; i <= i1; ++i) {
268 for (j = 1; j <= i1; ++j) {
269 if (a_ref(j, k) != 0.) {
272 for (i = 1; i <= i2; ++i) {
273 b_ref(i, j) -= temp * b_ref(i, k);
281 for (i = 1; i <= i1; ++i) {
282 b_ref(i, k) *= alpha2;
290 for (k = 1; k <= i1; ++k) {
292 temp = 1. / a_ref(k, k);
294 for (i = 1; i <= i2; ++i) {
300 for (j = k + 1; j <= i2; ++j) {
301 if (a_ref(j, k) != 0.) {
304 for (i = 1; i <= i3; ++i) {
305 b_ref(i, j) -= temp * b_ref(i, k);
313 for (i = 1; i <= i2; ++i) {
314 b_ref(i, k) *= alpha2;
DSDP uses BLAS and LAPACK for many of its operations.