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