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 / sspr.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  5.9 KB  |  202 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE SSPR  ( UPLO, N, ALPHA, X, INCX, AP )
  5. *     .. Scalar Arguments ..
  6.       REAL               ALPHA
  7.       INTEGER            INCX, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       REAL               AP( * ), X( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  SSPR    performs the symmetric rank 1 operation
  17. *
  18. *     A := alpha*x*x' + A,
  19. *
  20. *  where alpha is a real scalar, x is an n element vector and A is an
  21. *  n by n symmetric 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  - REAL            .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - REAL             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     - REAL             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 symmetric 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 symmetric 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. *
  76. *
  77. *  Level 2 Blas routine.
  78. *
  79. *  -- Written on 22-October-1986.
  80. *     Jack Dongarra, Argonne National Lab.
  81. *     Jeremy Du Croz, Nag Central Office.
  82. *     Sven Hammarling, Nag Central Office.
  83. *     Richard Hanson, Sandia National Labs.
  84. *
  85. *
  86. *     .. Parameters ..
  87.       REAL               ZERO
  88.       PARAMETER        ( ZERO = 0.0E+0 )
  89. *     .. Local Scalars ..
  90.       REAL               TEMP
  91.       INTEGER            I, INFO, IX, J, JX, K, KK, KX
  92. *     .. External Functions ..
  93.       LOGICAL            LSAME
  94.       EXTERNAL           LSAME
  95. *     .. External Subroutines ..
  96.       EXTERNAL           XERBLA
  97. *     ..
  98. *     .. Executable Statements ..
  99. *
  100. *     Test the input parameters.
  101. *
  102.       INFO = 0
  103.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  104.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  105.          INFO = 1
  106.       ELSE IF( N.LT.0 )THEN
  107.          INFO = 2
  108.       ELSE IF( INCX.EQ.0 )THEN
  109.          INFO = 5
  110.       END IF
  111.       IF( INFO.NE.0 )THEN
  112.          CALL XERBLA( 'SSPR  ', INFO )
  113.          RETURN
  114.       END IF
  115. *
  116. *     Quick return if possible.
  117. *
  118.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  119.      $   RETURN
  120. *
  121. *     Set the start point in X if the increment is not unity.
  122. *
  123.       IF( INCX.LE.0 )THEN
  124.          KX = 1 - ( N - 1 )*INCX
  125.       ELSE IF( INCX.NE.1 )THEN
  126.          KX = 1
  127.       END IF
  128. *
  129. *     Start the operations. In this version the elements of the array AP
  130. *     are accessed sequentially with one pass through AP.
  131. *
  132.       KK = 1
  133.       IF( LSAME( UPLO, 'U' ) )THEN
  134. *
  135. *        Form  A  when upper triangle is stored in AP.
  136. *
  137.          IF( INCX.EQ.1 )THEN
  138.             DO 20, J = 1, N
  139.                IF( X( J ).NE.ZERO )THEN
  140.                   TEMP = ALPHA*X( J )
  141.                   K    = KK
  142.                   DO 10, I = 1, J
  143.                      AP( K ) = AP( K ) + X( I )*TEMP
  144.                      K       = K       + 1
  145.    10             CONTINUE
  146.                END IF
  147.                KK = KK + J
  148.    20       CONTINUE
  149.          ELSE
  150.             JX = KX
  151.             DO 40, J = 1, N
  152.                IF( X( JX ).NE.ZERO )THEN
  153.                   TEMP = ALPHA*X( JX )
  154.                   IX   = KX
  155.                   DO 30, K = KK, KK + J - 1
  156.                      AP( K ) = AP( K ) + X( IX )*TEMP
  157.                      IX      = IX      + INCX
  158.    30             CONTINUE
  159.                END IF
  160.                JX = JX + INCX
  161.                KK = KK + J
  162.    40       CONTINUE
  163.          END IF
  164.       ELSE
  165. *
  166. *        Form  A  when lower triangle is stored in AP.
  167. *
  168.          IF( INCX.EQ.1 )THEN
  169.             DO 60, J = 1, N
  170.                IF( X( J ).NE.ZERO )THEN
  171.                   TEMP = ALPHA*X( J )
  172.                   K    = KK
  173.                   DO 50, I = J, N
  174.                      AP( K ) = AP( K ) + X( I )*TEMP
  175.                      K       = K       + 1
  176.    50             CONTINUE
  177.                END IF
  178.                KK = KK + N - J + 1
  179.    60       CONTINUE
  180.          ELSE
  181.             JX = KX
  182.             DO 80, J = 1, N
  183.                IF( X( JX ).NE.ZERO )THEN
  184.                   TEMP = ALPHA*X( JX )
  185.                   IX   = JX
  186.                   DO 70, K = KK, KK + N - J
  187.                      AP( K ) = AP( K ) + X( IX )*TEMP
  188.                      IX      = IX      + INCX
  189.    70             CONTINUE
  190.                END IF
  191.                JX = JX + INCX
  192.                KK = KK + N - J + 1
  193.    80       CONTINUE
  194.          END IF
  195.       END IF
  196. *
  197.       RETURN
  198. *
  199. *     End of SSPR  .
  200. *
  201.       END
  202.