LAPACK: BLAS/SRC/sspr2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
00002 *     .. Scalar Arguments ..
00003       REAL ALPHA
00004       INTEGER INCX,INCY,N
00005       CHARACTER UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       REAL AP(*),X(*),Y(*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  SSPR2  performs the symmetric rank 2 operation
00015 *
00016 *     A := alpha*x*y' + alpha*y*x' + A,
00017 *
00018 *  where alpha is a scalar, x and y are n element vectors and A is an
00019 *  n by n symmetric matrix, supplied in packed form.
00020 *
00021 *  Arguments
00022 *  ==========
00023 *
00024 *  UPLO   - CHARACTER*1.
00025 *           On entry, UPLO specifies whether the upper or lower
00026 *           triangular part of the matrix A is supplied in the packed
00027 *           array AP as follows:
00028 *
00029 *              UPLO = 'U' or 'u'   The upper triangular part of A is
00030 *                                  supplied in AP.
00031 *
00032 *              UPLO = 'L' or 'l'   The lower triangular part of A is
00033 *                                  supplied in AP.
00034 *
00035 *           Unchanged on exit.
00036 *
00037 *  N      - INTEGER.
00038 *           On entry, N specifies the order of the matrix A.
00039 *           N must be at least zero.
00040 *           Unchanged on exit.
00041 *
00042 *  ALPHA  - REAL            .
00043 *           On entry, ALPHA specifies the scalar alpha.
00044 *           Unchanged on exit.
00045 *
00046 *  X      - REAL             array of dimension at least
00047 *           ( 1 + ( n - 1 )*abs( INCX ) ).
00048 *           Before entry, the incremented array X must contain the n
00049 *           element vector x.
00050 *           Unchanged on exit.
00051 *
00052 *  INCX   - INTEGER.
00053 *           On entry, INCX specifies the increment for the elements of
00054 *           X. INCX must not be zero.
00055 *           Unchanged on exit.
00056 *
00057 *  Y      - REAL             array of dimension at least
00058 *           ( 1 + ( n - 1 )*abs( INCY ) ).
00059 *           Before entry, the incremented array Y must contain the n
00060 *           element vector y.
00061 *           Unchanged on exit.
00062 *
00063 *  INCY   - INTEGER.
00064 *           On entry, INCY specifies the increment for the elements of
00065 *           Y. INCY must not be zero.
00066 *           Unchanged on exit.
00067 *
00068 *  AP     - REAL             array of DIMENSION at least
00069 *           ( ( n*( n + 1 ) )/2 ).
00070 *           Before entry with  UPLO = 'U' or 'u', the array AP must
00071 *           contain the upper triangular part of the symmetric matrix
00072 *           packed sequentially, column by column, so that AP( 1 )
00073 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
00074 *           and a( 2, 2 ) respectively, and so on. On exit, the array
00075 *           AP is overwritten by the upper triangular part of the
00076 *           updated matrix.
00077 *           Before entry with UPLO = 'L' or 'l', the array AP must
00078 *           contain the lower triangular part of the symmetric matrix
00079 *           packed sequentially, column by column, so that AP( 1 )
00080 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
00081 *           and a( 3, 1 ) respectively, and so on. On exit, the array
00082 *           AP is overwritten by the lower triangular part of the
00083 *           updated matrix.
00084 *
00085 *  Further Details
00086 *  ===============
00087 *
00088 *  Level 2 Blas routine.
00089 *
00090 *  -- Written on 22-October-1986.
00091 *     Jack Dongarra, Argonne National Lab.
00092 *     Jeremy Du Croz, Nag Central Office.
00093 *     Sven Hammarling, Nag Central Office.
00094 *     Richard Hanson, Sandia National Labs.
00095 *
00096 *  =====================================================================
00097 *
00098 *     .. Parameters ..
00099       REAL ZERO
00100       PARAMETER (ZERO=0.0E+0)
00101 *     ..
00102 *     .. Local Scalars ..
00103       REAL TEMP1,TEMP2
00104       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
00105 *     ..
00106 *     .. External Functions ..
00107       LOGICAL LSAME
00108       EXTERNAL LSAME
00109 *     ..
00110 *     .. External Subroutines ..
00111       EXTERNAL XERBLA
00112 *     ..
00113 *
00114 *     Test the input parameters.
00115 *
00116       INFO = 0
00117       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
00118           INFO = 1
00119       ELSE IF (N.LT.0) THEN
00120           INFO = 2
00121       ELSE IF (INCX.EQ.0) THEN
00122           INFO = 5
00123       ELSE IF (INCY.EQ.0) THEN
00124           INFO = 7
00125       END IF
00126       IF (INFO.NE.0) THEN
00127           CALL XERBLA('SSPR2 ',INFO)
00128           RETURN
00129       END IF
00130 *
00131 *     Quick return if possible.
00132 *
00133       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
00134 *
00135 *     Set up the start points in X and Y if the increments are not both
00136 *     unity.
00137 *
00138       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
00139           IF (INCX.GT.0) THEN
00140               KX = 1
00141           ELSE
00142               KX = 1 - (N-1)*INCX
00143           END IF
00144           IF (INCY.GT.0) THEN
00145               KY = 1
00146           ELSE
00147               KY = 1 - (N-1)*INCY
00148           END IF
00149           JX = KX
00150           JY = KY
00151       END IF
00152 *
00153 *     Start the operations. In this version the elements of the array AP
00154 *     are accessed sequentially with one pass through AP.
00155 *
00156       KK = 1
00157       IF (LSAME(UPLO,'U')) THEN
00158 *
00159 *        Form  A  when upper triangle is stored in AP.
00160 *
00161           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00162               DO 20 J = 1,N
00163                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
00164                       TEMP1 = ALPHA*Y(J)
00165                       TEMP2 = ALPHA*X(J)
00166                       K = KK
00167                       DO 10 I = 1,J
00168                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
00169                           K = K + 1
00170    10                 CONTINUE
00171                   END IF
00172                   KK = KK + J
00173    20         CONTINUE
00174           ELSE
00175               DO 40 J = 1,N
00176                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
00177                       TEMP1 = ALPHA*Y(JY)
00178                       TEMP2 = ALPHA*X(JX)
00179                       IX = KX
00180                       IY = KY
00181                       DO 30 K = KK,KK + J - 1
00182                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
00183                           IX = IX + INCX
00184                           IY = IY + INCY
00185    30                 CONTINUE
00186                   END IF
00187                   JX = JX + INCX
00188                   JY = JY + INCY
00189                   KK = KK + J
00190    40         CONTINUE
00191           END IF
00192       ELSE
00193 *
00194 *        Form  A  when lower triangle is stored in AP.
00195 *
00196           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
00197               DO 60 J = 1,N
00198                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
00199                       TEMP1 = ALPHA*Y(J)
00200                       TEMP2 = ALPHA*X(J)
00201                       K = KK
00202                       DO 50 I = J,N
00203                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
00204                           K = K + 1
00205    50                 CONTINUE
00206                   END IF
00207                   KK = KK + N - J + 1
00208    60         CONTINUE
00209           ELSE
00210               DO 80 J = 1,N
00211                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
00212                       TEMP1 = ALPHA*Y(JY)
00213                       TEMP2 = ALPHA*X(JX)
00214                       IX = JX
00215                       IY = JY
00216                       DO 70 K = KK,KK + N - J
00217                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
00218                           IX = IX + INCX
00219                           IY = IY + INCY
00220    70                 CONTINUE
00221                   END IF
00222                   JX = JX + INCX
00223                   JY = JY + INCY
00224                   KK = KK + N - J + 1
00225    80         CONTINUE
00226           END IF
00227       END IF
00228 *
00229       RETURN
00230 *
00231 *     End of SSPR2 .
00232 *
00233       END