trunk/blas/dgemm.f

      SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
     $                   BETA, C, LDC )
*     .. Scalar Arguments ..
      CHARACTER*1        TRANSA, TRANSB
      INTEGER            M, N, K, LDA, LDB, LDC
      DOUBLE PRECISION   ALPHA, BETA
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
*     ..
*
*  Purpose
*  =======
*
*  DGEMM  performs one of the matrix-matrix operations
*
*     C := alpha*op( A )*op( B ) + beta*C,
*
*  where  op( X ) is one of
*
*     op( X ) = X   or   op( X ) = X',
*
*  alpha and beta are scalars, and A, B and C are matrices, with op( A )
*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
*
*  Parameters
*  ==========
*
*  TRANSA - CHARACTER*1.
*           On entry, TRANSA specifies the form of op( A ) to be used in
*           the matrix multiplication as follows:
*
*              TRANSA = 'N' or 'n',  op( A ) = A.
*
*              TRANSA = 'T' or 't',  op( A ) = A'.
*
*              TRANSA = 'C' or 'c',  op( A ) = A'.
*
*           Unchanged on exit.
*
*  TRANSB - CHARACTER*1.
*           On entry, TRANSB specifies the form of op( B ) to be used in
*           the matrix multiplication as follows:
*
*              TRANSB = 'N' or 'n',  op( B ) = B.
*
*              TRANSB = 'T' or 't',  op( B ) = B'.
*
*              TRANSB = 'C' or 'c',  op( B ) = B'.
*
*           Unchanged on exit.
*
*  M      - INTEGER.
*           On entry,  M  specifies  the number  of rows  of the  matrix
*           op( A )  and of the  matrix  C.  M  must  be at least  zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry,  N  specifies the number  of columns of the matrix
*           op( B ) and the number of columns of the matrix C. N must be
*           at least zero.
*           Unchanged on exit.
*
*  K      - INTEGER.
*           On entry,  K  specifies  the number of columns of the matrix
*           op( A ) and the number of rows of the matrix op( B ). K must
*           be at least  zero.
*           Unchanged on exit.
*
*  ALPHA  - DOUBLE PRECISION.
*           On entry, ALPHA specifies the scalar alpha.
*           Unchanged on exit.
*
*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
*           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
*           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
*           part of the array  A  must contain the matrix  A,  otherwise
*           the leading  k by m  part of the array  A  must contain  the
*           matrix A.
*           Unchanged on exit.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
*           LDA must be at least  max( 1, m ), otherwise  LDA must be at
*           least  max( 1, k ).
*           Unchanged on exit.
*
*  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
*           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
*           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
*           part of the array  B  must contain the matrix  B,  otherwise
*           the leading  n by k  part of the array  B  must contain  the
*           matrix B.
*           Unchanged on exit.
*
*  LDB    - INTEGER.
*           On entry, LDB specifies the first dimension of B as declared
*           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
*           LDB must be at least  max( 1, k ), otherwise  LDB must be at
*           least  max( 1, n ).
*           Unchanged on exit.
*
*  BETA   - DOUBLE PRECISION.
*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
*           supplied as zero then C need not be set on input.
*           Unchanged on exit.
*
*  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
*           Before entry, the leading  m by n  part of the array  C must
*           contain the matrix  C,  except when  beta  is zero, in which
*           case C need not be set on entry.
*           On exit, the array  C  is overwritten by the  m by n  matrix
*           ( alpha*op( A )*op( B ) + beta*C ).
*
*  LDC    - INTEGER.
*           On entry, LDC specifies the first dimension of C as declared
*           in  the  calling  (sub)  program.   LDC  must  be  at  least
*           max( 1, m ).
*           Unchanged on exit.
*
*
*  Level 3 Blas routine.
*
*  -- Written on 8-February-1989.
*     Jack Dongarra, Argonne National Laboratory.
*     Iain Duff, AERE Harwell.
*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
*     Sven Hammarling, Numerical Algorithms Group Ltd.
*
*
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     .. External Subroutines ..
      EXTERNAL           XERBLA
*     .. Intrinsic Functions ..
      INTRINSIC          MAX
*     .. Local Scalars ..
      LOGICAL            NOTA, NOTB
      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
      DOUBLE PRECISION   TEMP
*     .. Parameters ..
      DOUBLE PRECISION   ONE         , ZERO
      PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Executable Statements ..
*
*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
*     transposed and set  NROWA, NCOLA and  NROWB  as the number of rows
*     and  columns of  A  and the  number of  rows  of  B  respectively.
*
      NOTA  = LSAME( TRANSA, 'N' )
      NOTB  = LSAME( TRANSB, 'N' )
      IF( NOTA )THEN
     NROWA = M
     NCOLA = K
      ELSE
     NROWA = K
     NCOLA = M
      END IF
      IF( NOTB )THEN
     NROWB = K
      ELSE
     NROWB = N
      END IF
