zhemm.f

Go to the documentation of this file.

*> \brief \b ZHEMM
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZHEMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
*       .. Scalar Arguments ..
*       COMPLEX*16 ALPHA,BETA
*       INTEGER LDA,LDB,LDC,M,N
*       CHARACTER SIDE,UPLO
*       ..
*       .. Array Arguments ..
*       COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZHEMM  performs one of the matrix-matrix operations
*>
*>    C := alpha*A*B + beta*C,
*>
*> or
*>
*>    C := alpha*B*A + beta*C,
*>
*> where alpha and beta are scalars, A is an hermitian matrix and  B and
*> C are m by n matrices.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] SIDE
*> \verbatim
*>          SIDE is CHARACTER*1
*>           On entry,  SIDE  specifies whether  the  hermitian matrix  A
*>           appears on the  left or right  in the  operation as follows:
*>
*>              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
*>
*>              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
*> \endverbatim
*>
*> \param[in] UPLO
*> \verbatim
*>          UPLO is CHARACTER*1
*>           On  entry,   UPLO  specifies  whether  the  upper  or  lower
*>           triangular  part  of  the  hermitian  matrix   A  is  to  be
*>           referenced as follows:
*>
*>              UPLO = 'U' or 'u'   Only the upper triangular part of the
*>                                  hermitian matrix is to be referenced.
*>
*>              UPLO = 'L' or 'l'   Only the lower triangular part of the
*>                                  hermitian matrix is to be referenced.
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>           On entry,  M  specifies the number of rows 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 C.
*>           N  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
*>           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
*>           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
*>           the array  A  must contain the  hermitian matrix,  such that
*>           when  UPLO = 'U' or 'u', the leading m by m upper triangular
*>           part of the array  A  must contain the upper triangular part
*>           of the  hermitian matrix and the  strictly  lower triangular
*>           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
*>           the leading  m by m  lower triangular part  of the  array  A
*>           must  contain  the  lower triangular part  of the  hermitian
*>           matrix and the  strictly upper triangular part of  A  is not
*>           referenced.
*>           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
*>           the array  A  must contain the  hermitian matrix,  such that
*>           when  UPLO = 'U' or 'u', the leading n by n upper triangular
*>           part of the array  A  must contain the upper triangular part
*>           of the  hermitian matrix and the  strictly  lower triangular
*>           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
*>           the leading  n by n  lower triangular part  of the  array  A
*>           must  contain  the  lower triangular part  of the  hermitian
*>           matrix and the  strictly upper triangular part of  A  is not
*>           referenced.
*>           Note that the imaginary parts  of the diagonal elements need
*>           not be set, they are assumed to be zero.
*> \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  SIDE = 'L' or 'l'  then
*>           LDA must be at least  max( 1, m ), otherwise  LDA must be at
*>           least max( 1, n ).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is COMPLEX*16 array, dimension ( LDB, N )
*>           Before entry, the leading  m by n 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.   LDB  must  be  at  least
*>           max( 1, m ).
*> \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 updated
*>           matrix.
*> \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 hemm
*
*> \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 zhemm(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
*
*  -- Reference BLAS level3 routine --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      COMPLEX*16 ALPHA,BETA
      INTEGER LDA,LDB,LDC,M,N
      CHARACTER SIDE,UPLO
*     ..
*     .. Array Arguments ..
      COMPLEX*16 A(LDA,*),B(LDB,*),C(LDC,*)
*     ..
*
*  =====================================================================
*
*     .. External Functions ..
      LOGICAL LSAME
      EXTERNAL lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC dble,dconjg,max
*     ..
*     .. Local Scalars ..
      COMPLEX*16 TEMP1,TEMP2
      INTEGER I,INFO,J,K,NROWA
      LOGICAL UPPER
*     ..
*     .. Parameters ..
      COMPLEX*16 ONE
      parameter(one= (1.0d+0,0.0d+0))
      COMPLEX*16 ZERO
      parameter(zero= (0.0d+0,0.0d+0))
*     ..
*
*     Set NROWA as the number of rows of A.
*
      IF (lsame(side,'L')) THEN
          nrowa = m
      ELSE
          nrowa = n
      END IF
      upper = lsame(uplo,'U')
*
*     Test the input parameters.
*
      info = 0
      IF ((.NOT.lsame(side,'L')) .AND.
     +    (.NOT.lsame(side,'R'))) THEN
          info = 1
      ELSE IF ((.NOT.upper) .AND.
     +         (.NOT.lsame(uplo,'L'))) THEN
          info = 2
      ELSE IF (m.LT.0) THEN
          info = 3
      ELSE IF (n.LT.0) THEN
          info = 4
      ELSE IF (lda.LT.max(1,nrowa)) THEN
          info = 7
      ELSE IF (ldb.LT.max(1,m)) THEN
          info = 9
      ELSE IF (ldc.LT.max(1,m)) THEN
          info = 12
      END IF
      IF (info.NE.0) THEN
          CALL xerbla('ZHEMM ',info)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((m.EQ.0) .OR. (n.EQ.0) .OR.
     +    ((alpha.EQ.zero).AND. (beta.EQ.one))) RETURN
*
*     And when  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 (lsame(side,'L')) THEN
*
*        Form  C := alpha*A*B + beta*C.
*
          IF (upper) THEN
              DO 70 j = 1,n
                  DO 60 i = 1,m
                      temp1 = alpha*b(i,j)
                      temp2 = zero
                      DO 50 k = 1,i - 1
                          c(k,j) = c(k,j) + temp1*a(k,i)
                          temp2 = temp2 + b(k,j)*dconjg(a(k,i))
   50                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = temp1*dble(a(i,i)) + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + temp1*dble(a(i,i)) +
     +                             alpha*temp2
                      END IF
   60             CONTINUE
   70         CONTINUE
          ELSE
              DO 100 j = 1,n
                  DO 90 i = m,1,-1
                      temp1 = alpha*b(i,j)
                      temp2 = zero
                      DO 80 k = i + 1,m
                          c(k,j) = c(k,j) + temp1*a(k,i)
                          temp2 = temp2 + b(k,j)*dconjg(a(k,i))
   80                 CONTINUE
                      IF (beta.EQ.zero) THEN
                          c(i,j) = temp1*dble(a(i,i)) + alpha*temp2
                      ELSE
                          c(i,j) = beta*c(i,j) + temp1*dble(a(i,i)) +
     +                             alpha*temp2
                      END IF
   90             CONTINUE
  100         CONTINUE
          END IF
      ELSE
*
*        Form  C := alpha*B*A + beta*C.
*
          DO 170 j = 1,n
              temp1 = alpha*dble(a(j,j))
              IF (beta.EQ.zero) THEN
                  DO 110 i = 1,m
                      c(i,j) = temp1*b(i,j)
  110             CONTINUE
              ELSE
                  DO 120 i = 1,m
                      c(i,j) = beta*c(i,j) + temp1*b(i,j)
  120             CONTINUE
              END IF
              DO 140 k = 1,j - 1
                  IF (upper) THEN
                      temp1 = alpha*a(k,j)
                  ELSE
                      temp1 = alpha*dconjg(a(j,k))
                  END IF
                  DO 130 i = 1,m
                      c(i,j) = c(i,j) + temp1*b(i,k)
  130             CONTINUE
  140         CONTINUE
              DO 160 k = j + 1,n
                  IF (upper) THEN
                      temp1 = alpha*dconjg(a(j,k))
                  ELSE
                      temp1 = alpha*a(k,j)
                  END IF
                  DO 150 i = 1,m
                      c(i,j) = c(i,j) + temp1*b(i,k)
  150             CONTINUE
  160         CONTINUE
  170     CONTINUE
      END IF
*
      RETURN
*
*     End of ZHEMM
*
      END