MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
dtrsm.f
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**T.
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**T.
55 *
56 * TRANSA = 'C' or 'c' op( A ) = A**T.
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**T )*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**T ).
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