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