--- rpl/lapack/blas/zgemm.f 2010/08/13 21:03:41 1.5 +++ rpl/lapack/blas/zgemm.f 2015/11/26 11:44:15 1.12 @@ -1,133 +1,205 @@ - SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) -* .. Scalar Arguments .. - DOUBLE COMPLEX ALPHA,BETA - INTEGER K,LDA,LDB,LDC,M,N - CHARACTER TRANSA,TRANSB -* .. -* .. Array Arguments .. - DOUBLE COMPLEX A(LDA,*),B(LDB,*),C(LDC,*) -* .. +*> \brief \b ZGEMM * -* Purpose -* ======= +* =========== DOCUMENTATION =========== * -* ZGEMM performs one of the matrix-matrix operations +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ * -* C := alpha*op( A )*op( B ) + beta*C, +* Definition: +* =========== +* +* SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) +* +* .. Scalar Arguments .. +* COMPLEX*16 ALPHA,BETA +* INTEGER K,LDA,LDB,LDC,M,N +* CHARACTER TRANSA,TRANSB +* .. +* .. Array Arguments .. +* COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZGEMM 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**T or op( X ) = X**H, +*> +*> 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. +*> \endverbatim * -* where op( X ) is one of +* Arguments: +* ========== * -* op( X ) = X or op( X ) = X' or op( X ) = conjg( X' ), +*> \param[in] TRANSA +*> \verbatim +*> TRANSA is 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**T. +*> +*> TRANSA = 'C' or 'c', op( A ) = A**H. +*> \endverbatim +*> +*> \param[in] TRANSB +*> \verbatim +*> TRANSB is 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**T. +*> +*> TRANSB = 'C' or 'c', op( B ) = B**H. +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is 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. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is 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. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is 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. +*> \endverbatim +*> +*> \param[in] ALPHA +*> \verbatim +*> ALPHA is COMPLEX*16 +*> On entry, ALPHA specifies the scalar alpha. +*> \endverbatim +*> +*> \param[in] A +*> \verbatim +*> A is COMPLEX*16 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. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is 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 ). +*> \endverbatim +*> +*> \param[in] B +*> \verbatim +*> B is COMPLEX*16 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. +*> \endverbatim +*> +*> \param[in] LDB +*> \verbatim +*> LDB is 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 ). +*> \endverbatim +*> +*> \param[in] BETA +*> \verbatim +*> BETA is COMPLEX*16 +*> On entry, BETA specifies the scalar beta. When BETA is +*> supplied as zero then C need not be set on input. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is COMPLEX*16 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 ). +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is 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 ). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2015 +* +*> \ingroup complex16_blas_level3 +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> 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. +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC) * -* 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. +* -- Reference BLAS level3 routine (version 3.6.0) -- +* -- Reference BLAS is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2015 * -* Arguments -* ========== -* -* 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 ) = conjg( 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 ) = conjg( 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 - COMPLEX*16 . -* On entry, ALPHA specifies the scalar alpha. -* Unchanged on exit. -* -* A - COMPLEX*16 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 - COMPLEX*16 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 - COMPLEX*16 . -* 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 - COMPLEX*16 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. -* -* Further Details -* =============== -* -* 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. +* .. Scalar Arguments .. + COMPLEX*16 ALPHA,BETA + INTEGER K,LDA,LDB,LDC,M,N + CHARACTER TRANSA,TRANSB +* .. +* .. Array Arguments .. + COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*) +* .. * * ===================================================================== * @@ -142,14 +214,14 @@ INTRINSIC DCONJG,MAX * .. * .. Local Scalars .. - DOUBLE COMPLEX TEMP + COMPLEX*16 TEMP INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB LOGICAL CONJA,CONJB,NOTA,NOTB * .. * .. Parameters .. - DOUBLE COMPLEX ONE + COMPLEX*16 ONE PARAMETER (ONE= (1.0D+0,0.0D+0)) - DOUBLE COMPLEX ZERO + COMPLEX*16 ZERO PARAMETER (ZERO= (0.0D+0,0.0D+0)) * .. * @@ -245,17 +317,15 @@ 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 + TEMP = ALPHA*B(L,J) + DO 70 I = 1,M + C(I,J) = C(I,J) + TEMP*A(I,L) + 70 CONTINUE 80 CONTINUE 90 CONTINUE ELSE IF (CONJA) THEN * -* Form C := alpha*conjg( A' )*B + beta*C. +* Form C := alpha*A**H*B + beta*C. * DO 120 J = 1,N DO 110 I = 1,M @@ -272,7 +342,7 @@ 120 CONTINUE ELSE * -* Form C := alpha*A'*B + beta*C +* Form C := alpha*A**T*B + beta*C * DO 150 J = 1,N DO 140 I = 1,M @@ -291,7 +361,7 @@ ELSE IF (NOTA) THEN IF (CONJB) THEN * -* Form C := alpha*A*conjg( B' ) + beta*C. +* Form C := alpha*A*B**H + beta*C. * DO 200 J = 1,N IF (BETA.EQ.ZERO) THEN @@ -304,17 +374,15 @@ 170 CONTINUE END IF DO 190 L = 1,K - IF (B(J,L).NE.ZERO) THEN - TEMP = ALPHA*DCONJG(B(J,L)) - DO 180 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 180 CONTINUE - END IF + TEMP = ALPHA*DCONJG(B(J,L)) + DO 180 I = 1,M + C(I,J) = C(I,J) + TEMP*A(I,L) + 180 CONTINUE 190 CONTINUE 200 CONTINUE ELSE * -* Form C := alpha*A*B' + beta*C +* Form C := alpha*A*B**T + beta*C * DO 250 J = 1,N IF (BETA.EQ.ZERO) THEN @@ -327,19 +395,17 @@ 220 CONTINUE END IF DO 240 L = 1,K - IF (B(J,L).NE.ZERO) THEN - TEMP = ALPHA*B(J,L) - DO 230 I = 1,M - C(I,J) = C(I,J) + TEMP*A(I,L) - 230 CONTINUE - END IF + TEMP = ALPHA*B(J,L) + DO 230 I = 1,M + C(I,J) = C(I,J) + TEMP*A(I,L) + 230 CONTINUE 240 CONTINUE 250 CONTINUE END IF ELSE IF (CONJA) THEN IF (CONJB) THEN * -* Form C := alpha*conjg( A' )*conjg( B' ) + beta*C. +* Form C := alpha*A**H*B**H + beta*C. * DO 280 J = 1,N DO 270 I = 1,M @@ -356,7 +422,7 @@ 280 CONTINUE ELSE * -* Form C := alpha*conjg( A' )*B' + beta*C +* Form C := alpha*A**H*B**T + beta*C * DO 310 J = 1,N DO 300 I = 1,M @@ -375,7 +441,7 @@ ELSE IF (CONJB) THEN * -* Form C := alpha*A'*conjg( B' ) + beta*C +* Form C := alpha*A**T*B**H + beta*C * DO 340 J = 1,N DO 330 I = 1,M @@ -392,7 +458,7 @@ 340 CONTINUE ELSE * -* Form C := alpha*A'*B' + beta*C +* Form C := alpha*A**T*B**T + beta*C * DO 370 J = 1,N DO 360 I = 1,M