share/lapack/blas/fortran/zhemm.f

   1       SUBROUTINE ZHEMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB,
   2      $                   BETA, C, LDC )
   3 *     .. Scalar Arguments ..
   4       CHARACTER*1        SIDE, UPLO
   5       INTEGER            M, N, LDA, LDB, LDC
   6       COMPLEX*16         ALPHA, BETA
   7 *     .. Array Arguments ..
   8       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  ZHEMM  performs one of the matrix-matrix operations
  15 *
  16 *     C := alpha*A*B + beta*C,
  17 *
  18 *  or
  19 *
  20 *     C := alpha*B*A + beta*C,
  21 *
  22 *  where alpha and beta are scalars, A is an hermitian matrix and  B and
  23 *  C are m by n matrices.
  24 *
  25 *  Parameters
  26 *  ==========
  27 *
  28 *  SIDE   - CHARACTER*1.
  29 *           On entry,  SIDE  specifies whether  the  hermitian matrix  A
  30 *           appears on the  left or right  in the  operation as follows:
  31 *
  32 *              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
  33 *
  34 *              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
  35 *
  36 *           Unchanged on exit.
  37 *
  38 *  UPLO   - CHARACTER*1.
  39 *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
  40 *           triangular  part  of  the  hermitian  matrix   A  is  to  be
  41 *           referenced as follows:
  42 *
  43 *              UPLO = 'U' or 'u'   Only the upper triangular part of the
  44 *                                  hermitian matrix is to be referenced.
  45 *
  46 *              UPLO = 'L' or 'l'   Only the lower triangular part of the
  47 *                                  hermitian matrix is to be referenced.
  48 *
  49 *           Unchanged on exit.
  50 *
  51 *  M      - INTEGER.
  52 *           On entry,  M  specifies the number of rows of the matrix  C.
  53 *           M  must be at least zero.
  54 *           Unchanged on exit.
  55 *
  56 *  N      - INTEGER.
  57 *           On entry, N specifies the number of columns of the matrix C.
  58 *           N  must be at least zero.
  59 *           Unchanged on exit.
  60 *
  61 *  ALPHA  - COMPLEX*16      .
  62 *           On entry, ALPHA specifies the scalar alpha.
  63 *           Unchanged on exit.
  64 *
  65 *  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
  66 *           m  when  SIDE = 'L' or 'l'  and is n  otherwise.
  67 *           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
  68 *           the array  A  must contain the  hermitian matrix,  such that
  69 *           when  UPLO = 'U' or 'u', the leading m by m upper triangular
  70 *           part of the array  A  must contain the upper triangular part
  71 *           of the  hermitian matrix and the  strictly  lower triangular
  72 *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  73 *           the leading  m by m  lower triangular part  of the  array  A
  74 *           must  contain  the  lower triangular part  of the  hermitian
  75 *           matrix and the  strictly upper triangular part of  A  is not
  76 *           referenced.
  77 *           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
  78 *           the array  A  must contain the  hermitian matrix,  such that
  79 *           when  UPLO = 'U' or 'u', the leading n by n upper triangular
  80 *           part of the array  A  must contain the upper triangular part
  81 *           of the  hermitian matrix and the  strictly  lower triangular
  82 *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
  83 *           the leading  n by n  lower triangular part  of the  array  A
  84 *           must  contain  the  lower triangular part  of the  hermitian
  85 *           matrix and the  strictly upper triangular part of  A  is not
  86 *           referenced.
  87 *           Note that the imaginary parts  of the diagonal elements need
  88 *           not be set, they are assumed to be zero.
  89 *           Unchanged on exit.
  90 *
  91 *  LDA    - INTEGER.
  92 *           On entry, LDA specifies the first dimension of A as declared
  93 *           in the  calling (sub) program. When  SIDE = 'L' or 'l'  then
  94 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
  95 *           least max( 1, n ).
  96 *           Unchanged on exit.
  97 *
  98 *  B      - COMPLEX*16       array of DIMENSION ( LDB, n ).
  99 *           Before entry, the leading  m by n part of the array  B  must
 100 *           contain the matrix B.
 101 *           Unchanged on exit.
 102 *
 103 *  LDB    - INTEGER.
 104 *           On entry, LDB specifies the first dimension of B as declared
 105 *           in  the  calling  (sub)  program.   LDB  must  be  at  least
 106 *           max( 1, m ).
 107 *           Unchanged on exit.
 108 *
 109 *  BETA   - COMPLEX*16      .
 110 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
 111 *           supplied as zero then C need not be set on input.
 112 *           Unchanged on exit.
 113 *
 114 *  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
 115 *           Before entry, the leading  m by n  part of the array  C must
 116 *           contain the matrix  C,  except when  beta  is zero, in which
 117 *           case C need not be set on entry.
 118 *           On exit, the array  C  is overwritten by the  m by n updated
 119 *           matrix.
 120 *
 121 *  LDC    - INTEGER.
 122 *           On entry, LDC specifies the first dimension of C as declared
 123 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
 124 *           max( 1, m ).
 125 *           Unchanged on exit.
 126 *
 127 *
 128 *  Level 3 Blas routine.
 129 *
 130 *  -- Written on 8-February-1989.
 131 *     Jack Dongarra, Argonne National Laboratory.
 132 *     Iain Duff, AERE Harwell.
 133 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 134 *     Sven Hammarling, Numerical Algorithms Group Ltd.
 135 *
 136 *
 137 *     .. External Functions ..
 138       LOGICAL            LSAME
 139       EXTERNAL           LSAME
 140 *     .. External Subroutines ..
 141       EXTERNAL           XERBLA
 142 *     .. Intrinsic Functions ..
 143       INTRINSIC          DCONJG, MAX, DBLE
 144 *     .. Local Scalars ..
 145       LOGICAL            UPPER
 146       INTEGER            I, INFO, J, K, NROWA
 147       COMPLEX*16         TEMP1, TEMP2
 148 *     .. Parameters ..
 149       COMPLEX*16         ONE
 150       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
 151       COMPLEX*16         ZERO
 152       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
 153 *     ..
 154 *     .. Executable Statements ..
 155 *
 156 *     Set NROWA as the number of rows of A.
 157 *
 158       IF( LSAME( SIDE, 'L' ) )THEN
 159          NROWA = M
 160       ELSE
 161          NROWA = N
 162       END IF
 163       UPPER = LSAME( UPLO, 'U' )
 164 *
 165 *     Test the input parameters.
 166 *
 167       INFO = 0
 168       IF(      ( .NOT.LSAME( SIDE, 'L' ) ).AND.
 169      $         ( .NOT.LSAME( SIDE, 'R' ) )      )THEN
 170          INFO = 1
 171       ELSE IF( ( .NOT.UPPER              ).AND.
 172      $         ( .NOT.LSAME( UPLO, 'L' ) )      )THEN
 173          INFO = 2
 174       ELSE IF( M  .LT.0               )THEN
 175          INFO = 3
 176       ELSE IF( N  .LT.0               )THEN
 177          INFO = 4
 178       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
 179          INFO = 7
 180       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
 181          INFO = 9
 182       ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
 183          INFO = 12
 184       END IF
 185       IF( INFO.NE.0 )THEN
 186          CALL XERBLA( 'ZHEMM ', INFO )
 187          RETURN
 188       END IF
 189 *
 190 *     Quick return if possible.
 191 *
 192       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
 193      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
 194      $   RETURN
 195 *
 196 *     And when  alpha.eq.zero.
 197 *
 198       IF( ALPHA.EQ.ZERO )THEN
 199          IF( BETA.EQ.ZERO )THEN
 200             DO 20, J = 1, N
 201                DO 10, I = 1, M
 202                   C( I, J ) = ZERO
 203    10          CONTINUE
 204    20       CONTINUE
 205          ELSE
 206             DO 40, J = 1, N
 207                DO 30, I = 1, M
 208                   C( I, J ) = BETA*C( I, J )
 209    30          CONTINUE
 210    40       CONTINUE
 211          END IF
 212          RETURN
 213       END IF
 214 *
 215 *     Start the operations.
 216 *
 217       IF( LSAME( SIDE, 'L' ) )THEN
 218 *
 219 *        Form  C := alpha*A*B + beta*C.
 220 *
 221          IF( UPPER )THEN
 222             DO 70, J = 1, N
 223                DO 60, I = 1, M
 224                   TEMP1 = ALPHA*B( I, J )
 225                   TEMP2 = ZERO
 226                   DO 50, K = 1, I - 1
 227                      C( K, J ) = C( K, J ) + TEMP1*A( K, I )
 228                      TEMP2     = TEMP2     +
 229      $                           B( K, J )*DCONJG( A( K, I ) )
 230    50             CONTINUE
 231                   IF( BETA.EQ.ZERO )THEN
 232                      C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
 233      $                           ALPHA*TEMP2
 234                   ELSE
 235                      C( I, J ) = BETA *C( I, J )         +
 236      $                           TEMP1*DBLE( A( I, I ) ) +
 237      $                           ALPHA*TEMP2
 238                   END IF
 239    60          CONTINUE
 240    70       CONTINUE
 241          ELSE
 242             DO 100, J = 1, N
 243                DO 90, I = M, 1, -1
 244                   TEMP1 = ALPHA*B( I, J )
 245                   TEMP2 = ZERO
 246                   DO 80, K = I + 1, M
 247                      C( K, J ) = C( K, J ) + TEMP1*A( K, I )
 248                      TEMP2     = TEMP2     +
 249      $                           B( K, J )*DCONJG( A( K, I ) )
 250    80             CONTINUE
 251                   IF( BETA.EQ.ZERO )THEN
 252                      C( I, J ) = TEMP1*DBLE( A( I, I ) ) +
 253      $                           ALPHA*TEMP2
 254                   ELSE
 255                      C( I, J ) = BETA *C( I, J )         +
 256      $                           TEMP1*DBLE( A( I, I ) ) +
 257      $                           ALPHA*TEMP2
 258                   END IF
 259    90          CONTINUE
 260   100       CONTINUE
 261          END IF
 262       ELSE
 263 *
 264 *        Form  C := alpha*B*A + beta*C.
 265 *
 266          DO 170, J = 1, N
 267             TEMP1 = ALPHA*DBLE( A( J, J ) )
 268             IF( BETA.EQ.ZERO )THEN
 269                DO 110, I = 1, M
 270                   C( I, J ) = TEMP1*B( I, J )
 271   110          CONTINUE
 272             ELSE
 273                DO 120, I = 1, M
 274                   C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J )
 275   120          CONTINUE
 276             END IF
 277             DO 140, K = 1, J - 1
 278                IF( UPPER )THEN
 279                   TEMP1 = ALPHA*A( K, J )
 280                ELSE
 281                   TEMP1 = ALPHA*DCONJG( A( J, K ) )
 282                END IF
 283                DO 130, I = 1, M
 284                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
 285   130          CONTINUE
 286   140       CONTINUE
 287             DO 160, K = J + 1, N
 288                IF( UPPER )THEN
 289                   TEMP1 = ALPHA*DCONJG( A( J, K ) )
 290                ELSE
 291                   TEMP1 = ALPHA*A( K, J )
 292                END IF
 293                DO 150, I = 1, M
 294                   C( I, J ) = C( I, J ) + TEMP1*B( I, K )
 295   150          CONTINUE
 296   160       CONTINUE
 297   170    CONTINUE
 298       END IF
 299 *
 300       RETURN
 301 *
 302 *     End of ZHEMM .
 303 *
 304       END