Annotation of rpl/lapack/blas/dspr2.f, revision 1.16

1.8       bertrand    1: *> \brief \b DSPR2
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.13      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *  Definition:
                      9: *  ===========
                     10: *
                     11: *       SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
1.13      bertrand   12: *
1.8       bertrand   13: *       .. Scalar Arguments ..
                     14: *       DOUBLE PRECISION ALPHA
                     15: *       INTEGER INCX,INCY,N
                     16: *       CHARACTER UPLO
                     17: *       ..
                     18: *       .. Array Arguments ..
                     19: *       DOUBLE PRECISION AP(*),X(*),Y(*)
                     20: *       ..
1.13      bertrand   21: *
1.8       bertrand   22: *
                     23: *> \par Purpose:
                     24: *  =============
                     25: *>
                     26: *> \verbatim
                     27: *>
                     28: *> DSPR2  performs the symmetric rank 2 operation
                     29: *>
                     30: *>    A := alpha*x*y**T + alpha*y*x**T + A,
                     31: *>
                     32: *> where alpha is a scalar, x and y are n element vectors and A is an
                     33: *> n by n symmetric matrix, supplied in packed form.
                     34: *> \endverbatim
                     35: *
                     36: *  Arguments:
                     37: *  ==========
                     38: *
                     39: *> \param[in] UPLO
                     40: *> \verbatim
                     41: *>          UPLO is CHARACTER*1
                     42: *>           On entry, UPLO specifies whether the upper or lower
                     43: *>           triangular part of the matrix A is supplied in the packed
                     44: *>           array AP as follows:
                     45: *>
                     46: *>              UPLO = 'U' or 'u'   The upper triangular part of A is
                     47: *>                                  supplied in AP.
                     48: *>
                     49: *>              UPLO = 'L' or 'l'   The lower triangular part of A is
                     50: *>                                  supplied in AP.
                     51: *> \endverbatim
                     52: *>
                     53: *> \param[in] N
                     54: *> \verbatim
                     55: *>          N is INTEGER
                     56: *>           On entry, N specifies the order of the matrix A.
                     57: *>           N must be at least zero.
                     58: *> \endverbatim
                     59: *>
                     60: *> \param[in] ALPHA
                     61: *> \verbatim
                     62: *>          ALPHA is DOUBLE PRECISION.
                     63: *>           On entry, ALPHA specifies the scalar alpha.
                     64: *> \endverbatim
                     65: *>
                     66: *> \param[in] X
                     67: *> \verbatim
1.14      bertrand   68: *>          X is DOUBLE PRECISION array, dimension at least
1.8       bertrand   69: *>           ( 1 + ( n - 1 )*abs( INCX ) ).
                     70: *>           Before entry, the incremented array X must contain the n
                     71: *>           element vector x.
                     72: *> \endverbatim
                     73: *>
                     74: *> \param[in] INCX
                     75: *> \verbatim
                     76: *>          INCX is INTEGER
                     77: *>           On entry, INCX specifies the increment for the elements of
                     78: *>           X. INCX must not be zero.
                     79: *> \endverbatim
                     80: *>
                     81: *> \param[in] Y
                     82: *> \verbatim
1.14      bertrand   83: *>          Y is DOUBLE PRECISION array, dimension at least
1.8       bertrand   84: *>           ( 1 + ( n - 1 )*abs( INCY ) ).
                     85: *>           Before entry, the incremented array Y must contain the n
                     86: *>           element vector y.
                     87: *> \endverbatim
                     88: *>
                     89: *> \param[in] INCY
                     90: *> \verbatim
                     91: *>          INCY is INTEGER
                     92: *>           On entry, INCY specifies the increment for the elements of
                     93: *>           Y. INCY must not be zero.
                     94: *> \endverbatim
                     95: *>
                     96: *> \param[in,out] AP
                     97: *> \verbatim
1.14      bertrand   98: *>          AP is DOUBLE PRECISION array, dimension at least
1.8       bertrand   99: *>           ( ( n*( n + 1 ) )/2 ).
                    100: *>           Before entry with  UPLO = 'U' or 'u', the array AP must
                    101: *>           contain the upper triangular part of the symmetric matrix
                    102: *>           packed sequentially, column by column, so that AP( 1 )
                    103: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
                    104: *>           and a( 2, 2 ) respectively, and so on. On exit, the array
                    105: *>           AP is overwritten by the upper triangular part of the
                    106: *>           updated matrix.
                    107: *>           Before entry with UPLO = 'L' or 'l', the array AP must
                    108: *>           contain the lower triangular part of the symmetric matrix
                    109: *>           packed sequentially, column by column, so that AP( 1 )
                    110: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
                    111: *>           and a( 3, 1 ) respectively, and so on. On exit, the array
                    112: *>           AP is overwritten by the lower triangular part of the
                    113: *>           updated matrix.
                    114: *> \endverbatim
                    115: *
                    116: *  Authors:
                    117: *  ========
                    118: *
1.13      bertrand  119: *> \author Univ. of Tennessee
                    120: *> \author Univ. of California Berkeley
                    121: *> \author Univ. of Colorado Denver
                    122: *> \author NAG Ltd.
1.8       bertrand  123: *
                    124: *> \ingroup double_blas_level2
                    125: *
                    126: *> \par Further Details:
                    127: *  =====================
                    128: *>
                    129: *> \verbatim
                    130: *>
                    131: *>  Level 2 Blas routine.
                    132: *>
                    133: *>  -- Written on 22-October-1986.
                    134: *>     Jack Dongarra, Argonne National Lab.
                    135: *>     Jeremy Du Croz, Nag Central Office.
                    136: *>     Sven Hammarling, Nag Central Office.
                    137: *>     Richard Hanson, Sandia National Labs.
                    138: *> \endverbatim
                    139: *>
                    140: *  =====================================================================
