← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/blas/dtrmm.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

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