![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas|
Return to dtrsm.f CVS log | Up to [local] / rpl / lapack / blas |
1.1 bertrand 1: SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA
4: INTEGER LDA,LDB,M,N
5: CHARACTER DIAG,SIDE,TRANSA,UPLO
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),B(LDB,*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DTRSM solves one of the matrix equations
15: *
16: * op( A )*X = alpha*B, or X*op( A ) = alpha*B,
17: *
18: * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
19: * non-unit, upper or lower triangular matrix and op( A ) is one of
20: *
21: * op( A ) = A or op( A ) = A'.
22: *
23: * The matrix X is overwritten on B.
24: *
25: * Arguments
26: * ==========
27: *
28: * SIDE - CHARACTER*1.
29: * On entry, SIDE specifies whether op( A ) appears on the left
30: * or right of X as follows:
31: *
32: * SIDE = 'L' or 'l' op( A )*X = alpha*B.
33: *
34: * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
35: *
36: * Unchanged on exit.
37: *
38: * UPLO - CHARACTER*1.
39: * On entry, UPLO specifies whether the matrix A is an upper or
40: * lower triangular matrix as follows:
41: *
42: * UPLO = 'U' or 'u' A is an upper triangular matrix.
43: *
44: * UPLO = 'L' or 'l' A is a lower triangular matrix.
45: *
46: * Unchanged on exit.
47: *
48: * TRANSA - CHARACTER*1.
49: * On entry, TRANSA specifies the form of op( A ) to be used in
50: * the matrix multiplication as follows:
51: *
52: * TRANSA = 'N' or 'n' op( A ) = A.
53: *
54: * TRANSA = 'T' or 't' op( A ) = A'.
55: *
56: * TRANSA = 'C' or 'c' op( A ) = A'.
57: *
58: * Unchanged on exit.
59: *
60: * DIAG - CHARACTER*1.
61: * On entry, DIAG specifies whether or not A is unit triangular
62: * as follows:
63: *
64: * DIAG = 'U' or 'u' A is assumed to be unit triangular.
65: *
66: * DIAG = 'N' or 'n' A is not assumed to be unit
67: * triangular.
68: *
69: * Unchanged on exit.
70: *
71: * M - INTEGER.
72: * On entry, M specifies the number of rows of B. M must be at
73: * least zero.
74: * Unchanged on exit.
75: *
76: * N - INTEGER.
77: * On entry, N specifies the number of columns of B. N must be
78: * at least zero.
79: * Unchanged on exit.
80: *
81: * ALPHA - DOUBLE PRECISION.
82: * On entry, ALPHA specifies the scalar alpha. When alpha is
83: * zero then A is not referenced and B need not be set before
84: * entry.
85: * Unchanged on exit.
86: *
87: * A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
88: * when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'.
89: * Before entry with UPLO = 'U' or 'u', the leading k by k
90: * upper triangular part of the array A must contain the upper
91: * triangular matrix and the strictly lower triangular part of
92: * A is not referenced.
93: * Before entry with UPLO = 'L' or 'l', the leading k by k
94: * lower triangular part of the array A must contain the lower
95: * triangular matrix and the strictly upper triangular part of
96: * A is not referenced.
97: * Note that when DIAG = 'U' or 'u', the diagonal elements of
98: * A are not referenced either, but are assumed to be unity.
99: * Unchanged on exit.
100: *
101: * LDA - INTEGER.
102: * On entry, LDA specifies the first dimension of A as declared
103: * in the calling (sub) program. When SIDE = 'L' or 'l' then
104: * LDA must be at least max( 1, m ), when SIDE = 'R' or 'r'
105: * then LDA must be at least max( 1, n ).
106: * Unchanged on exit.
107: *
108: * B - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
109: * Before entry, the leading m by n part of the array B must
110: * contain the right-hand side matrix B, and on exit is
111: * overwritten by the solution matrix X.
112: *
113: * LDB - INTEGER.
114: * On entry, LDB specifies the first dimension of B as declared
115: * in the calling (sub) program. LDB must be at least
116: * max( 1, m ).
117: * Unchanged on exit.
118: *
119: * Further Details
120: * ===============
121: *
122: * Level 3 Blas routine.
123: *
124: *
125: * -- Written on 8-February-1989.
126: * Jack Dongarra, Argonne National Laboratory.
127: * Iain Duff, AERE Harwell.
128: * Jeremy Du Croz, Numerical Algorithms Group Ltd.
129: * Sven Hammarling, Numerical Algorithms Group Ltd.
130: *
131: * =====================================================================
132: *
133: * .. External Functions ..
134: LOGICAL LSAME
135: EXTERNAL LSAME
136: * ..
137: * .. External Subroutines ..
138: EXTERNAL XERBLA
139: * ..
140: * .. Intrinsic Functions ..
141: INTRINSIC MAX
142: * ..
143: * .. Local Scalars ..
144: DOUBLE PRECISION TEMP
145: INTEGER I,INFO,J,K,NROWA
146: LOGICAL LSIDE,NOUNIT,UPPER
147: * ..
148: * .. Parameters ..
149: DOUBLE PRECISION ONE,ZERO
150: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
151: * ..
152: *
153: * Test the input parameters.
154: *
155: LSIDE = LSAME(SIDE,'L')
156: IF (LSIDE) THEN
157: NROWA = M
158: ELSE
159: NROWA = N
160: END IF
161: NOUNIT = LSAME(DIAG,'N')
162: UPPER = LSAME(UPLO,'U')
163: *
164: INFO = 0
165: IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
166: INFO = 1
167: ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
168: INFO = 2
169: ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
170: + (.NOT.LSAME(TRANSA,'T')) .AND.
171: + (.NOT.LSAME(TRANSA,'C'))) THEN
172: INFO = 3
173: ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
174: INFO = 4
175: ELSE IF (M.LT.0) THEN
176: INFO = 5
177: ELSE IF (N.LT.0) THEN
178: INFO = 6
179: ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
180: INFO = 9
181: ELSE IF (LDB.LT.MAX(1,M)) THEN
182: INFO = 11
183: END IF
184: IF (INFO.NE.0) THEN
185: CALL XERBLA('DTRSM ',INFO)
186: RETURN
187: END IF
188: *
189: * Quick return if possible.
190: *
191: IF (M.EQ.0 .OR. N.EQ.0) RETURN
192: *
193: * And when alpha.eq.zero.
194: *
195: IF (ALPHA.EQ.ZERO) THEN
196: DO 20 J = 1,N
197: DO 10 I = 1,M
198: B(I,J) = ZERO
199: 10 CONTINUE
200: 20 CONTINUE
201: RETURN
202: END IF
203: *
204: * Start the operations.
205: *
206: IF (LSIDE) THEN
207: IF (LSAME(TRANSA,'N')) THEN
208: *
209: * Form B := alpha*inv( A )*B.
210: *
211: IF (UPPER) THEN
212: DO 60 J = 1,N
213: IF (ALPHA.NE.ONE) THEN
214: DO 30 I = 1,M
215: B(I,J) = ALPHA*B(I,J)
216: 30 CONTINUE
217: END IF
218: DO 50 K = M,1,-1
219: IF (B(K,J).NE.ZERO) THEN
220: IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
221: DO 40 I = 1,K - 1
222: B(I,J) = B(I,J) - B(K,J)*A(I,K)
223: 40 CONTINUE
224: END IF
225: 50 CONTINUE
226: 60 CONTINUE
227: ELSE
228: DO 100 J = 1,N
229: IF (ALPHA.NE.ONE) THEN
230: DO 70 I = 1,M
231: B(I,J) = ALPHA*B(I,J)
232: 70 CONTINUE
233: END IF
234: DO 90 K = 1,M
235: IF (B(K,J).NE.ZERO) THEN
236: IF (NOUNIT) B(K,J) = B(K,J)/A(K,K)
237: DO 80 I = K + 1,M
238: B(I,J) = B(I,J) - B(K,J)*A(I,K)
239: 80 CONTINUE
240: END IF
241: 90 CONTINUE
242: 100 CONTINUE
243: END IF
244: ELSE
245: *
246: * Form B := alpha*inv( A' )*B.
247: *
248: IF (UPPER) THEN
249: DO 130 J = 1,N
250: DO 120 I = 1,M
251: TEMP = ALPHA*B(I,J)
252: DO 110 K = 1,I - 1
253: TEMP = TEMP - A(K,I)*B(K,J)
254: 110 CONTINUE
255: IF (NOUNIT) TEMP = TEMP/A(I,I)
256: B(I,J) = TEMP
257: 120 CONTINUE
258: 130 CONTINUE
259: ELSE
260: DO 160 J = 1,N
261: DO 150 I = M,1,-1
262: TEMP = ALPHA*B(I,J)
263: DO 140 K = I + 1,M
264: TEMP = TEMP - A(K,I)*B(K,J)
265: 140 CONTINUE
266: IF (NOUNIT) TEMP = TEMP/A(I,I)
267: B(I,J) = TEMP
268: 150 CONTINUE
269: 160 CONTINUE
270: END IF
271: END IF
272: ELSE
273: IF (LSAME(TRANSA,'N')) THEN
274: *
275: * Form B := alpha*B*inv( A ).
276: *
277: IF (UPPER) THEN
278: DO 210 J = 1,N
279: IF (ALPHA.NE.ONE) THEN
280: DO 170 I = 1,M
281: B(I,J) = ALPHA*B(I,J)
282: 170 CONTINUE
283: END IF
284: DO 190 K = 1,J - 1
285: IF (A(K,J).NE.ZERO) THEN
286: DO 180 I = 1,M
287: B(I,J) = B(I,J) - A(K,J)*B(I,K)
288: 180 CONTINUE
289: END IF
290: 190 CONTINUE
291: IF (NOUNIT) THEN
292: TEMP = ONE/A(J,J)
293: DO 200 I = 1,M
294: B(I,J) = TEMP*B(I,J)
295: 200 CONTINUE
296: END IF
297: 210 CONTINUE
298: ELSE
299: DO 260 J = N,1,-1
300: IF (ALPHA.NE.ONE) THEN
301: DO 220 I = 1,M
302: B(I,J) = ALPHA*B(I,J)
303: 220 CONTINUE
304: END IF
305: DO 240 K = J + 1,N
306: IF (A(K,J).NE.ZERO) THEN
307: DO 230 I = 1,M
308: B(I,J) = B(I,J) - A(K,J)*B(I,K)
309: 230 CONTINUE
310: END IF
311: 240 CONTINUE
312: IF (NOUNIT) THEN
313: TEMP = ONE/A(J,J)
314: DO 250 I = 1,M
315: B(I,J) = TEMP*B(I,J)
316: 250 CONTINUE
317: END IF
318: 260 CONTINUE
319: END IF
320: ELSE
321: *
322: * Form B := alpha*B*inv( A' ).
323: *
324: IF (UPPER) THEN
325: DO 310 K = N,1,-1
326: IF (NOUNIT) THEN
327: TEMP = ONE/A(K,K)
328: DO 270 I = 1,M
329: B(I,K) = TEMP*B(I,K)
330: 270 CONTINUE
331: END IF
332: DO 290 J = 1,K - 1
333: IF (A(J,K).NE.ZERO) THEN
334: TEMP = A(J,K)
335: DO 280 I = 1,M
336: B(I,J) = B(I,J) - TEMP*B(I,K)
337: 280 CONTINUE
338: END IF
339: 290 CONTINUE
340: IF (ALPHA.NE.ONE) THEN
341: DO 300 I = 1,M
342: B(I,K) = ALPHA*B(I,K)
343: 300 CONTINUE
344: END IF
345: 310 CONTINUE
346: ELSE
347: DO 360 K = 1,N
348: IF (NOUNIT) THEN
349: TEMP = ONE/A(K,K)
350: DO 320 I = 1,M
351: B(I,K) = TEMP*B(I,K)
352: 320 CONTINUE
353: END IF
354: DO 340 J = K + 1,N
355: IF (A(J,K).NE.ZERO) THEN
356: TEMP = A(J,K)
357: DO 330 I = 1,M
358: B(I,J) = B(I,J) - TEMP*B(I,K)
359: 330 CONTINUE
360: END IF
361: 340 CONTINUE
362: IF (ALPHA.NE.ONE) THEN
363: DO 350 I = 1,M
364: B(I,K) = ALPHA*B(I,K)
365: 350 CONTINUE
366: END IF
367: 360 CONTINUE
368: END IF
369: END IF
370: END IF
371: *
372: RETURN
373: *
374: * End of DTRSM .
375: *
376: END