1.1       bertrand  141:       SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
1.8       bertrand  142: *
1.16    ! bertrand  143: *  -- Reference BLAS level2 routine --
1.8       bertrand  144: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
                    145: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    146: *
1.1       bertrand  147: *     .. Scalar Arguments ..
                    148:       DOUBLE PRECISION ALPHA
                    149:       INTEGER INCX,INCY,N
                    150:       CHARACTER UPLO
                    151: *     ..
                    152: *     .. Array Arguments ..
                    153:       DOUBLE PRECISION AP(*),X(*),Y(*)
                    154: *     ..
                    155: *
                    156: *  =====================================================================
                    157: *
                    158: *     .. Parameters ..
                    159:       DOUBLE PRECISION ZERO
                    160:       PARAMETER (ZERO=0.0D+0)
                    161: *     ..
                    162: *     .. Local Scalars ..
                    163:       DOUBLE PRECISION TEMP1,TEMP2
                    164:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
                    165: *     ..
                    166: *     .. External Functions ..
                    167:       LOGICAL LSAME
                    168:       EXTERNAL LSAME
                    169: *     ..
                    170: *     .. External Subroutines ..
                    171:       EXTERNAL XERBLA
                    172: *     ..
                    173: *
                    174: *     Test the input parameters.
                    175: *
                    176:       INFO = 0
                    177:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
                    178:           INFO = 1
                    179:       ELSE IF (N.LT.0) THEN
                    180:           INFO = 2
                    181:       ELSE IF (INCX.EQ.0) THEN
                    182:           INFO = 5
                    183:       ELSE IF (INCY.EQ.0) THEN
                    184:           INFO = 7
                    185:       END IF
                    186:       IF (INFO.NE.0) THEN
                    187:           CALL XERBLA('DSPR2 ',INFO)
                    188:           RETURN
                    189:       END IF
                    190: *
                    191: *     Quick return if possible.
                    192: *
                    193:       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
                    194: *
                    195: *     Set up the start points in X and Y if the increments are not both
                    196: *     unity.
                    197: *
                    198:       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
                    199:           IF (INCX.GT.0) THEN
                    200:               KX = 1
                    201:           ELSE
                    202:               KX = 1 - (N-1)*INCX
                    203:           END IF
                    204:           IF (INCY.GT.0) THEN
                    205:               KY = 1
                    206:           ELSE
                    207:               KY = 1 - (N-1)*INCY
                    208:           END IF
                    209:           JX = KX
                    210:           JY = KY
                    211:       END IF
                    212: *
                    213: *     Start the operations. In this version the elements of the array AP
                    214: *     are accessed sequentially with one pass through AP.
                    215: *
                    216:       KK = 1
                    217:       IF (LSAME(UPLO,'U')) THEN
                    218: *
                    219: *        Form  A  when upper triangle is stored in AP.
                    220: *
                    221:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
                    222:               DO 20 J = 1,N
                    223:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
                    224:                       TEMP1 = ALPHA*Y(J)
                    225:                       TEMP2 = ALPHA*X(J)
                    226:                       K = KK
                    227:                       DO 10 I = 1,J
                    228:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
                    229:                           K = K + 1
                    230:    10                 CONTINUE
                    231:                   END IF
                    232:                   KK = KK + J
                    233:    20         CONTINUE
                    234:           ELSE
                    235:               DO 40 J = 1,N
                    236:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
                    237:                       TEMP1 = ALPHA*Y(JY)
                    238:                       TEMP2 = ALPHA*X(JX)
                    239:                       IX = KX
                    240:                       IY = KY
                    241:                       DO 30 K = KK,KK + J - 1
                    242:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
                    243:                           IX = IX + INCX
                    244:                           IY = IY + INCY
                    245:    30                 CONTINUE
                    246:                   END IF
                    247:                   JX = JX + INCX
                    248:                   JY = JY + INCY
                    249:                   KK = KK + J
                    250:    40         CONTINUE
                    251:           END IF
                    252:       ELSE
                    253: *
                    254: *        Form  A  when lower triangle is stored in AP.
                    255: *
                    256:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
                    257:               DO 60 J = 1,N
                    258:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
                    259:                       TEMP1 = ALPHA*Y(J)
                    260:                       TEMP2 = ALPHA*X(J)
                    261:                       K = KK
                    262:                       DO 50 I = J,N
                    263:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
                    264:                           K = K + 1
                    265:    50                 CONTINUE
                    266:                   END IF
                    267:                   KK = KK + N - J + 1
                    268:    60         CONTINUE
                    269:           ELSE
                    270:               DO 80 J = 1,N
                    271:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
                    272:                       TEMP1 = ALPHA*Y(JY)
                    273:                       TEMP2 = ALPHA*X(JX)
                    274:                       IX = JX
                    275:                       IY = JY
                    276:                       DO 70 K = KK,KK + N - J
                    277:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
                    278:                           IX = IX + INCX
                    279:                           IY = IY + INCY
                    280:    70                 CONTINUE
                    281:                   END IF
                    282:                   JX = JX + INCX
                    283:                   JY = JY + INCY
                    284:                   KK = KK + N - J + 1
                    285:    80         CONTINUE
                    286:           END IF
                    287:       END IF
                    288: *
                    289:       RETURN
                    290: *
1.16    ! bertrand  291: *     End of DSPR2
1.1       bertrand  292: *
                    293:       END

CVSweb interface <joel.bertrand@systella.fr>