3 /* Subroutine */ int dtrsm_(char *side, char *uplo, char *transa, char *diag,
4 integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
5 lda, doublereal *b, integer *ldb)
7 /* System generated locals */
8 integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
11 integer i__, j, k, info;
14 extern logical lsame_(char *, char *);
17 extern /* Subroutine */ int xerbla_(char *, integer *);
20 /* .. Scalar Arguments .. */
22 /* .. Array Arguments .. */
28 /* DTRSM solves one of the matrix equations */
30 /* op( A )*X = alpha*B, or X*op( A ) = alpha*B, */
32 /* where alpha is a scalar, X and B are m by n matrices, A is a unit, or */
33 /* non-unit, upper or lower triangular matrix and op( A ) is one of */
35 /* op( A ) = A or op( A ) = A'. */
37 /* The matrix X is overwritten on B. */
42 /* SIDE - CHARACTER*1. */
43 /* On entry, SIDE specifies whether op( A ) appears on the left */
44 /* or right of X as follows: */
46 /* SIDE = 'L' or 'l' op( A )*X = alpha*B. */
48 /* SIDE = 'R' or 'r' X*op( A ) = alpha*B. */
50 /* Unchanged on exit. */
52 /* UPLO - CHARACTER*1. */
53 /* On entry, UPLO specifies whether the matrix A is an upper or */
54 /* lower triangular matrix as follows: */
56 /* UPLO = 'U' or 'u' A is an upper triangular matrix. */
58 /* UPLO = 'L' or 'l' A is a lower triangular matrix. */
60 /* Unchanged on exit. */
62 /* TRANSA - CHARACTER*1. */
63 /* On entry, TRANSA specifies the form of op( A ) to be used in */
64 /* the matrix multiplication as follows: */
66 /* TRANSA = 'N' or 'n' op( A ) = A. */
68 /* TRANSA = 'T' or 't' op( A ) = A'. */
70 /* TRANSA = 'C' or 'c' op( A ) = A'. */
72 /* Unchanged on exit. */
74 /* DIAG - CHARACTER*1. */
75 /* On entry, DIAG specifies whether or not A is unit triangular */
78 /* DIAG = 'U' or 'u' A is assumed to be unit triangular. */
80 /* DIAG = 'N' or 'n' A is not assumed to be unit */
83 /* Unchanged on exit. */
86 /* On entry, M specifies the number of rows of B. M must be at */
88 /* Unchanged on exit. */
91 /* On entry, N specifies the number of columns of B. N must be */
93 /* Unchanged on exit. */
95 /* ALPHA - DOUBLE PRECISION. */
96 /* On entry, ALPHA specifies the scalar alpha. When alpha is */
97 /* zero then A is not referenced and B need not be set before */
99 /* Unchanged on exit. */
101 /* A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m */
102 /* when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. */
103 /* Before entry with UPLO = 'U' or 'u', the leading k by k */
104 /* upper triangular part of the array A must contain the upper */
105 /* triangular matrix and the strictly lower triangular part of */
106 /* A is not referenced. */
107 /* Before entry with UPLO = 'L' or 'l', the leading k by k */
108 /* lower triangular part of the array A must contain the lower */
109 /* triangular matrix and the strictly upper triangular part of */
110 /* A is not referenced. */
111 /* Note that when DIAG = 'U' or 'u', the diagonal elements of */
112 /* A are not referenced either, but are assumed to be unity. */
113 /* Unchanged on exit. */
116 /* On entry, LDA specifies the first dimension of A as declared */
117 /* in the calling (sub) program. When SIDE = 'L' or 'l' then */
118 /* LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' */
119 /* then LDA must be at least max( 1, n ). */
120 /* Unchanged on exit. */
122 /* B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
123 /* Before entry, the leading m by n part of the array B must */
124 /* contain the right-hand side matrix B, and on exit is */
125 /* overwritten by the solution matrix X. */
128 /* On entry, LDB specifies the first dimension of B as declared */
129 /* in the calling (sub) program. LDB must be at least */
131 /* Unchanged on exit. */
134 /* Level 3 Blas routine. */
137 /* -- Written on 8-February-1989. */
138 /* Jack Dongarra, Argonne National Laboratory. */
139 /* Iain Duff, AERE Harwell. */
140 /* Jeremy Du Croz, Numerical Algorithms Group Ltd. */
141 /* Sven Hammarling, Numerical Algorithms Group Ltd. */
144 /* .. External Functions .. */
146 /* .. External Subroutines .. */
148 /* .. Intrinsic Functions .. */
150 /* .. Local Scalars .. */
152 /* .. Parameters .. */
155 /* Test the input parameters. */
157 /* Parameter adjustments */
159 a_offset = 1 + a_dim1;
162 b_offset = 1 + b_dim1;
166 lside = lsame_(side, "L");
172 nounit = lsame_(diag, "N");
173 upper = lsame_(uplo, "U");
176 if (! lside && ! lsame_(side, "R")) {
178 } else if (! upper && ! lsame_(uplo, "L")) {
180 } else if (! lsame_(transa, "N") && ! lsame_(transa,
181 "T") && ! lsame_(transa, "C")) {
183 } else if (! lsame_(diag, "U") && ! lsame_(diag,
190 } else if (*lda < max(1,nrowa)) {
192 } else if (*ldb < max(1,*m)) {
196 xerbla_("DTRSM ", &info);
200 /* Quick return if possible. */
206 /* And when alpha.eq.zero. */
210 for (j = 1; j <= i__1; ++j) {
212 for (i__ = 1; i__ <= i__2; ++i__) {
213 b[i__ + j * b_dim1] = 0.;
221 /* Start the operations. */
224 if (lsame_(transa, "N")) {
226 /* Form B := alpha*inv( A )*B. */
230 for (j = 1; j <= i__1; ++j) {
233 for (i__ = 1; i__ <= i__2; ++i__) {
234 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
239 for (k = *m; k >= 1; --k) {
240 if (b[k + j * b_dim1] != 0.) {
242 b[k + j * b_dim1] /= a[k + k * a_dim1];
245 for (i__ = 1; i__ <= i__2; ++i__) {
246 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
257 for (j = 1; j <= i__1; ++j) {
260 for (i__ = 1; i__ <= i__2; ++i__) {
261 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
267 for (k = 1; k <= i__2; ++k) {
268 if (b[k + j * b_dim1] != 0.) {
270 b[k + j * b_dim1] /= a[k + k * a_dim1];
273 for (i__ = k + 1; i__ <= i__3; ++i__) {
274 b[i__ + j * b_dim1] -= b[k + j * b_dim1] * a[
286 /* Form B := alpha*inv( A' )*B. */
290 for (j = 1; j <= i__1; ++j) {
292 for (i__ = 1; i__ <= i__2; ++i__) {
293 temp = *alpha * b[i__ + j * b_dim1];
295 for (k = 1; k <= i__3; ++k) {
296 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
300 temp /= a[i__ + i__ * a_dim1];
302 b[i__ + j * b_dim1] = temp;
309 for (j = 1; j <= i__1; ++j) {
310 for (i__ = *m; i__ >= 1; --i__) {
311 temp = *alpha * b[i__ + j * b_dim1];
313 for (k = i__ + 1; k <= i__2; ++k) {
314 temp -= a[k + i__ * a_dim1] * b[k + j * b_dim1];
318 temp /= a[i__ + i__ * a_dim1];
320 b[i__ + j * b_dim1] = temp;
328 if (lsame_(transa, "N")) {
330 /* Form B := alpha*B*inv( A ). */
334 for (j = 1; j <= i__1; ++j) {
337 for (i__ = 1; i__ <= i__2; ++i__) {
338 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
344 for (k = 1; k <= i__2; ++k) {
345 if (a[k + j * a_dim1] != 0.) {
347 for (i__ = 1; i__ <= i__3; ++i__) {
348 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
356 temp = 1. / a[j + j * a_dim1];
358 for (i__ = 1; i__ <= i__2; ++i__) {
359 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
366 for (j = *n; j >= 1; --j) {
369 for (i__ = 1; i__ <= i__1; ++i__) {
370 b[i__ + j * b_dim1] = *alpha * b[i__ + j * b_dim1]
376 for (k = j + 1; k <= i__1; ++k) {
377 if (a[k + j * a_dim1] != 0.) {
379 for (i__ = 1; i__ <= i__2; ++i__) {
380 b[i__ + j * b_dim1] -= a[k + j * a_dim1] * b[
388 temp = 1. / a[j + j * a_dim1];
390 for (i__ = 1; i__ <= i__1; ++i__) {
391 b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
400 /* Form B := alpha*B*inv( A' ). */
403 for (k = *n; k >= 1; --k) {
405 temp = 1. / a[k + k * a_dim1];
407 for (i__ = 1; i__ <= i__1; ++i__) {
408 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
413 for (j = 1; j <= i__1; ++j) {
414 if (a[j + k * a_dim1] != 0.) {
415 temp = a[j + k * a_dim1];
417 for (i__ = 1; i__ <= i__2; ++i__) {
418 b[i__ + j * b_dim1] -= temp * b[i__ + k *
427 for (i__ = 1; i__ <= i__1; ++i__) {
428 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]
437 for (k = 1; k <= i__1; ++k) {
439 temp = 1. / a[k + k * a_dim1];
441 for (i__ = 1; i__ <= i__2; ++i__) {
442 b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
447 for (j = k + 1; j <= i__2; ++j) {
448 if (a[j + k * a_dim1] != 0.) {
449 temp = a[j + k * a_dim1];
451 for (i__ = 1; i__ <= i__3; ++i__) {
452 b[i__ + j * b_dim1] -= temp * b[i__ + k *
461 for (i__ = 1; i__ <= i__2; ++i__) {
462 b[i__ + k * b_dim1] = *alpha * b[i__ + k * b_dim1]