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

Annotation of /trunk/blas/dtrsm.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (15 years, 5 months ago) by aw0a
File size: 10615 byte(s)
Setting up web subdirectory in repository
1 aw0a 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