1 /* dtrsm.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14 #include "blaswrap.h"
15
dtrsm_(char * side,char * uplo,char * transa,char * diag,integer * m,integer * n,doublereal * alpha,doublereal * a,integer * lda,doublereal * b,integer * ldb)16 /* Subroutine */ int dtrsm_(char *side, char *uplo, char *transa, char *diag,
17 integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
18 lda, doublereal *b, integer *ldb)
19 {
20 /* System generated locals */
21 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
22
23 /* Local variables */
24 integer i__, j, k, info;
25 doublereal temp;
26 logical lside;
27 extern logical lsame_(char *, char *);
28 integer nrowa;
29 logical upper;
30 extern /* Subroutine */ int xerbla_(char *, integer *);
31 logical nounit;
32
33 /* .. Scalar Arguments .. */
34 /* .. */
35 /* .. Array Arguments .. */
36 /* .. */
37
38 /* Purpose */
39 /* ======= */
40
41 /* DTRSM solves one of the matrix equations */
42
43 /* op( A )*X = alpha*B, or X*op( A ) = alpha*B, */
44
45 /* where alpha is a scalar, X and B are m by n matrices, A is a unit, or */
46 /* non-unit, upper or lower triangular matrix and op( A ) is one of */
47
48 /* op( A ) = A or op( A ) = A'. */
49
50 /* The matrix X is overwritten on B. */
51
52 /* Arguments */
53 /* ========== */
54
55 /* SIDE - CHARACTER*1. */
56 /* On entry, SIDE specifies whether op( A ) appears on the left */
57 /* or right of X as follows: */
58
59 /* SIDE = 'L' or 'l' op( A )*X = alpha*B. */
60
61 /* SIDE = 'R' or 'r' X*op( A ) = alpha*B. */
62
63 /* Unchanged on exit. */
64
65 /* UPLO - CHARACTER*1. */
66 /* On entry, UPLO specifies whether the matrix A is an upper or */
67 /* lower triangular matrix as follows: */
68
69 /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
70
71 /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
72
73 /* Unchanged on exit. */
74
75 /* TRANSA - CHARACTER*1. */
76 /* On entry, TRANSA specifies the form of op( A ) to be used in */
77 /* the matrix multiplication as follows: */
78
79 /* TRANSA = 'N' or 'n' op( A ) = A. */
80
81 /* TRANSA = 'T' or 't' op( A ) = A'. */
82
83 /* TRANSA = 'C' or 'c' op( A ) = A'. */
84
85 /* Unchanged on exit. */
86
87 /* DIAG - CHARACTER*1. */
88 /* On entry, DIAG specifies whether or not A is unit triangular */
89 /* as follows: */
90
91 /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
92
93 /* DIAG = 'N' or 'n' A is not assumed to be unit */
94 /* triangular. */
95
96 /* Unchanged on exit. */
97
98 /* M - INTEGER. */
99 /* On entry, M specifies the number of rows of B. M must be at */
100 /* least zero. */
101 /* Unchanged on exit. */
102
103 /* N - INTEGER. */
104 /* On entry, N specifies the number of columns of B. N must be */
105 /* at least zero. */
106 /* Unchanged on exit. */
107
108 /* ALPHA - DOUBLE PRECISION. */
109 /* On entry, ALPHA specifies the scalar alpha. When alpha is */
110 /* zero then A is not referenced and B need not be set before */
111 /* entry. */
112 /* Unchanged on exit. */
113
114 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m */
115 /* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. */
116 /* Before entry with UPLO = 'U' or 'u', the leading k by k */
117 /* upper triangular part of the array A must contain the upper */
118 /* triangular matrix and the strictly lower triangular part of */
119 /* A is not referenced. */
120 /* Before entry with UPLO = 'L' or 'l', the leading k by k */
121 /* lower triangular part of the array A must contain the lower */
122 /* triangular matrix and the strictly upper triangular part of */
123 /* A is not referenced. */
124 /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
125 /* A are not referenced either, but are assumed to be unity. */
126 /* Unchanged on exit. */
127
128 /* LDA - INTEGER. */
129 /* On entry, LDA specifies the first dimension of A as declared */
130 /* in the calling (sub) program. When SIDE = 'L' or 'l' then */
131 /* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' */
132 /* then LDA must be at least max( 1, n ). */
133 /* Unchanged on exit. */
134
135 /* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
136 /* Before entry, the leading m by n part of the array B must */
137 /* contain the right-hand side matrix B, and on exit is */
138 /* overwritten by the solution matrix X. */
139
140 /* LDB - INTEGER. */
141 /* On entry, LDB specifies the first dimension of B as declared */
142 /* in the calling (sub) program. LDB must be at least */
143 /* max( 1, m ). */
144 /* Unchanged on exit. */
145
146
147 /* Level 3 Blas routine. */
148
149
150 /* -- Written on 8-February-1989. */
151 /* Jack Dongarra, Argonne National Laboratory. */
152 /* Iain Duff, AERE Harwell. */
153 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
154 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
155
156
157 /* .. External Functions .. */
158 /* .. */
159 /* .. External Subroutines .. */
160 /* .. */
161 /* .. Intrinsic Functions .. */
162 /* .. */
163 /* .. Local Scalars .. */
164 /* .. */
165 /* .. Parameters .. */
166 /* .. */
167
168 /* Test the input parameters. */
169
170 /* Parameter adjustments */
171 a_dim1 = *lda;
172 a_offset = 1 + a_dim1;
173 a -= a_offset;
174 b_dim1 = *ldb;
175 b_offset = 1 + b_dim1;
176 b -= b_offset;
177
178 /* Function Body */
179 lside = lsame_(side, "L");
180 if (lside) {
181 nrowa = *m;
182 } else {
183 nrowa = *n;
184 }
185 nounit = lsame_(diag, "N");
186 upper = lsame_(uplo, "U");
187
188 info = 0;
189 if (! lside && ! lsame_(side, "R")) {
190 info = 1;
191 } else if (! upper && ! lsame_(uplo, "L")) {
192 info = 2;
193 } else if (! lsame_(transa, "N") && ! lsame_(transa,
194 "T") && ! lsame_(transa, "C")) {
195 info = 3;
196 } else if (! lsame_(diag, "U") && ! lsame_(diag,
197 "N")) {
198 info = 4;
199 } else if (*m < 0) {
200 info = 5;
201 } else if (*n < 0) {
202 info = 6;
203 } else if (*lda < max(1,nrowa)) {
204 info = 9;
205 } else if (*ldb < max(1,*m)) {
206 info = 11;
207 }
208 if (info != 0) {
209 xerbla_("DTRSM ", &info);
210 return 0;
211 }
212
213 /* Quick return if possible. */
214
215 if (*m == 0 || *n == 0) {
216 return 0;
217 }
218
219 /* And when alpha.eq.zero. */
220
221 if (*alpha == 0.) {
222 i__1 = *n;
223 for (j = 1; j <= i__1; ++j) {
224 i__2 = *m;
225 for (i__ = 1; i__ <= i__2; ++i__) {
226 b[i__ + j * b_dim1] = 0.;
227 /* L10: */
228 }
229 /* L20: */
230 }
231 return 0;
232 }
233
234 /* Start the operations. */
235
236 if (lside) {
237 if (lsame_(transa, "N")) {
238
239 /* Form B := alpha*inv( A )*B. */
240
241 if (upper) {
242 i__1 = *n;
243 for (j = 1; j <= i__1; ++j) {
244 if (*alpha != 1.) {
245 i__2 = *m;
246 for (i__ = 1; i__ <= i__2; ++i__) {
247 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
248 ;
249 /* L30: */
250 }
251 }
252 for (k = *m; k >= 1; --k) {
253 if (b[k + j * b_dim1] != 0.) {
254 if (nounit) {
255 b[k + j * b_dim1] /= a[k + k * a_dim1];
256 }
257 i__2 = k - 1;
258 for (i__ = 1; i__ <= i__2; ++i__) {
259 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
260 i__ + k * a_dim1];
261 /* L40: */
262 }
263 }
264 /* L50: */
265 }
266 /* L60: */
267 }
268 } else {
269 i__1 = *n;
270 for (j = 1; j <= i__1; ++j) {
271 if (*alpha != 1.) {
272 i__2 = *m;
273 for (i__ = 1; i__ <= i__2; ++i__) {
274 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
275 ;
276 /* L70: */
277 }
278 }
279 i__2 = *m;
280 for (k = 1; k <= i__2; ++k) {
281 if (b[k + j * b_dim1] != 0.) {
282 if (nounit) {
283 b[k + j * b_dim1] /= a[k + k * a_dim1];
284 }
285 i__3 = *m;
286 for (i__ = k + 1; i__ <= i__3; ++i__) {
287 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
288 i__ + k * a_dim1];
289 /* L80: */
290 }
291 }
292 /* L90: */
293 }
294 /* L100: */
295 }
296 }
297 } else {
298
299 /* Form B := alpha*inv( A' )*B. */
300
301 if (upper) {
302 i__1 = *n;
303 for (j = 1; j <= i__1; ++j) {
304 i__2 = *m;
305 for (i__ = 1; i__ <= i__2; ++i__) {
306 temp = *alpha * b[i__ + j * b_dim1];
307 i__3 = i__ - 1;
308 for (k = 1; k <= i__3; ++k) {
309 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
310 /* L110: */
311 }
312 if (nounit) {
313 temp /= a[i__ + i__ * a_dim1];
314 }
315 b[i__ + j * b_dim1] = temp;
316 /* L120: */
317 }
318 /* L130: */
319 }
320 } else {
321 i__1 = *n;
322 for (j = 1; j <= i__1; ++j) {
323 for (i__ = *m; i__ >= 1; --i__) {
324 temp = *alpha * b[i__ + j * b_dim1];
325 i__2 = *m;
326 for (k = i__ + 1; k <= i__2; ++k) {
327 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
328 /* L140: */
329 }
330 if (nounit) {
331 temp /= a[i__ + i__ * a_dim1];
332 }
333 b[i__ + j * b_dim1] = temp;
334 /* L150: */
335 }
336 /* L160: */
337 }
338 }
339 }
340 } else {
341 if (lsame_(transa, "N")) {
342
343 /* Form B := alpha*B*inv( A ). */
344
345 if (upper) {
346 i__1 = *n;
347 for (j = 1; j <= i__1; ++j) {
348 if (*alpha != 1.) {
349 i__2 = *m;
350 for (i__ = 1; i__ <= i__2; ++i__) {
351 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
352 ;
353 /* L170: */
354 }
355 }
356 i__2 = j - 1;
357 for (k = 1; k <= i__2; ++k) {
358 if (a[k + j * a_dim1] != 0.) {
359 i__3 = *m;
360 for (i__ = 1; i__ <= i__3; ++i__) {
361 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
362 i__ + k * b_dim1];
363 /* L180: */
364 }
365 }
366 /* L190: */
367 }
368 if (nounit) {
369 temp = 1. / a[j + j * a_dim1];
370 i__2 = *m;
371 for (i__ = 1; i__ <= i__2; ++i__) {
372 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
373 /* L200: */
374 }
375 }
376 /* L210: */
377 }
378 } else {
379 for (j = *n; j >= 1; --j) {
380 if (*alpha != 1.) {
381 i__1 = *m;
382 for (i__ = 1; i__ <= i__1; ++i__) {
383 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
384 ;
385 /* L220: */
386 }
387 }
388 i__1 = *n;
389 for (k = j + 1; k <= i__1; ++k) {
390 if (a[k + j * a_dim1] != 0.) {
391 i__2 = *m;
392 for (i__ = 1; i__ <= i__2; ++i__) {
393 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
394 i__ + k * b_dim1];
395 /* L230: */
396 }
397 }
398 /* L240: */
399 }
400 if (nounit) {
401 temp = 1. / a[j + j * a_dim1];
402 i__1 = *m;
403 for (i__ = 1; i__ <= i__1; ++i__) {
404 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
405 /* L250: */
406 }
407 }
408 /* L260: */
409 }
410 }
411 } else {
412
413 /* Form B := alpha*B*inv( A' ). */
414
415 if (upper) {
416 for (k = *n; k >= 1; --k) {
417 if (nounit) {
418 temp = 1. / a[k + k * a_dim1];
419 i__1 = *m;
420 for (i__ = 1; i__ <= i__1; ++i__) {
421 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
422 /* L270: */
423 }
424 }
425 i__1 = k - 1;
426 for (j = 1; j <= i__1; ++j) {
427 if (a[j + k * a_dim1] != 0.) {
428 temp = a[j + k * a_dim1];
429 i__2 = *m;
430 for (i__ = 1; i__ <= i__2; ++i__) {
431 b[i__ + j * b_dim1] -= temp * b[i__ + k *
432 b_dim1];
433 /* L280: */
434 }
435 }
436 /* L290: */
437 }
438 if (*alpha != 1.) {
439 i__1 = *m;
440 for (i__ = 1; i__ <= i__1; ++i__) {
441 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
442 ;
443 /* L300: */
444 }
445 }
446 /* L310: */
447 }
448 } else {
449 i__1 = *n;
450 for (k = 1; k <= i__1; ++k) {
451 if (nounit) {
452 temp = 1. / a[k + k * a_dim1];
453 i__2 = *m;
454 for (i__ = 1; i__ <= i__2; ++i__) {
455 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
456 /* L320: */
457 }
458 }
459 i__2 = *n;
460 for (j = k + 1; j <= i__2; ++j) {
461 if (a[j + k * a_dim1] != 0.) {
462 temp = a[j + k * a_dim1];
463 i__3 = *m;
464 for (i__ = 1; i__ <= i__3; ++i__) {
465 b[i__ + j * b_dim1] -= temp * b[i__ + k *
466 b_dim1];
467 /* L330: */
468 }
469 }
470 /* L340: */
471 }
472 if (*alpha != 1.) {
473 i__2 = *m;
474 for (i__ = 1; i__ <= i__2; ++i__) {
475 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
476 ;
477 /* L350: */
478 }
479 }
480 /* L360: */
481 }
482 }
483 }
484 }
485
486 return 0;
487
488 /* End of DTRSM . */
489
490 } /* dtrsm_ */
491