Diff for /rpl/lapack/blas/zgemm.f between versions 1.6 and 1.17

version 1.6, 2010/12/21 13:51:26 version 1.17, 2023/08/07 08:38:45
Line 1 Line 1
       SUBROUTINE ZGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)  *> \brief \b ZGEMM
 *     .. 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,*)  
 *     ..  
 *  *
 *  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, 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, 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, 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.
   *
   *> \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 )  *  -- Reference BLAS level3 routine --
 *  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.  *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
   *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *  *
 *  Arguments  *     .. Scalar Arguments ..
 *  ==========        COMPLEX*16 ALPHA,BETA
 *        INTEGER K,LDA,LDB,LDC,M,N
 *  TRANSA - CHARACTER*1.        CHARACTER TRANSA,TRANSB
 *           On entry, TRANSA specifies the form of op( A ) to be used in  *     ..
 *           the matrix multiplication as follows:  *     .. Array Arguments ..
 *        COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
 *              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.  
 *  *
 *  =====================================================================  *  =====================================================================
 *  *
Line 142 Line 211
       INTRINSIC DCONJG,MAX        INTRINSIC DCONJG,MAX
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       DOUBLE COMPLEX TEMP        COMPLEX*16 TEMP
       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB        INTEGER I,INFO,J,L,NROWA,NROWB
       LOGICAL CONJA,CONJB,NOTA,NOTB        LOGICAL CONJA,CONJB,NOTA,NOTB
 *     ..  *     ..
 *     .. Parameters ..  *     .. Parameters ..
       DOUBLE COMPLEX ONE        COMPLEX*16 ONE
       PARAMETER (ONE= (1.0D+0,0.0D+0))        PARAMETER (ONE= (1.0D+0,0.0D+0))
       DOUBLE COMPLEX ZERO        COMPLEX*16 ZERO
       PARAMETER (ZERO= (0.0D+0,0.0D+0))        PARAMETER (ZERO= (0.0D+0,0.0D+0))
 *     ..  *     ..
 *  *
 *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not  *     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
 *     conjugated or transposed, set  CONJA and CONJB  as true if  A  and  *     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
 *     B  respectively are to be  transposed but  not conjugated  and set  *     B  respectively are to be  transposed but  not conjugated  and set
 *     NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A  *     NROWA and NROWB  as the number of rows  of  A  and  B  respectively.
 *     and the number of rows of  B  respectively.  
 *  *
       NOTA = LSAME(TRANSA,'N')        NOTA = LSAME(TRANSA,'N')
       NOTB = LSAME(TRANSB,'N')        NOTB = LSAME(TRANSB,'N')
Line 165 Line 233
       CONJB = LSAME(TRANSB,'C')        CONJB = LSAME(TRANSB,'C')
       IF (NOTA) THEN        IF (NOTA) THEN
           NROWA = M            NROWA = M
           NCOLA = K  
       ELSE        ELSE
           NROWA = K            NROWA = K
           NCOLA = M  
       END IF        END IF
       IF (NOTB) THEN        IF (NOTB) THEN
           NROWB = K            NROWB = K
Line 245 Line 311
    60                 CONTINUE     60                 CONTINUE
                   END IF                    END IF
                   DO 80 L = 1,K                    DO 80 L = 1,K
                       IF (B(L,J).NE.ZERO) THEN                        TEMP = ALPHA*B(L,J)
                           TEMP = ALPHA*B(L,J)                        DO 70 I = 1,M
                           DO 70 I = 1,M                            C(I,J) = C(I,J) + TEMP*A(I,L)
                               C(I,J) = C(I,J) + TEMP*A(I,L)     70                 CONTINUE
    70                     CONTINUE  
                       END IF  
    80             CONTINUE     80             CONTINUE
    90         CONTINUE     90         CONTINUE
           ELSE IF (CONJA) THEN            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 120 J = 1,N
                   DO 110 I = 1,M                    DO 110 I = 1,M
