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