DSDP
dtrsm2.c
1#include "dsdplapack.h"
2
3typedef int integer;
4typedef int logical;
5
6#define max(a,b) ((a) >= (b) ? (a) : (b))
7#define dmax(a,b) (double)max(a,b)
8
9static int lsame2(const char c1[], const char c2[]){
10 if (c1[0]==c2[0]) return 1;
11 else return 0;
12}
13
14/* Subroutine */ 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)
17{
18 /* System generated locals */
19 integer a_dim1, a_offset, b_dim1, b_offset, i1, i2, i3;
20 /* Local variables */
21 static integer info;
22 static double temp,alpha2,aref;
23 static integer i, j, k;
24 static logical lside;
25 static integer nrowa;
26 static logical upper;
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]
30
31 a_dim1 = *lda;
32 a_offset = 1 + a_dim1 * 1;
33 a -= a_offset;
34 b_dim1 = *ldb;
35 b_offset = 1 + b_dim1 * 1;
36 b -= b_offset;
37 /* Function Body */
38 lside = lsame2(side, "L");
39 if (lside) {
40 nrowa = *m;
41 } else {
42 nrowa = *n;
43 }
44 nounit = lsame2(diag, "N");
45 upper = lsame2(uplo, "U");
46 info = 0;
47 if (! lside && ! lsame2(side, "R")) {
48 info = 1;
49 } else if (! upper && ! lsame2(uplo, "L")) {
50 info = 2;
51 } else if (! lsame2(transa, "N") && ! lsame2(transa,
52 "T") && ! lsame2(transa, "C")) {
53 info = 3;
54 } else if (! lsame2(diag, "U") && ! lsame2(diag,
55 "N")) {
56 info = 4;
57 } else if (*m < 0) {
58 info = 5;
59 } else if (*n < 0) {
60 info = 6;
61 } else if (*lda < max(1,nrowa)) {
62 info = 9;
63 } else if (*ldb < max(1,*m)) {
64 info = 11;
65 }
66 if (info != 0) {
67 return info;
68 }
69/* Quick return if possible. */
70 if (*n == 0) {
71 return 0;
72 }
73/* And when alpha.eq.zero. */
74 if (*alpha == 0.) {
75 i1 = *n;
76 for (j = 1; j <= i1; ++j) {
77 i2 = *m;
78 for (i = 1; i <= i2; ++i) {
79 b_ref(i, j) = 0.;
80/* L10: */
81 }
82/* L20: */
83 }
84 return 0;
85 }
86/* Start the operations. */
87 if (lside) {
88 if (lsame2(transa, "N")) {
89 /* Form B := alpha*inv( A )*B. */
90 if (upper) {
91 i1 = *n;
92 for (j = 1; j <= i1; ++j) {
93 if (*alpha != 1.) {
94 i2 = *m;
95 alpha2=*alpha;
96 for (i = 1; i <= i2; ++i) {
97 b_ref(i, j) *= alpha2;
98 /* L30: */
99 }
100 }
101 for (k = *m; k >= 1; --k) {
102 if (b_ref(k, j) != 0.) {
103 if (nounit) {
104 b_ref(k, j) /= a_ref(k, k);
105 }
106 i2 = k - 1;
107 aref=-b_ref(k, j);
108 for (i = 1; i <= i2; ++i) {
109 b_ref(i, j) += aref * a_ref(i, k);
110/* L40: */
111 }
112 }
113 /* L50: */
114 }
115 /* L60: */
116 }
117 } else {
118 i1 = *n;
119 for (j = 1; j <= i1; ++j) {
120 if (*alpha != 1.) {
121 i2 = *m;
122 alpha2=*alpha;
123 for (i = 1; i <= i2; ++i) {
124 b_ref(i, j) *= alpha2;
125 /* L70: */
126 }
127 }
128 i2 = *m;
129 for (k = 1; k <= i2; ++k) {
130 if (b_ref(k, j) != 0.) {
131 if (nounit) {
132 b_ref(k, j) /= a_ref(k, k);
133 }
134 i3 = *m;
135 aref= -b_ref(k, j);
136 for (i = k + 1; i <= i3; ++i) {
137 b_ref(i, j) += aref * a_ref(i, k);
138 /* L80: */
139 }
140 }
141 /* L90: */
142 }
143 /* L100: */
144 }
145 }
146 } else {
147 /* Form B := alpha*inv( A' )*B. */
148 if (upper) {
149 i1 = *n;
150 for (j = 1; j <= i1; ++j) {
151 i2 = *m;
152 alpha2=*alpha;
153 for (i = 1; i <= i2; ++i) {
154 temp = alpha2 * b_ref(i, j);
155 i3 = i - 1;
156 for (k = 1; k <= i3; ++k) {
157 temp -= a_ref(k, i) * b_ref(k, j);
158 /* L110: */
159 }
160 if (nounit) {
161 temp /= a_ref(i, i);
162 }
163 b_ref(i, j) = temp;
164/* L120: */
165 }
166 /* L130: */
167 }
168 } else {
169 i1 = *n;
170 for (j = 1; j <= i1; ++j) {
171 for (i = *m; i >= 1; --i) {
172 temp = alpha2 * b_ref(i, j);
173 i2 = *m;
174 for (k = i + 1; k <= i2; ++k) {
175 temp -= a_ref(k, i) * b_ref(k, j);
176 /* L140: */
177 }
178 if (nounit) {
179 temp /= a_ref(i, i);
180 }
181 b_ref(i, j) = temp;
182 /* L150: */
183 }
184/* L160: */
185 }
186 }
187 }
188 } else {
189 if (lsame2(transa, "N")) {
190 /* Form B := alpha*B*inv( A ). */
191 if (upper) {
192 i1 = *n;
193 for (j = 1; j <= i1; ++j) {
194 if (*alpha != 1.) {
195 i2 = *m;
196 for (i = 1; i <= i2; ++i) {
197 b_ref(i, j) *= alpha2;
198 /* L170: */
199 }
200 }
201 i2 = j - 1;
202 for (k = 1; k <= i2; ++k) {
203 if (a_ref(k, j) != 0.) {
204 i3 = *m;
205 aref=-a_ref(k, j);
206 for (i = 1; i <= i3; ++i) {
207 b_ref(i, j) += aref * b_ref(i, k);
208 /* L180: */
209 }
210 }
211 /* L190: */
212 }
213 if (nounit) {
214 temp = 1. / a_ref(j, j);
215 i2 = *m;
216 for (i = 1; i <= i2; ++i) {
217 b_ref(i, j) *= temp;
218 /* L200: */
219 }
220 }
221 /* L210: */
222 }
223 } else {
224 for (j = *n; j >= 1; --j) {
225 if (*alpha != 1.) {
226 i1 = *m;
227 for (i = 1; i <= i1; ++i) {
228 b_ref(i, j) *= alpha2;
229 /* L220: */
230 }
231 }
232 i1 = *n;
233 for (k = j + 1; k <= i1; ++k) {
234 if (a_ref(k, j) != 0.) {
235 i2 = *m;
236 aref=-a_ref(k, j);
237 for (i = 1; i <= i2; ++i) {
238 b_ref(i, j) += aref * b_ref(i, k);
239 /* L230: */
240 }
241 }
242 /* L240: */
243 }
244 if (nounit) {
245 temp = 1. / a_ref(j, j);
246 i1 = *m;
247 for (i = 1; i <= i1; ++i) {
248 b_ref(i, j) *= temp;
249 /* L250: */
250 }
251 }
252 /* L260: */
253 }
254 }
255 } else {
256 /* Form B := alpha*B*inv( A' ). */
257 if (upper) {
258 for (k = *n; k >= 1; --k) {
259 if (nounit) {
260 temp = 1. / a_ref(k, k);
261 i1 = *m;
262 for (i = 1; i <= i1; ++i) {
263 b_ref(i, k) *= temp;
264 /* L270: */
265 }
266 }
267 i1 = k - 1;
268 for (j = 1; j <= i1; ++j) {
269 if (a_ref(j, k) != 0.) {
270 temp = a_ref(j, k);
271 i2 = *m;
272 for (i = 1; i <= i2; ++i) {
273 b_ref(i, j) -= temp * b_ref(i, k);
274 /* L280: */
275 }
276 }
277 /* L290: */
278 }
279 if (*alpha != 1.) {
280 i1 = *m;
281 for (i = 1; i <= i1; ++i) {
282 b_ref(i, k) *= alpha2;
283 /* L300: */
284 }
285 }
286 /* L310: */
287 }
288 } else {
289 i1 = *n;
290 for (k = 1; k <= i1; ++k) {
291 if (nounit) {
292 temp = 1. / a_ref(k, k);
293 i2 = *m;
294 for (i = 1; i <= i2; ++i) {
295 b_ref(i, k) *= temp;
296 /* L320: */
297 }
298 }
299 i2 = *n;
300 for (j = k + 1; j <= i2; ++j) {
301 if (a_ref(j, k) != 0.) {
302 temp = a_ref(j, k);
303 i3 = *m;
304 for (i = 1; i <= i3; ++i) {
305 b_ref(i, j) -= temp * b_ref(i, k);
306 /* L330: */
307 }
308 }
309 /* L340: */
310 }
311 if (*alpha != 1.) {
312 i2 = *m;
313 for (i = 1; i <= i2; ++i) {
314 b_ref(i, k) *= alpha2;
315 /* L350: */
316 }
317 }
318 /* L360: */
319 }
320 }
321 }
322 }
323 return 0;
324 /* End of DTRSM . */
325} /* dtrsm_ */
326#undef b_ref
327#undef a_ref
328
DSDP uses BLAS and LAPACK for many of its operations.