updated top-level README and version_decl for V4.5 (#1847)
[WRF.git] / var / external / blas / dgemm.inc
blob968459e476faf355872db0f5b0e3e255d2cae219
1       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
2 !     .. Scalar Arguments ..
3       DOUBLE PRECISION ALPHA,BETA
4       INTEGER K,LDA,LDB,LDC,M,N
5       CHARACTER TRANSA,TRANSB
6 !     ..
7 !     .. Array Arguments ..
8       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
9 !     ..
11 !  Purpose
12 !  =======
14 !  DGEMM  performs one of the matrix-matrix operations
16 !     C := alpha*op( A )*op( B ) + beta*C,
18 !  where  op( X ) is one of
20 !     op( X ) = X   or   op( X ) = X',
22 !  alpha and beta are scalars, and A, B and C are matrices, with op( A )
23 !  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
25 !  Arguments
26 !  ==========
28 !  TRANSA - CHARACTER*1.
29 !           On entry, TRANSA specifies the form of op( A ) to be used in
30 !           the matrix multiplication as follows:
32 !              TRANSA = 'N' or 'n',  op( A ) = A.
34 !              TRANSA = 'T' or 't',  op( A ) = A'.
36 !              TRANSA = 'C' or 'c',  op( A ) = A'.
38 !           Unchanged on exit.
40 !  TRANSB - CHARACTER*1.
41 !           On entry, TRANSB specifies the form of op( B ) to be used in
42 !           the matrix multiplication as follows:
44 !              TRANSB = 'N' or 'n',  op( B ) = B.
46 !              TRANSB = 'T' or 't',  op( B ) = B'.
48 !              TRANSB = 'C' or 'c',  op( B ) = B'.
50 !           Unchanged on exit.
52 !  M      - INTEGER.
53 !           On entry,  M  specifies  the number  of rows  of the  matrix
54 !           op( A )  and of the  matrix  C.  M  must  be at least  zero.
55 !           Unchanged on exit.
57 !  N      - INTEGER.
58 !           On entry,  N  specifies the number  of columns of the matrix
59 !           op( B ) and the number of columns of the matrix C. N must be
60 !           at least zero.
61 !           Unchanged on exit.
63 !  K      - INTEGER.
64 !           On entry,  K  specifies  the number of columns of the matrix
65 !           op( A ) and the number of rows of the matrix op( B ). K must
66 !           be at least  zero.
67 !           Unchanged on exit.
69 !  ALPHA  - DOUBLE PRECISION.
70 !           On entry, ALPHA specifies the scalar alpha.
71 !           Unchanged on exit.
73 !  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
74 !           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
75 !           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
76 !           part of the array  A  must contain the matrix  A,  otherwise
77 !           the leading  k by m  part of the array  A  must contain  the
78 !           matrix A.
79 !           Unchanged on exit.
81 !  LDA    - INTEGER.
82 !           On entry, LDA specifies the first dimension of A as declared
83 !           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
84 !           LDA must be at least  max( 1, m ), otherwise  LDA must be at
85 !           least  max( 1, k ).
86 !           Unchanged on exit.
88 !  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
89 !           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
90 !           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
91 !           part of the array  B  must contain the matrix  B,  otherwise
92 !           the leading  n by k  part of the array  B  must contain  the
93 !           matrix B.
94 !           Unchanged on exit.
96 !  LDB    - INTEGER.
97 !           On entry, LDB specifies the first dimension of B as declared
98 !           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
99 !           LDB must be at least  max( 1, k ), otherwise  LDB must be at
100 !           least  max( 1, n ).
101 !           Unchanged on exit.
103 !  BETA   - DOUBLE PRECISION.
104 !           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
105 !           supplied as zero then C need not be set on input.
106 !           Unchanged on exit.
108 !  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
109 !           Before entry, the leading  m by n  part of the array  C must
110 !           contain the matrix  C,  except when  beta  is zero, in which
111 !           case C need not be set on entry.
112 !           On exit, the array  C  is overwritten by the  m by n  matrix
113 !           ( alpha*op( A )*op( B ) + beta*C ).
115 !  LDC    - INTEGER.
116 !           On entry, LDC specifies the first dimension of C as declared
117 !           in  the  calling  (sub)  program.   LDC  must  be  at  least
118 !           max( 1, m ).
119 !           Unchanged on exit.
122 !  Level 3 Blas routine.
124 !  -- Written on 8-February-1989.
125 !     Jack Dongarra, Argonne National Laboratory.
126 !     Iain Duff, AERE Harwell.
127 !     Jeremy Du Croz, Numerical Algorithms Group Ltd.
128 !     Sven Hammarling, Numerical Algorithms Group Ltd.
131 !     .. External Functions ..
132 !     LOGICAL LSAME
133 !     EXTERNAL LSAME
134 !     ..
135 !     .. External Subroutines ..
136 !     EXTERNAL XERBLA
137 !     ..
138 !     .. Intrinsic Functions ..
139       INTRINSIC MAX
140 !     ..
141 !     .. Local Scalars ..
142       DOUBLE PRECISION TEMP
143       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
144       LOGICAL NOTA,NOTB
145 !     ..
146 !     .. Parameters ..
147       DOUBLE PRECISION ONE,ZERO
148       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
149 !     ..
151 !     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
152 !     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
153 !     and  columns of  A  and the  number of  rows  of  B  respectively.
155       NOTA = LSAME(TRANSA,'N')
156       NOTB = LSAME(TRANSB,'N')
157       IF (NOTA) THEN
158           NROWA = M
159           NCOLA = K
160       ELSE
161           NROWA = K
162           NCOLA = M
163       END IF
164       IF (NOTB) THEN
165           NROWB = K
166       ELSE
167           NROWB = N
168       END IF
170 !     Test the input parameters.
172       INFO = 0
173       IF ((.NOT.NOTA) .AND. (.NOT.LSAME(TRANSA,'C')) .AND. &
174           (.NOT.LSAME(TRANSA,'T'))) THEN
175           INFO = 1
176       ELSE IF ((.NOT.NOTB) .AND. (.NOT.LSAME(TRANSB,'C')) .AND. &
177                (.NOT.LSAME(TRANSB,'T'))) THEN
178           INFO = 2
179       ELSE IF (M.LT.0) THEN
180           INFO = 3
181       ELSE IF (N.LT.0) THEN
182           INFO = 4
183       ELSE IF (K.LT.0) THEN
184           INFO = 5
185       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
186           INFO = 8
187       ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
188           INFO = 10
189       ELSE IF (LDC.LT.MAX(1,M)) THEN
190           INFO = 13
191       END IF
192       IF (INFO.NE.0) THEN
193           CALL XERBLA('DGEMM ',INFO)
194           RETURN
195       END IF
197 !     Quick return if possible.
199       IF ((M.EQ.0) .OR. (N.EQ.0) .OR. &
200           (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
202 !     And if  alpha.eq.zero.
204       IF (ALPHA.EQ.ZERO) THEN
205           IF (BETA.EQ.ZERO) THEN
206               DO 20 J = 1,N
207                   DO 10 I = 1,M
208                       C(I,J) = ZERO
209    10             CONTINUE
210    20         CONTINUE
211           ELSE
212               DO 40 J = 1,N
213                   DO 30 I = 1,M
214                       C(I,J) = BETA*C(I,J)
215    30             CONTINUE
216    40         CONTINUE
217           END IF
218           RETURN
219       END IF
221 !     Start the operations.
223       IF (NOTB) THEN
224           IF (NOTA) THEN
226 !           Form  C := alpha*A*B + beta*C.
228               DO 90 J = 1,N
229                   IF (BETA.EQ.ZERO) THEN
230                       DO 50 I = 1,M
231                           C(I,J) = ZERO
232    50                 CONTINUE
233                   ELSE IF (BETA.NE.ONE) THEN
234                       DO 60 I = 1,M
235                           C(I,J) = BETA*C(I,J)
236    60                 CONTINUE
237                   END IF
238                   DO 80 L = 1,K
239                       IF (B(L,J).NE.ZERO) THEN
240                           TEMP = ALPHA*B(L,J)
241                           DO 70 I = 1,M
242                               C(I,J) = C(I,J) + TEMP*A(I,L)
243    70                     CONTINUE
244                       END IF
245    80             CONTINUE
246    90         CONTINUE
247           ELSE
249 !           Form  C := alpha*A'*B + beta*C
251               DO 120 J = 1,N
252                   DO 110 I = 1,M
253                       TEMP = ZERO
254                       DO 100 L = 1,K
255                           TEMP = TEMP + A(L,I)*B(L,J)
256   100                 CONTINUE
257                       IF (BETA.EQ.ZERO) THEN
258                           C(I,J) = ALPHA*TEMP
259                       ELSE
260                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
261                       END IF
262   110             CONTINUE
263   120         CONTINUE
264           END IF
265       ELSE
266           IF (NOTA) THEN
268 !           Form  C := alpha*A*B' + beta*C
270               DO 170 J = 1,N
271                   IF (BETA.EQ.ZERO) THEN
272                       DO 130 I = 1,M
273                           C(I,J) = ZERO
274   130                 CONTINUE
275                   ELSE IF (BETA.NE.ONE) THEN
276                       DO 140 I = 1,M
277                           C(I,J) = BETA*C(I,J)
278   140                 CONTINUE
279                   END IF
280                   DO 160 L = 1,K
281                       IF (B(J,L).NE.ZERO) THEN
282                           TEMP = ALPHA*B(J,L)
283                           DO 150 I = 1,M
284                               C(I,J) = C(I,J) + TEMP*A(I,L)
285   150                     CONTINUE
286                       END IF
287   160             CONTINUE
288   170         CONTINUE
289           ELSE
291 !           Form  C := alpha*A'*B' + beta*C
293               DO 200 J = 1,N
294                   DO 190 I = 1,M
295                       TEMP = ZERO
296                       DO 180 L = 1,K
297                           TEMP = TEMP + A(L,I)*B(J,L)
298   180                 CONTINUE
299                       IF (BETA.EQ.ZERO) THEN
300                           C(I,J) = ALPHA*TEMP
301                       ELSE
302                           C(I,J) = ALPHA*TEMP + BETA*C(I,J)
303                       END IF
304   190             CONTINUE
305   200         CONTINUE
306           END IF
307       END IF
309       RETURN
311 !     End of DGEMM .
313       END SUBROUTINE DGEMM