 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZSBMV( UPLO, N, K, ALPHA, A, LDA, X, INCX, BETA, Y, 00002 $ INCY ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.1) -- 00005 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 CHARACTER UPLO 00010 INTEGER INCX, INCY, K, LDA, N 00011 COMPLEX*16 ALPHA, BETA 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX*16 A( LDA, * ), X( * ), Y( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * ZSBMV performs the matrix-vector operation 00021 * 00022 * y := alpha*A*x + beta*y, 00023 * 00024 * where alpha and beta are scalars, x and y are n element vectors and 00025 * A is an n by n symmetric band matrix, with k super-diagonals. 00026 * 00027 * Arguments 00028 * ========== 00029 * 00030 * UPLO - CHARACTER*1 00031 * On entry, UPLO specifies whether the upper or lower 00032 * triangular part of the band matrix A is being supplied as 00033 * follows: 00034 * 00035 * UPLO = 'U' or 'u' The upper triangular part of A is 00036 * being supplied. 00037 * 00038 * UPLO = 'L' or 'l' The lower triangular part of A is 00039 * being supplied. 00040 * 00041 * Unchanged on exit. 00042 * 00043 * N - INTEGER 00044 * On entry, N specifies the order of the matrix A. 00045 * N must be at least zero. 00046 * Unchanged on exit. 00047 * 00048 * K - INTEGER 00049 * On entry, K specifies the number of super-diagonals of the 00050 * matrix A. K must satisfy 0 .le. K. 00051 * Unchanged on exit. 00052 * 00053 * ALPHA - COMPLEX*16 00054 * On entry, ALPHA specifies the scalar alpha. 00055 * Unchanged on exit. 00056 * 00057 * A - COMPLEX*16 array, dimension( LDA, N ) 00058 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) 00059 * by n part of the array A must contain the upper triangular 00060 * band part of the symmetric matrix, supplied column by 00061 * column, with the leading diagonal of the matrix in row 00062 * ( k + 1 ) of the array, the first super-diagonal starting at 00063 * position 2 in row k, and so on. The top left k by k triangle 00064 * of the array A is not referenced. 00065 * The following program segment will transfer the upper 00066 * triangular part of a symmetric band matrix from conventional 00067 * full matrix storage to band storage: 00068 * 00069 * DO 20, J = 1, N 00070 * M = K + 1 - J 00071 * DO 10, I = MAX( 1, J - K ), J 00072 * A( M + I, J ) = matrix( I, J ) 00073 * 10 CONTINUE 00074 * 20 CONTINUE 00075 * 00076 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) 00077 * by n part of the array A must contain the lower triangular 00078 * band part of the symmetric matrix, supplied column by 00079 * column, with the leading diagonal of the matrix in row 1 of 00080 * the array, the first sub-diagonal starting at position 1 in 00081 * row 2, and so on. The bottom right k by k triangle of the 00082 * array A is not referenced. 00083 * The following program segment will transfer the lower 00084 * triangular part of a symmetric band matrix from conventional 00085 * full matrix storage to band storage: 00086 * 00087 * DO 20, J = 1, N 00088 * M = 1 - J 00089 * DO 10, I = J, MIN( N, J + K ) 00090 * A( M + I, J ) = matrix( I, J ) 00091 * 10 CONTINUE 00092 * 20 CONTINUE 00093 * 00094 * Unchanged on exit. 00095 * 00096 * LDA - INTEGER 00097 * On entry, LDA specifies the first dimension of A as declared 00098 * in the calling (sub) program. LDA must be at least 00099 * ( k + 1 ). 00100 * Unchanged on exit. 00101 * 00102 * X - COMPLEX*16 array, dimension at least 00103 * ( 1 + ( N - 1 )*abs( INCX ) ). 00104 * Before entry, the incremented array X must contain the 00105 * vector x. 00106 * Unchanged on exit. 00107 * 00108 * INCX - INTEGER 00109 * On entry, INCX specifies the increment for the elements of 00110 * X. INCX must not be zero. 00111 * Unchanged on exit. 00112 * 00113 * BETA - COMPLEX*16 00114 * On entry, BETA specifies the scalar beta. 00115 * Unchanged on exit. 00116 * 00117 * Y - COMPLEX*16 array, dimension at least 00118 * ( 1 + ( N - 1 )*abs( INCY ) ). 00119 * Before entry, the incremented array Y must contain the 00120 * vector y. On exit, Y is overwritten by the updated vector y. 00121 * 00122 * INCY - INTEGER 00123 * On entry, INCY specifies the increment for the elements of 00124 * Y. INCY must not be zero. 00125 * Unchanged on exit. 00126 * 00127 * ===================================================================== 00128 * 00129 * .. Parameters .. 00130 COMPLEX*16 ONE 00131 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) 00132 COMPLEX*16 ZERO 00133 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 00134 * .. 00135 * .. Local Scalars .. 00136 INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L 00137 COMPLEX*16 TEMP1, TEMP2 00138 * .. 00139 * .. External Functions .. 00140 LOGICAL LSAME 00141 EXTERNAL LSAME 00142 * .. 00143 * .. External Subroutines .. 00144 EXTERNAL XERBLA 00145 * .. 00146 * .. Intrinsic Functions .. 00147 INTRINSIC MAX, MIN 00148 * .. 00149 * .. Executable Statements .. 00150 * 00151 * Test the input parameters. 00152 * 00153 INFO = 0 00154 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00155 INFO = 1 00156 ELSE IF( N.LT.0 ) THEN 00157 INFO = 2 00158 ELSE IF( K.LT.0 ) THEN 00159 INFO = 3 00160 ELSE IF( LDA.LT.( K+1 ) ) THEN 00161 INFO = 6 00162 ELSE IF( INCX.EQ.0 ) THEN 00163 INFO = 8 00164 ELSE IF( INCY.EQ.0 ) THEN 00165 INFO = 11 00166 END IF 00167 IF( INFO.NE.0 ) THEN 00168 CALL XERBLA( 'ZSBMV ', INFO ) 00169 RETURN 00170 END IF 00171 * 00172 * Quick return if possible. 00173 * 00174 IF( ( N.EQ.0 ) .OR. ( ( ALPHA.EQ.ZERO ) .AND. ( BETA.EQ.ONE ) ) ) 00175 $ RETURN 00176 * 00177 * Set up the start points in X and Y. 00178 * 00179 IF( INCX.GT.0 ) THEN 00180 KX = 1 00181 ELSE 00182 KX = 1 - ( N-1 )*INCX 00183 END IF 00184 IF( INCY.GT.0 ) THEN 00185 KY = 1 00186 ELSE 00187 KY = 1 - ( N-1 )*INCY 00188 END IF 00189 * 00190 * Start the operations. In this version the elements of the array A 00191 * are accessed sequentially with one pass through A. 00192 * 00193 * First form y := beta*y. 00194 * 00195 IF( BETA.NE.ONE ) THEN 00196 IF( INCY.EQ.1 ) THEN 00197 IF( BETA.EQ.ZERO ) THEN 00198 DO 10 I = 1, N 00199 Y( I ) = ZERO 00200 10 CONTINUE 00201 ELSE 00202 DO 20 I = 1, N 00203 Y( I ) = BETA*Y( I ) 00204 20 CONTINUE 00205 END IF 00206 ELSE 00207 IY = KY 00208 IF( BETA.EQ.ZERO ) THEN 00209 DO 30 I = 1, N 00210 Y( IY ) = ZERO 00211 IY = IY + INCY 00212 30 CONTINUE 00213 ELSE 00214 DO 40 I = 1, N 00215 Y( IY ) = BETA*Y( IY ) 00216 IY = IY + INCY 00217 40 CONTINUE 00218 END IF 00219 END IF 00220 END IF 00221 IF( ALPHA.EQ.ZERO ) 00222 $ RETURN 00223 IF( LSAME( UPLO, 'U' ) ) THEN 00224 * 00225 * Form y when upper triangle of A is stored. 00226 * 00227 KPLUS1 = K + 1 00228 IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN 00229 DO 60 J = 1, N 00230 TEMP1 = ALPHA*X( J ) 00231 TEMP2 = ZERO 00232 L = KPLUS1 - J 00233 DO 50 I = MAX( 1, J-K ), J - 1 00234 Y( I ) = Y( I ) + TEMP1*A( L+I, J ) 00235 TEMP2 = TEMP2 + A( L+I, J )*X( I ) 00236 50 CONTINUE 00237 Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 00238 60 CONTINUE 00239 ELSE 00240 JX = KX 00241 JY = KY 00242 DO 80 J = 1, N 00243 TEMP1 = ALPHA*X( JX ) 00244 TEMP2 = ZERO 00245 IX = KX 00246 IY = KY 00247 L = KPLUS1 - J 00248 DO 70 I = MAX( 1, J-K ), J - 1 00249 Y( IY ) = Y( IY ) + TEMP1*A( L+I, J ) 00250 TEMP2 = TEMP2 + A( L+I, J )*X( IX ) 00251 IX = IX + INCX 00252 IY = IY + INCY 00253 70 CONTINUE 00254 Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 00255 JX = JX + INCX 00256 JY = JY + INCY 00257 IF( J.GT.K ) THEN 00258 KX = KX + INCX 00259 KY = KY + INCY 00260 END IF 00261 80 CONTINUE 00262 END IF 00263 ELSE 00264 * 00265 * Form y when lower triangle of A is stored. 00266 * 00267 IF( ( INCX.EQ.1 ) .AND. ( INCY.EQ.1 ) ) THEN 00268 DO 100 J = 1, N 00269 TEMP1 = ALPHA*X( J ) 00270 TEMP2 = ZERO 00271 Y( J ) = Y( J ) + TEMP1*A( 1, J ) 00272 L = 1 - J 00273 DO 90 I = J + 1, MIN( N, J+K ) 00274 Y( I ) = Y( I ) + TEMP1*A( L+I, J ) 00275 TEMP2 = TEMP2 + A( L+I, J )*X( I ) 00276 90 CONTINUE 00277 Y( J ) = Y( J ) + ALPHA*TEMP2 00278 100 CONTINUE 00279 ELSE 00280 JX = KX 00281 JY = KY 00282 DO 120 J = 1, N 00283 TEMP1 = ALPHA*X( JX ) 00284 TEMP2 = ZERO 00285 Y( JY ) = Y( JY ) + TEMP1*A( 1, J ) 00286 L = 1 - J 00287 IX = JX 00288 IY = JY 00289 DO 110 I = J + 1, MIN( N, J+K ) 00290 IX = IX + INCX 00291 IY = IY + INCY 00292 Y( IY ) = Y( IY ) + TEMP1*A( L+I, J ) 00293 TEMP2 = TEMP2 + A( L+I, J )*X( IX ) 00294 110 CONTINUE 00295 Y( JY ) = Y( JY ) + ALPHA*TEMP2 00296 JX = JX + INCX 00297 JY = JY + INCY 00298 120 CONTINUE 00299 END IF 00300 END IF 00301 * 00302 RETURN 00303 * 00304 * End of ZSBMV 00305 * 00306 END
1.7.3 
