DSDP
dtrsm2.c
1 #include "dsdplapack.h"
2 
3 typedef long int integer;
4 typedef long int logical;
5 
6 #define max(a,b) ((a) >= (b) ? (a) : (b))
7 #define dmax(a,b) (double)max(a,b)
8 
9 static 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.