Télécharger dgemm.eso

Retour à la liste

Numérotation des lignes :

dgemm
  1. C DGEMM SOURCE BP208322 22/09/16 21:15:02 11454
  2. *> \brief \b DGEMM
  3. *
  4. * =========== DOCUMENTATION ===========
  5. *
  6. * Online html documentation available at
  7. * http://www.netlib.org/lapack/explore-html/
  8. *
  9. * Definition:
  10. * ===========
  11. *
  12. * SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  13. *
  14. * .. Scalar Arguments ..
  15. * REAL ALPHA,BETA
  16. * INTEGER K,LDA,LDB,LDC,M,N
  17. * CHARACTER TRANSA,TRANSB
  18. * ..
  19. * .. Array Arguments ..
  20. * REAL A(LDA,*),B(LDB,*),C(LDC,*)
  21. * ..
  22. *
  23. *
  24. *> \par Purpose:
  25. * =============
  26. *>
  27. *> \verbatim
  28. *>
  29. *> DGEMM performs one of the matrix-matrix operations
  30. *>
  31. *> C := alpha*op( A )*op( B ) + beta*C,
  32. *>
  33. *> where op( X ) is one of
  34. *>
  35. *> op( X ) = X or op( X ) = X**T,
  36. *>
  37. *> alpha and beta are scalars, and A, B and C are matrices, with op( A )
  38. *> an m by k matrix, op( B ) a k by n matrix and C an m by n matrix.
  39. *> \endverbatim
  40. *
  41. * Arguments:
  42. * ==========
  43. *
  44. *> \param[in] TRANSA
  45. *> \verbatim
  46. *> TRANSA is 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. *> \endverbatim
  56. *>
  57. *> \param[in] TRANSB
  58. *> \verbatim
  59. *> TRANSB is CHARACTER*1
  60. *> On entry, TRANSB specifies the form of op( B ) to be used in
  61. *> the matrix multiplication as follows:
  62. *>
  63. *> TRANSB = 'N' or 'n', op( B ) = B.
  64. *>
  65. *> TRANSB = 'T' or 't', op( B ) = B**T.
  66. *>
  67. *> TRANSB = 'C' or 'c', op( B ) = B**T.
  68. *> \endverbatim
  69. *>
  70. *> \param[in] M
  71. *> \verbatim
  72. *> M is INTEGER
  73. *> On entry, M specifies the number of rows of the matrix
  74. *> op( A ) and of the matrix C. M must be at least zero.
  75. *> \endverbatim
  76. *>
  77. *> \param[in] N
  78. *> \verbatim
  79. *> N is INTEGER
  80. *> On entry, N specifies the number of columns of the matrix
  81. *> op( B ) and the number of columns of the matrix C. N must be
  82. *> at least zero.
  83. *> \endverbatim
  84. *>
  85. *> \param[in] K
  86. *> \verbatim
  87. *> K is INTEGER
  88. *> On entry, K specifies the number of columns of the matrix
  89. *> op( A ) and the number of rows of the matrix op( B ). K must
  90. *> be at least zero.
  91. *> \endverbatim
  92. *>
  93. *> \param[in] ALPHA
  94. *> \verbatim
  95. *> ALPHA is REAL
  96. *> On entry, ALPHA specifies the scalar alpha.
  97. *> \endverbatim
  98. *>
  99. *> \param[in] A
  100. *> \verbatim
  101. *> A is REAL array, dimension ( LDA, ka ), where ka is
  102. *> k when TRANSA = 'N' or 'n', and is m otherwise.
  103. *> Before entry with TRANSA = 'N' or 'n', the leading m by k
  104. *> part of the array A must contain the matrix A, otherwise
  105. *> the leading k by m part of the array A must contain the
  106. *> matrix A.
  107. *> \endverbatim
  108. *>
  109. *> \param[in] LDA
  110. *> \verbatim
  111. *> LDA is INTEGER
  112. *> On entry, LDA specifies the first dimension of A as declared
  113. *> in the calling (sub) program. When TRANSA = 'N' or 'n' then
  114. *> LDA must be at least max( 1, m ), otherwise LDA must be at
  115. *> least max( 1, k ).
  116. *> \endverbatim
  117. *>
  118. *> \param[in] B
  119. *> \verbatim
  120. *> B is REAL array, dimension ( LDB, kb ), where kb is
  121. *> n when TRANSB = 'N' or 'n', and is k otherwise.
  122. *> Before entry with TRANSB = 'N' or 'n', the leading k by n
  123. *> part of the array B must contain the matrix B, otherwise
  124. *> the leading n by k part of the array B must contain the
  125. *> matrix B.
  126. *> \endverbatim
  127. *>
  128. *> \param[in] LDB
  129. *> \verbatim
  130. *> LDB is INTEGER
  131. *> On entry, LDB specifies the first dimension of B as declared
  132. *> in the calling (sub) program. When TRANSB = 'N' or 'n' then
  133. *> LDB must be at least max( 1, k ), otherwise LDB must be at
  134. *> least max( 1, n ).
  135. *> \endverbatim
  136. *>
  137. *> \param[in] BETA
  138. *> \verbatim
  139. *> BETA is REAL
  140. *> On entry, BETA specifies the scalar beta. When BETA is
  141. *> supplied as zero then C need not be set on input.
  142. *> \endverbatim
  143. *>
  144. *> \param[in,out] C
  145. *> \verbatim
  146. *> C is REAL array, dimension ( LDC, N )
  147. *> Before entry, the leading m by n part of the array C must
  148. *> contain the matrix C, except when beta is zero, in which
  149. *> case C need not be set on entry.
  150. *> On exit, the array C is overwritten by the m by n matrix
  151. *> ( alpha*op( A )*op( B ) + beta*C ).
  152. *> \endverbatim
  153. *>
  154. *> \param[in] LDC
  155. *> \verbatim
  156. *> LDC is INTEGER
  157. *> On entry, LDC specifies the first dimension of C as declared
  158. *> in the calling (sub) program. LDC must be at least
  159. *> max( 1, m ).
  160. *> \endverbatim
  161. *
  162. * Authors:
  163. * ========
  164. *
  165. *> \author Univ. of Tennessee
  166. *> \author Univ. of California Berkeley
  167. *> \author Univ. of Colorado Denver
  168. *> \author NAG Ltd.
  169. *
  170. *> \ingroup double_blas_level3
  171. *
  172. *> \par Further Details:
  173. * =====================
  174. *>
  175. *> \verbatim
  176. *>
  177. *> Level 3 Blas routine.
  178. *>
  179. *> -- Written on 8-February-1989.
  180. *> Jack Dongarra, Argonne National Laboratory.
  181. *> Iain Duff, AERE Harwell.
  182. *> Jeremy Du Croz, Numerical Algorithms Group Ltd.
  183. *> Sven Hammarling, Numerical Algorithms Group Ltd.
  184. *> \endverbatim
  185. *>
  186. * =====================================================================
  187. SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
  188. *
  189. * -- Reference BLAS level3 routine --
  190. * -- Reference BLAS is a software package provided by Univ. of Tennessee, --
  191. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  192.  
  193. IMPLICIT INTEGER(I-N)
  194. IMPLICIT REAL*8(A-H,O-Z)
  195. *
  196. * .. Scalar Arguments ..
  197. REAL*8 ALPHA,BETA
  198. INTEGER K,LDA,LDB,LDC,M,N
  199. CHARACTER TRANSA,TRANSB
  200. * ..
  201. * .. Array Arguments ..
  202. REAL*8 A(LDA,*),B(LDB,*),C(LDC,*)
  203. * ..
  204. *
  205. * =====================================================================
  206. * ..
  207. * .. External Subroutines ..
  208. * EXTERNAL XERBLA
  209. * .. Parameters ..
  210. REAL*8 ONE, ZERO
  211. PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  212. * .. Local Scalars ..
  213. REAL*8 TEMP
  214. INTEGER I,INFO,J,L,NROWA,NROWB
  215. LOGICAL NOTA,NOTB
  216. *
  217. * Set NOTA and NOTB as true if A and B respectively are not
  218. * transposed and set NROWA and NROWB as the number of rows of A
  219. * and B respectively.
  220. *
  221. NOTA = (TRANSA.EQ.'N')
  222. NOTB = (TRANSB.EQ.'N')
  223. IF (NOTA) THEN
  224. NROWA = M
  225. ELSE
  226. NROWA = K
  227. END IF
  228. IF (NOTB) THEN
  229. NROWB = K
  230. ELSE
  231. NROWB = N
  232. END IF
  233. *
  234. * Test the input parameters.
  235. *
  236. INFO = 0
  237. IF ((.NOT.NOTA) .AND. (.NOT.(TRANSA.EQ.'C')) .AND.
  238. + (.NOT.(TRANSA.EQ.'T'))) THEN
  239. INFO = 1
  240. ELSE IF ((.NOT.NOTB) .AND. (.NOT.(TRANSB.EQ.'C')) .AND.
  241. + (.NOT.(TRANSB.EQ.'T'))) THEN
  242. INFO = 2
  243. ELSE IF (M.LT.0) THEN
  244. INFO = 3
  245. ELSE IF (N.LT.0) THEN
  246. INFO = 4
  247. ELSE IF (K.LT.0) THEN
  248. INFO = 5
  249. ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
  250. INFO = 8
  251. ELSE IF (LDB.LT.MAX(1,NROWB)) THEN
  252. INFO = 10
  253. ELSE IF (LDC.LT.MAX(1,M)) THEN
  254. INFO = 13
  255. END IF
  256. IF (INFO.NE.0) THEN
  257. CALL XERBLA('DGEMM ',INFO)
  258. RETURN
  259. END IF
  260. *
  261. * Quick return if possible.
  262. *
  263. IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
  264. + (((ALPHA.EQ.ZERO).OR. (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
  265. *
  266. * And if alpha.eq.zero.
  267. *
  268. IF (ALPHA.EQ.ZERO) THEN
  269. IF (BETA.EQ.ZERO) THEN
  270. DO 20 J = 1,N
  271. DO 10 I = 1,M
  272. C(I,J) = ZERO
  273. 10 CONTINUE
  274. 20 CONTINUE
  275. ELSE
  276. DO 40 J = 1,N
  277. DO 30 I = 1,M
  278. C(I,J) = BETA*C(I,J)
  279. 30 CONTINUE
  280. 40 CONTINUE
  281. END IF
  282. RETURN
  283. END IF
  284. *
  285. * Start the operations.
  286. *
  287. IF (NOTB) THEN
  288. IF (NOTA) THEN
  289. *
  290. * Form C := alpha*A*B + beta*C.
  291. *
  292. DO 90 J = 1,N
  293. IF (BETA.EQ.ZERO) THEN
  294. DO 50 I = 1,M
  295. C(I,J) = ZERO
  296. 50 CONTINUE
  297. ELSE IF (BETA.NE.ONE) THEN
  298. DO 60 I = 1,M
  299. C(I,J) = BETA*C(I,J)
  300. 60 CONTINUE
  301. END IF
  302. DO 80 L = 1,K
  303. TEMP = ALPHA*B(L,J)
  304. DO 70 I = 1,M
  305. C(I,J) = C(I,J) + TEMP*A(I,L)
  306. 70 CONTINUE
  307. 80 CONTINUE
  308. 90 CONTINUE
  309. ELSE
  310. *
  311. * Form C := alpha*A**T*B + beta*C
  312. *
  313. DO 120 J = 1,N
  314. DO 110 I = 1,M
  315. TEMP = ZERO
  316. DO 100 L = 1,K
  317. TEMP = TEMP + A(L,I)*B(L,J)
  318. 100 CONTINUE
  319. IF (BETA.EQ.ZERO) THEN
  320. C(I,J) = ALPHA*TEMP
  321. ELSE
  322. C(I,J) = ALPHA*TEMP + BETA*C(I,J)
  323. END IF
  324. 110 CONTINUE
  325. 120 CONTINUE
  326. END IF
  327. ELSE
  328. IF (NOTA) THEN
  329. *
  330. * Form C := alpha*A*B**T + beta*C
  331. *
  332. DO 170 J = 1,N
  333. IF (BETA.EQ.ZERO) THEN
  334. DO 130 I = 1,M
  335. C(I,J) = ZERO
  336. 130 CONTINUE
  337. ELSE IF (BETA.NE.ONE) THEN
  338. DO 140 I = 1,M
  339. C(I,J) = BETA*C(I,J)
  340. 140 CONTINUE
  341. END IF
  342. DO 160 L = 1,K
  343. TEMP = ALPHA*B(J,L)
  344. DO 150 I = 1,M
  345. C(I,J) = C(I,J) + TEMP*A(I,L)
  346. 150 CONTINUE
  347. 160 CONTINUE
  348. 170 CONTINUE
  349. ELSE
  350. *
  351. * Form C := alpha*A**T*B**T + beta*C
  352. *
  353. DO 200 J = 1,N
  354. DO 190 I = 1,M
  355. TEMP = ZERO
  356. DO 180 L = 1,K
  357. TEMP = TEMP + A(L,I)*B(J,L)
  358. 180 CONTINUE
  359. IF (BETA.EQ.ZERO) THEN
  360. C(I,J) = ALPHA*TEMP
  361. ELSE
  362. C(I,J) = ALPHA*TEMP + BETA*C(I,J)
  363. END IF
  364. 190 CONTINUE
  365. 200 CONTINUE
  366. END IF
  367. END IF
  368. *
  369. RETURN
  370. *
  371. * End of DGEMM
  372. *
  373. END
  374.  
  375.  
  376.  
  377.  
  378.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales