source: trunk/SOURCES/BLAS/strmm.f @ 22

Last change on this file since 22 was 22, checked in by roche, 9 years ago

Petites adaptations diverses du code pour compilation en gfortran. Ajout d un Makefile flexible a option pour choisir ifort ou gfortran.
File size: 11.0 KB
Line 
1      SUBROUTINE STRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
2*     .. Scalar Arguments ..
3      REAL ALPHA
4      INTEGER LDA,LDB,M,N
5      CHARACTER DIAG,SIDE,TRANSA,UPLO
6*     ..
7*     .. Array Arguments ..
8      REAL A(LDA,*),B(LDB,*)
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  STRMM  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**T.
22*
23*  Arguments
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**T.
53*
54*              TRANSA = 'C' or 'c'   op( A ) = A**T.
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  - REAL            .
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      - REAL             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      - REAL             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*  Further Details
118*  ===============
119*
120*  Level 3 Blas routine.
121*
122*  -- Written on 8-February-1989.
123*     Jack Dongarra, Argonne National Laboratory.
124*     Iain Duff, AERE Harwell.
125*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
126*     Sven Hammarling, Numerical Algorithms Group Ltd.
127*
128*  =====================================================================
129*
130*     .. External Functions ..
131      LOGICAL LSAME
132      EXTERNAL LSAME
133*     ..
134*     .. External Subroutines ..
135      EXTERNAL XERBLA
136*     ..
137*     .. Intrinsic Functions ..
138      INTRINSIC MAX
139*     ..
140*     .. Local Scalars ..
141      REAL TEMP
142      INTEGER I,INFO,J,K,NROWA
143      LOGICAL LSIDE,NOUNIT,UPPER
144*     ..
145*     .. Parameters ..
146      REAL ONE,ZERO
147      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
148*     ..
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. (.NOT.LSAME(SIDE,'R'))) THEN
163          INFO = 1
164      ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
165          INFO = 2
166      ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
167     +         (.NOT.LSAME(TRANSA,'T')) .AND.
168     +         (.NOT.LSAME(TRANSA,'C'))) THEN
169          INFO = 3
170      ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
171          INFO = 4
172      ELSE IF (M.LT.0) THEN
173          INFO = 5
174      ELSE IF (N.LT.0) THEN
175          INFO = 6
176      ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
177          INFO = 9
178      ELSE IF (LDB.LT.MAX(1,M)) THEN
179          INFO = 11
180      END IF
181      IF (INFO.NE.0) THEN
182          CALL XERBLA('STRMM ',INFO)
183          RETURN
184      END IF
185*
186*     Quick return if possible.
187*
188      IF (M.EQ.0 .OR. N.EQ.0) RETURN
189*
190*     And when  alpha.eq.zero.
191*
192      IF (ALPHA.EQ.ZERO) THEN
193          DO 20 J = 1,N
194              DO 10 I = 1,M
195                  B(I,J) = ZERO
196   10         CONTINUE
197   20     CONTINUE
198          RETURN
199      END IF
200*
201*     Start the operations.
202*
203      IF (LSIDE) THEN
204          IF (LSAME(TRANSA,'N')) THEN
205*
206*           Form  B := alpha*A*B.
207*
208              IF (UPPER) THEN
209                  DO 50 J = 1,N
210                      DO 40 K = 1,M
211                          IF (B(K,J).NE.ZERO) THEN
212                              TEMP = ALPHA*B(K,J)
213                              DO 30 I = 1,K - 1
214                                  B(I,J) = B(I,J) + TEMP*A(I,K)
215   30                         CONTINUE
216                              IF (NOUNIT) TEMP = TEMP*A(K,K)
217                              B(K,J) = TEMP
218                          END IF
219   40                 CONTINUE
220   50             CONTINUE
221              ELSE
222                  DO 80 J = 1,N
223                      DO 70 K = M,1,-1
224                          IF (B(K,J).NE.ZERO) THEN
225                              TEMP = ALPHA*B(K,J)
226                              B(K,J) = TEMP
227                              IF (NOUNIT) 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*A**T*B.
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) TEMP = TEMP*A(I,I)
244                          DO 90 K = 1,I - 1
245                              TEMP = TEMP + A(K,I)*B(K,J)
246   90                     CONTINUE
247                          B(I,J) = ALPHA*TEMP
248  100                 CONTINUE
249  110             CONTINUE
250              ELSE
251                  DO 140 J = 1,N
252                      DO 130 I = 1,M
253                          TEMP = B(I,J)
254                          IF (NOUNIT) TEMP = TEMP*A(I,I)
255                          DO 120 K = I + 1,M
256                              TEMP = TEMP + A(K,I)*B(K,J)
257  120                     CONTINUE
258                          B(I,J) = ALPHA*TEMP
259  130                 CONTINUE
260  140             CONTINUE
261              END IF
262          END IF
263      ELSE
264          IF (LSAME(TRANSA,'N')) THEN
265*
266*           Form  B := alpha*B*A.
267*
268              IF (UPPER) THEN
269                  DO 180 J = N,1,-1
270                      TEMP = ALPHA
271                      IF (NOUNIT) TEMP = TEMP*A(J,J)
272                      DO 150 I = 1,M
273                          B(I,J) = TEMP*B(I,J)
274  150                 CONTINUE
275                      DO 170 K = 1,J - 1
276                          IF (A(K,J).NE.ZERO) THEN
277                              TEMP = ALPHA*A(K,J)
278                              DO 160 I = 1,M
279                                  B(I,J) = B(I,J) + TEMP*B(I,K)
280  160                         CONTINUE
281                          END IF
282  170                 CONTINUE
283  180             CONTINUE
284              ELSE
285                  DO 220 J = 1,N
286                      TEMP = ALPHA
287                      IF (NOUNIT) TEMP = TEMP*A(J,J)
288                      DO 190 I = 1,M
289                          B(I,J) = TEMP*B(I,J)
290  190                 CONTINUE
291                      DO 210 K = J + 1,N
292                          IF (A(K,J).NE.ZERO) THEN
293                              TEMP = ALPHA*A(K,J)
294                              DO 200 I = 1,M
295                                  B(I,J) = B(I,J) + TEMP*B(I,K)
296  200                         CONTINUE
297                          END IF
298  210                 CONTINUE
299  220             CONTINUE
300              END IF
301          ELSE
302*
303*           Form  B := alpha*B*A**T.
304*
305              IF (UPPER) THEN
306                  DO 260 K = 1,N
307                      DO 240 J = 1,K - 1
308                          IF (A(J,K).NE.ZERO) THEN
309                              TEMP = ALPHA*A(J,K)
310                              DO 230 I = 1,M
311                                  B(I,J) = B(I,J) + TEMP*B(I,K)
312  230                         CONTINUE
313                          END IF
314  240                 CONTINUE
315                      TEMP = ALPHA
316                      IF (NOUNIT) TEMP = TEMP*A(K,K)
317                      IF (TEMP.NE.ONE) THEN
318                          DO 250 I = 1,M
319                              B(I,K) = TEMP*B(I,K)
320  250                     CONTINUE
321                      END IF
322  260             CONTINUE
323              ELSE
324                  DO 300 K = N,1,-1
325                      DO 280 J = K + 1,N
326                          IF (A(J,K).NE.ZERO) THEN
327                              TEMP = ALPHA*A(J,K)
328                              DO 270 I = 1,M
329                                  B(I,J) = B(I,J) + TEMP*B(I,K)
330  270                         CONTINUE
331                          END IF
332  280                 CONTINUE
333                      TEMP = ALPHA
334                      IF (NOUNIT) TEMP = TEMP*A(K,K)
335                      IF (TEMP.NE.ONE) THEN
336                          DO 290 I = 1,M
337                              B(I,K) = TEMP*B(I,K)
338  290                     CONTINUE
339                      END IF
340  300             CONTINUE
341              END IF
342          END IF
343      END IF
344*
345      RETURN
346*
347*     End of STRMM .
348*
349      END
Note: See TracBrowser for help on using the repository browser.