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 / zhpr.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  6.9 KB  |  221 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZHPR  ( UPLO, N, ALPHA, X, INCX, AP )
  5. *     .. Scalar Arguments ..
  6.       DOUBLE PRECISION   ALPHA
  7.       INTEGER            INCX, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       COMPLEX*16         AP( * ), X( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  ZHPR    performs the hermitian rank 1 operation
  17. *
  18. *     A := alpha*x*conjg( x' ) + A,
  19. *
  20. *  where alpha is a real scalar, x is an n element vector and A is an
  21. *  n by n hermitian matrix, supplied in packed form.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the matrix A is supplied in the packed
  29. *           array AP as follows:
  30. *
  31. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  32. *                                  supplied in AP.
  33. *
  34. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  35. *                                  supplied in AP.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - DOUBLE PRECISION.
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - COMPLEX*16       array of dimension at least
  49. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  50. *           Before entry, the incremented array X must contain the n
  51. *           element vector x.
  52. *           Unchanged on exit.
  53. *
  54. *  INCX   - INTEGER.
  55. *           On entry, INCX specifies the increment for the elements of
  56. *           X. INCX must not be zero.
  57. *           Unchanged on exit.
  58. *
  59. *  AP     - COMPLEX*16       array of DIMENSION at least
  60. *           ( ( n*( n + 1 ) )/2 ).
  61. *           Before entry with  UPLO = 'U' or 'u', the array AP must
  62. *           contain the upper triangular part of the hermitian matrix
  63. *           packed sequentially, column by column, so that AP( 1 )
  64. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  65. *           and a( 2, 2 ) respectively, and so on. On exit, the array
  66. *           AP is overwritten by the upper triangular part of the
  67. *           updated matrix.
  68. *           Before entry with UPLO = 'L' or 'l', the array AP must
  69. *           contain the lower triangular part of the hermitian matrix
  70. *           packed sequentially, column by column, so that AP( 1 )
  71. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  72. *           and a( 3, 1 ) respectively, and so on. On exit, the array
  73. *           AP is overwritten by the lower triangular part of the
  74. *           updated matrix.
  75. *           Note that the imaginary parts of the diagonal elements need
  76. *           not be set, they are assumed to be zero, and on exit they
  77. *           are set to zero.
  78. *
  79. *
  80. *  Level 2 Blas routine.
  81. *
  82. *  -- Written on 22-October-1986.
  83. *     Jack Dongarra, Argonne National Lab.
  84. *     Jeremy Du Croz, Nag Central Office.
  85. *     Sven Hammarling, Nag Central Office.
  86. *     Richard Hanson, Sandia National Labs.
  87. *
  88. *
  89. *     .. Parameters ..
  90.       COMPLEX*16         ZERO
  91.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  92. *     .. Local Scalars ..
  93.       COMPLEX*16         TEMP
  94.       INTEGER            I, INFO, IX, J, JX, K, KK, KX
  95. *     .. External Functions ..
  96.       LOGICAL            LSAME
  97.       EXTERNAL           LSAME
  98. *     .. External Subroutines ..
  99.       EXTERNAL           XERBLA
  100. *     .. Intrinsic Functions ..
  101.       INTRINSIC          DCONJG, DBLE
  102. *     ..
  103. *     .. Executable Statements ..
  104. *
  105. *     Test the input parameters.
  106. *
  107.       INFO = 0
  108.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  109.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  110.          INFO = 1
  111.       ELSE IF( N.LT.0 )THEN
  112.          INFO = 2
  113.       ELSE IF( INCX.EQ.0 )THEN
  114.          INFO = 5
  115.       END IF
  116.       IF( INFO.NE.0 )THEN
  117.          CALL XERBLA( 'ZHPR  ', INFO )
  118.          RETURN
  119.       END IF
  120. *
  121. *     Quick return if possible.
  122. *
  123.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) )
  124.      $   RETURN
  125. *
  126. *     Set the start point in X if the increment is not unity.
  127. *
  128.       IF( INCX.LE.0 )THEN
  129.          KX = 1 - ( N - 1 )*INCX
  130.       ELSE IF( INCX.NE.1 )THEN
  131.          KX = 1
  132.       END IF
  133. *
  134. *     Start the operations. In this version the elements of the array AP
  135. *     are accessed sequentially with one pass through AP.
  136. *
  137.       KK = 1
  138.       IF( LSAME( UPLO, 'U' ) )THEN
  139. *
  140. *        Form  A  when upper triangle is stored in AP.
  141. *
  142.          IF( INCX.EQ.1 )THEN
  143.             DO 20, J = 1, N
  144.                IF( X( J ).NE.ZERO )THEN
  145.                   TEMP = ALPHA*DCONJG( X( J ) )
  146.                   K    = KK
  147.                   DO 10, I = 1, J - 1
  148.                      AP( K ) = AP( K ) + X( I )*TEMP
  149.                      K       = K       + 1
  150.    10             CONTINUE
  151.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  152.      $                               + DBLE( X( J )*TEMP )
  153.                ELSE
  154.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  155.                END IF
  156.                KK = KK + J
  157.    20       CONTINUE
  158.          ELSE
  159.             JX = KX
  160.             DO 40, J = 1, N
  161.                IF( X( JX ).NE.ZERO )THEN
  162.                   TEMP = ALPHA*DCONJG( X( JX ) )
  163.                   IX   = KX
  164.                   DO 30, K = KK, KK + J - 2
  165.                      AP( K ) = AP( K ) + X( IX )*TEMP
  166.                      IX      = IX      + INCX
  167.    30             CONTINUE
  168.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  169.      $                               + DBLE( X( JX )*TEMP )
  170.                ELSE
  171.                   AP( KK + J - 1 ) = DBLE( AP( KK + J - 1 ) )
  172.                END IF
  173.                JX = JX + INCX
  174.                KK = KK + J
  175.    40       CONTINUE
  176.          END IF
  177.       ELSE
  178. *
  179. *        Form  A  when lower triangle is stored in AP.
  180. *
  181.          IF( INCX.EQ.1 )THEN
  182.             DO 60, J = 1, N
  183.                IF( X( J ).NE.ZERO )THEN
  184.                   TEMP     = ALPHA*DCONJG( X( J ) )
  185.                   AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( J ) )
  186.                   K        = KK               + 1
  187.                   DO 50, I = J + 1, N
  188.                      AP( K ) = AP( K ) + X( I )*TEMP
  189.                      K       = K       + 1
  190.    50             CONTINUE
  191.                ELSE
  192.                   AP( KK ) = DBLE( AP( KK ) )
  193.                END IF
  194.                KK = KK + N - J + 1
  195.    60       CONTINUE
  196.          ELSE
  197.             JX = KX
  198.             DO 80, J = 1, N
  199.                IF( X( JX ).NE.ZERO )THEN
  200.                   TEMP    = ALPHA*DCONJG( X( JX ) )
  201.                   AP( KK ) = DBLE( AP( KK ) ) + DBLE( TEMP*X( JX ) )
  202.                   IX      = JX
  203.                   DO 70, K = KK + 1, KK + N - J
  204.                      IX      = IX      + INCX
  205.                      AP( K ) = AP( K ) + X( IX )*TEMP
  206.    70             CONTINUE
  207.                ELSE
  208.                   AP( KK ) = DBLE( AP( KK ) )
  209.                END IF
  210.                JX = JX + INCX
  211.                KK = KK + N - J + 1
  212.    80       CONTINUE
  213.          END IF
  214.       END IF
  215. *
  216.       RETURN
  217. *
  218. *     End of ZHPR  .
  219. *
  220.       END
  221.