*
*     Test the input parameters.
*
      INFO = 0
      IF(      ( .NOT.NOTA                 ).AND.
     $         ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
     INFO = 1
      ELSE IF( ( .NOT.NOTB                 ).AND.
     $         ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
     INFO = 2
      ELSE IF( M  .LT.0               )THEN
     INFO = 3
      ELSE IF( N  .LT.0               )THEN
     INFO = 4
      ELSE IF( K  .LT.0               )THEN
     INFO = 5
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
     INFO = 8
      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
     INFO = 10
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
     INFO = 13
      END IF
      IF( INFO.NE.0 )THEN
     CALL XERBLA( 'DGEMM ', INFO )
     RETURN
      END IF
*
*     Quick return if possible.
*
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
     $   RETURN
*
*     And if  alpha.eq.zero.
*
      IF( ALPHA.EQ.ZERO )THEN
     IF( BETA.EQ.ZERO )THEN
        DO 20, J = 1, N
           DO 10, I = 1, M
          C( I, J ) = ZERO
   10          CONTINUE
   20       CONTINUE
     ELSE
        DO 40, J = 1, N
           DO 30, I = 1, M
          C( I, J ) = BETA*C( I, J )
   30          CONTINUE
   40       CONTINUE
     END IF
     RETURN
      END IF
*
*     Start the operations.
*
      IF( NOTB )THEN
     IF( NOTA )THEN
*
*           Form  C := alpha*A*B + beta*C.
*
        DO 90, J = 1, N
           IF( BETA.EQ.ZERO )THEN
          DO 50, I = 1, M
             C( I, J ) = ZERO
   50             CONTINUE
           ELSE IF( BETA.NE.ONE )THEN
          DO 60, I = 1, M
             C( I, J ) = BETA*C( I, J )
   60             CONTINUE
           END IF
           DO 80, L = 1, K
          IF( B( L, J ).NE.ZERO )THEN
             TEMP = ALPHA*B( L, J )
             DO 70, I = 1, M
            C( I, J ) = C( I, J ) + TEMP*A( I, L )
   70                CONTINUE
          END IF
   80          CONTINUE
   90       CONTINUE
     ELSE
*
*           Form  C := alpha*A'*B + beta*C
*
        DO 120, J = 1, N
           DO 110, I = 1, M
          TEMP = ZERO
          DO 100, L = 1, K
             TEMP = TEMP + A( L, I )*B( L, J )
  100             CONTINUE
          IF( BETA.EQ.ZERO )THEN
             C( I, J ) = ALPHA*TEMP
          ELSE
             C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
          END IF
  110          CONTINUE
  120       CONTINUE
     END IF
      ELSE
     IF( NOTA )THEN
*
*           Form  C := alpha*A*B' + beta*C
*
        DO 170, J = 1, N
           IF( BETA.EQ.ZERO )THEN
          DO 130, I = 1, M
             C( I, J ) = ZERO
  130             CONTINUE
           ELSE IF( BETA.NE.ONE )THEN
          DO 140, I = 1, M
             C( I, J ) = BETA*C( I, J )
  140             CONTINUE
           END IF
           DO 160, L = 1, K
          IF( B( J, L ).NE.ZERO )THEN
             TEMP = ALPHA*B( J, L )
             DO 150, I = 1, M
            C( I, J ) = C( I, J ) + TEMP*A( I, L )
  150                CONTINUE
          END IF
  160          CONTINUE
  170       CONTINUE
     ELSE
*
*           Form  C := alpha*A'*B' + beta*C
*
        DO 200, J = 1, N
           DO 190, I = 1, M
          TEMP = ZERO
          DO 180, L = 1, K
             TEMP = TEMP + A( L, I )*B( J, L )
  180             CONTINUE
          IF( BETA.EQ.ZERO )THEN
             C( I, J ) = ALPHA*TEMP
          ELSE
             C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
          END IF
  190          CONTINUE
  200       CONTINUE
     END IF
      END IF
*
      RETURN
*
*     End of DGEMM .
*
      END

1	aw0a	1	SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
2			$ BETA, C, LDC )
3			* .. Scalar Arguments ..
4			CHARACTER*1 TRANSA, TRANSB
5			INTEGER M, N, K, LDA, LDB, LDC
6			DOUBLE PRECISION ALPHA, BETA
7			* .. Array Arguments ..
8			DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
9			* ..
10			*
11			* Purpose
12			* =======
13			*
14			* DGEMM performs one of the matrix-matrix operations
15			*
16			* C := alphaop( A )op( B ) + beta*C,
17			*
18			* where op( X ) is one of
19			*
20			* op( X ) = X or op( X ) = X',
21			*
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.
24			*
25			* Parameters
26			* ==========
27			*
28			* TRANSA - CHARACTER*1.
29			* On entry, TRANSA specifies the form of op( A ) to be used in
30			* the matrix multiplication as follows:
31			*
32			* TRANSA = 'N' or 'n', op( A ) = A.
33			*
34			* TRANSA = 'T' or 't', op( A ) = A'.
35			*
36			* TRANSA = 'C' or 'c', op( A ) = A'.
37			*
38			* Unchanged on exit.
39			*
40			* TRANSB - CHARACTER*1.
41			* On entry, TRANSB specifies the form of op( B ) to be used in
42			* the matrix multiplication as follows:
43			*
44			* TRANSB = 'N' or 'n', op( B ) = B.
45			*
46			* TRANSB = 'T' or 't', op( B ) = B'.
47			*
48			* TRANSB = 'C' or 'c', op( B ) = B'.
49			*
50			* Unchanged on exit.
51			*
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.
56			*
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.
62			*
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.
68			*
69			* ALPHA - DOUBLE PRECISION.
70			* On entry, ALPHA specifies the scalar alpha.
71			* Unchanged on exit.
72			*
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.
80			*
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.
87			*
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.
95			*
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.
102			*
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.
107			*
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			* ( alphaop( A )op( B ) + beta*C ).
114			*
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.
120			*
121			*
122			* Level 3 Blas routine.
123			*
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.
129			*
130			*
131			* .. External Functions ..
132			LOGICAL LSAME
133			EXTERNAL LSAME
134			* .. External Subroutines ..
135			EXTERNAL XERBLA
136			* .. Intrinsic Functions ..
137			INTRINSIC MAX
138			* .. Local Scalars ..
139			LOGICAL NOTA, NOTB
140			INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB
141			DOUBLE PRECISION TEMP
142			* .. Parameters ..
143			DOUBLE PRECISION ONE , ZERO
144			PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
145			* ..
146			* .. Executable Statements ..
147			*
148			* Set NOTA and NOTB as true if A and B respectively are not
149			* transposed and set NROWA, NCOLA and NROWB as the number of rows
150			* and columns of A and the number of rows of B respectively.
151			*
152			NOTA = LSAME( TRANSA, 'N' )
153			NOTB = LSAME( TRANSB, 'N' )
154			IF( NOTA )THEN
155			NROWA = M
156			NCOLA = K
157			ELSE
158			NROWA = K
159			NCOLA = M
160			END IF
161			IF( NOTB )THEN
162			NROWB = K
163			ELSE
164			NROWB = N
165			END IF
166			*
167			* Test the input parameters.
168			*
169			INFO = 0
170			IF( ( .NOT.NOTA ).AND.
171			$ ( .NOT.LSAME( TRANSA, 'C' ) ).AND.
172			$ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN
173			INFO = 1
174			ELSE IF( ( .NOT.NOTB ).AND.
175			$ ( .NOT.LSAME( TRANSB, 'C' ) ).AND.
176			$ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN
177			INFO = 2
178			ELSE IF( M .LT.0 )THEN
179			INFO = 3
180			ELSE IF( N .LT.0 )THEN
181			INFO = 4
182			ELSE IF( K .LT.0 )THEN
183			INFO = 5
184			ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
185			INFO = 8
186			ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
187			INFO = 10
188			ELSE IF( LDC.LT.MAX( 1, M ) )THEN
189			INFO = 13
190			END IF
191			IF( INFO.NE.0 )THEN
192			CALL XERBLA( 'DGEMM ', INFO )
193			RETURN
194			END IF
195			*
196			* Quick return if possible.
197			*
198			IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
199			$ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
200			$ RETURN
201			*
202			* And if alpha.eq.zero.
203			*
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
220			*
221			* Start the operations.
222			*
223			IF( NOTB )THEN
224			IF( NOTA )THEN
225			*
226			* Form C := alphaAB + beta*C.
227			*
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
248			*
249			* Form C := alphaA'B + beta*C
250			*
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 ) = ALPHATEMP + BETAC( I, J )
261			END IF
262			110 CONTINUE
263			120 CONTINUE
264			END IF
265			ELSE
266			IF( NOTA )THEN
267			*
268			* Form C := alphaAB' + beta*C
269			*
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
290			*
291			* Form C := alphaA'B' + beta*C
292			*
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 ) = ALPHATEMP + BETAC( I, J )
303			END IF
304			190 CONTINUE
305			200 CONTINUE
306			END IF
307			END IF
308			*
309			RETURN
310			*
311			* End of DGEMM .
312			*
313			END
314