Line 272 Line 336
   120         CONTINUE    120         CONTINUE
           ELSE            ELSE
 *  *
 *           Form  C := alpha*A'*B + beta*C  *           Form  C := alpha*A**T*B + beta*C
 *  *
               DO 150 J = 1,N                DO 150 J = 1,N
                   DO 140 I = 1,M                    DO 140 I = 1,M
Line 291 Line 355
       ELSE IF (NOTA) THEN        ELSE IF (NOTA) THEN
           IF (CONJB) 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                DO 200 J = 1,N
                   IF (BETA.EQ.ZERO) THEN                    IF (BETA.EQ.ZERO) THEN
Line 304 Line 368
   170                 CONTINUE    170                 CONTINUE
                   END IF                    END IF
                   DO 190 L = 1,K                    DO 190 L = 1,K
                       IF (B(J,L).NE.ZERO) THEN                        TEMP = ALPHA*DCONJG(B(J,L))
                           TEMP = ALPHA*DCONJG(B(J,L))                        DO 180 I = 1,M
                           DO 180 I = 1,M                            C(I,J) = C(I,J) + TEMP*A(I,L)
                               C(I,J) = C(I,J) + TEMP*A(I,L)    180                 CONTINUE
   180                     CONTINUE  
                       END IF  
   190             CONTINUE    190             CONTINUE
   200         CONTINUE    200         CONTINUE
           ELSE            ELSE
 *  *
 *           Form  C := alpha*A*B'          + beta*C  *           Form  C := alpha*A*B**T + beta*C
 *  *
               DO 250 J = 1,N                DO 250 J = 1,N
                   IF (BETA.EQ.ZERO) THEN                    IF (BETA.EQ.ZERO) THEN
Line 327 Line 389
   220                 CONTINUE    220                 CONTINUE
                   END IF                    END IF
                   DO 240 L = 1,K                    DO 240 L = 1,K
                       IF (B(J,L).NE.ZERO) THEN                        TEMP = ALPHA*B(J,L)
                           TEMP = ALPHA*B(J,L)                        DO 230 I = 1,M
                           DO 230 I = 1,M                            C(I,J) = C(I,J) + TEMP*A(I,L)
                               C(I,J) = C(I,J) + TEMP*A(I,L)    230                 CONTINUE
   230                     CONTINUE  
                       END IF  
   240             CONTINUE    240             CONTINUE
   250         CONTINUE    250         CONTINUE
           END IF            END IF
       ELSE IF (CONJA) THEN        ELSE IF (CONJA) THEN
           IF (CONJB) 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 280 J = 1,N
                   DO 270 I = 1,M                    DO 270 I = 1,M
Line 356 Line 416
   280         CONTINUE    280         CONTINUE
           ELSE            ELSE
 *  *
 *           Form  C := alpha*conjg( A' )*B' + beta*C  *           Form  C := alpha*A**H*B**T + beta*C
 *  *
               DO 310 J = 1,N                DO 310 J = 1,N
                   DO 300 I = 1,M                    DO 300 I = 1,M
Line 375 Line 435
       ELSE        ELSE
           IF (CONJB) THEN            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 340 J = 1,N
                   DO 330 I = 1,M                    DO 330 I = 1,M
Line 392 Line 452
   340         CONTINUE    340         CONTINUE
           ELSE            ELSE
 *  *
 *           Form  C := alpha*A'*B' + beta*C  *           Form  C := alpha*A**T*B**T + beta*C
 *  *
               DO 370 J = 1,N                DO 370 J = 1,N
                   DO 360 I = 1,M                    DO 360 I = 1,M
Line 412 Line 472
 *  *
       RETURN        RETURN
 *  *
 *     End of ZGEMM .  *     End of ZGEMM
 *  *
       END        END

Removed from v.1.6  
changed lines
  Added in v.1.17


CVSweb interface <joel.bertrand@systella.fr>