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

Annotation of /trunk/blas/dgemm.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: 9026 byte(s)
Setting up web subdirectory in repository
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 := alpha*op( 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     * ( alpha*op( 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 := alpha*A*B + 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 := alpha*A'*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 ) = ALPHA*TEMP + BETA*C( I, J )
261     END IF
262     110 CONTINUE
263     120 CONTINUE
264     END IF
265     ELSE
266     IF( NOTA )THEN
267     *
268     * Form C := alpha*A*B' + 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 := alpha*A'*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 ) = ALPHA*TEMP + BETA*C( 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    

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22