/[ascend]/trunk/blas/dtrsm.f
ViewVC logotype

Contents of /trunk/blas/dtrsm.f

Parent Directory Parent Directory | Revision Log Revision Log


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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22