|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER INFO, JOB, N 00010 DOUBLE PRECISION TOL 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IN( * ) 00014 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLAGTS may be used to solve one of the systems of equations 00021 * 00022 * (T - lambda*I)*x = y or (T - lambda*I)'*x = y, 00023 * 00024 * where T is an n by n tridiagonal matrix, for x, following the 00025 * factorization of (T - lambda*I) as 00026 * 00027 * (T - lambda*I) = P*L*U , 00028 * 00029 * by routine DLAGTF. The choice of equation to be solved is 00030 * controlled by the argument JOB, and in each case there is an option 00031 * to perturb zero or very small diagonal elements of U, this option 00032 * being intended for use in applications such as inverse iteration. 00033 * 00034 * Arguments 00035 * ========= 00036 * 00037 * JOB (input) INTEGER 00038 * Specifies the job to be performed by DLAGTS as follows: 00039 * = 1: The equations (T - lambda*I)x = y are to be solved, 00040 * but diagonal elements of U are not to be perturbed. 00041 * = -1: The equations (T - lambda*I)x = y are to be solved 00042 * and, if overflow would otherwise occur, the diagonal 00043 * elements of U are to be perturbed. See argument TOL 00044 * below. 00045 * = 2: The equations (T - lambda*I)'x = y are to be solved, 00046 * but diagonal elements of U are not to be perturbed. 00047 * = -2: The equations (T - lambda*I)'x = y are to be solved 00048 * and, if overflow would otherwise occur, the diagonal 00049 * elements of U are to be perturbed. See argument TOL 00050 * below. 00051 * 00052 * N (input) INTEGER 00053 * The order of the matrix T. 00054 * 00055 * A (input) DOUBLE PRECISION array, dimension (N) 00056 * On entry, A must contain the diagonal elements of U as 00057 * returned from DLAGTF. 00058 * 00059 * B (input) DOUBLE PRECISION array, dimension (N-1) 00060 * On entry, B must contain the first super-diagonal elements of 00061 * U as returned from DLAGTF. 00062 * 00063 * C (input) DOUBLE PRECISION array, dimension (N-1) 00064 * On entry, C must contain the sub-diagonal elements of L as 00065 * returned from DLAGTF. 00066 * 00067 * D (input) DOUBLE PRECISION array, dimension (N-2) 00068 * On entry, D must contain the second super-diagonal elements 00069 * of U as returned from DLAGTF. 00070 * 00071 * IN (input) INTEGER array, dimension (N) 00072 * On entry, IN must contain details of the matrix P as returned 00073 * from DLAGTF. 00074 * 00075 * Y (input/output) DOUBLE PRECISION array, dimension (N) 00076 * On entry, the right hand side vector y. 00077 * On exit, Y is overwritten by the solution vector x. 00078 * 00079 * TOL (input/output) DOUBLE PRECISION 00080 * On entry, with JOB .lt. 0, TOL should be the minimum 00081 * perturbation to be made to very small diagonal elements of U. 00082 * TOL should normally be chosen as about eps*norm(U), where eps 00083 * is the relative machine precision, but if TOL is supplied as 00084 * non-positive, then it is reset to eps*max( abs( u(i,j) ) ). 00085 * If JOB .gt. 0 then TOL is not referenced. 00086 * 00087 * On exit, TOL is changed as described above, only if TOL is 00088 * non-positive on entry. Otherwise TOL is unchanged. 00089 * 00090 * INFO (output) INTEGER 00091 * = 0 : successful exit 00092 * .lt. 0: if INFO = -i, the i-th argument had an illegal value 00093 * .gt. 0: overflow would occur when computing the INFO(th) 00094 * element of the solution vector x. This can only occur 00095 * when JOB is supplied as positive and either means 00096 * that a diagonal element of U is very small, or that 00097 * the elements of the right-hand side vector y are very 00098 * large. 00099 * 00100 * ===================================================================== 00101 * 00102 * .. Parameters .. 00103 DOUBLE PRECISION ONE, ZERO 00104 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00105 * .. 00106 * .. Local Scalars .. 00107 INTEGER K 00108 DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP 00109 * .. 00110 * .. Intrinsic Functions .. 00111 INTRINSIC ABS, MAX, SIGN 00112 * .. 00113 * .. External Functions .. 00114 DOUBLE PRECISION DLAMCH 00115 EXTERNAL DLAMCH 00116 * .. 00117 * .. External Subroutines .. 00118 EXTERNAL XERBLA 00119 * .. 00120 * .. Executable Statements .. 00121 * 00122 INFO = 0 00123 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN 00124 INFO = -1 00125 ELSE IF( N.LT.0 ) THEN 00126 INFO = -2 00127 END IF 00128 IF( INFO.NE.0 ) THEN 00129 CALL XERBLA( 'DLAGTS', -INFO ) 00130 RETURN 00131 END IF 00132 * 00133 IF( N.EQ.0 ) 00134 $ RETURN 00135 * 00136 EPS = DLAMCH( 'Epsilon' ) 00137 SFMIN = DLAMCH( 'Safe minimum' ) 00138 BIGNUM = ONE / SFMIN 00139 * 00140 IF( JOB.LT.0 ) THEN 00141 IF( TOL.LE.ZERO ) THEN 00142 TOL = ABS( A( 1 ) ) 00143 IF( N.GT.1 ) 00144 $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) ) 00145 DO 10 K = 3, N 00146 TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ), 00147 $ ABS( D( K-2 ) ) ) 00148 10 CONTINUE 00149 TOL = TOL*EPS 00150 IF( TOL.EQ.ZERO ) 00151 $ TOL = EPS 00152 END IF 00153 END IF 00154 * 00155 IF( ABS( JOB ).EQ.1 ) THEN 00156 DO 20 K = 2, N 00157 IF( IN( K-1 ).EQ.0 ) THEN 00158 Y( K ) = Y( K ) - C( K-1 )*Y( K-1 ) 00159 ELSE 00160 TEMP = Y( K-1 ) 00161 Y( K-1 ) = Y( K ) 00162 Y( K ) = TEMP - C( K-1 )*Y( K ) 00163 END IF 00164 20 CONTINUE 00165 IF( JOB.EQ.1 ) THEN 00166 DO 30 K = N, 1, -1 00167 IF( K.LE.N-2 ) THEN 00168 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 00169 ELSE IF( K.EQ.N-1 ) THEN 00170 TEMP = Y( K ) - B( K )*Y( K+1 ) 00171 ELSE 00172 TEMP = Y( K ) 00173 END IF 00174 AK = A( K ) 00175 ABSAK = ABS( AK ) 00176 IF( ABSAK.LT.ONE ) THEN 00177 IF( ABSAK.LT.SFMIN ) THEN 00178 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00179 $ THEN 00180 INFO = K 00181 RETURN 00182 ELSE 00183 TEMP = TEMP*BIGNUM 00184 AK = AK*BIGNUM 00185 END IF 00186 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00187 INFO = K 00188 RETURN 00189 END IF 00190 END IF 00191 Y( K ) = TEMP / AK 00192 30 CONTINUE 00193 ELSE 00194 DO 50 K = N, 1, -1 00195 IF( K.LE.N-2 ) THEN 00196 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 00197 ELSE IF( K.EQ.N-1 ) THEN 00198 TEMP = Y( K ) - B( K )*Y( K+1 ) 00199 ELSE 00200 TEMP = Y( K ) 00201 END IF 00202 AK = A( K ) 00203 PERT = SIGN( TOL, AK ) 00204 40 CONTINUE 00205 ABSAK = ABS( AK ) 00206 IF( ABSAK.LT.ONE ) THEN 00207 IF( ABSAK.LT.SFMIN ) THEN 00208 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00209 $ THEN 00210 AK = AK + PERT 00211 PERT = 2*PERT 00212 GO TO 40 00213 ELSE 00214 TEMP = TEMP*BIGNUM 00215 AK = AK*BIGNUM 00216 END IF 00217 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00218 AK = AK + PERT 00219 PERT = 2*PERT 00220 GO TO 40 00221 END IF 00222 END IF 00223 Y( K ) = TEMP / AK 00224 50 CONTINUE 00225 END IF 00226 ELSE 00227 * 00228 * Come to here if JOB = 2 or -2 00229 * 00230 IF( JOB.EQ.2 ) THEN 00231 DO 60 K = 1, N 00232 IF( K.GE.3 ) THEN 00233 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 00234 ELSE IF( K.EQ.2 ) THEN 00235 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 00236 ELSE 00237 TEMP = Y( K ) 00238 END IF 00239 AK = A( K ) 00240 ABSAK = ABS( AK ) 00241 IF( ABSAK.LT.ONE ) THEN 00242 IF( ABSAK.LT.SFMIN ) THEN 00243 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00244 $ THEN 00245 INFO = K 00246 RETURN 00247 ELSE 00248 TEMP = TEMP*BIGNUM 00249 AK = AK*BIGNUM 00250 END IF 00251 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00252 INFO = K 00253 RETURN 00254 END IF 00255 END IF 00256 Y( K ) = TEMP / AK 00257 60 CONTINUE 00258 ELSE 00259 DO 80 K = 1, N 00260 IF( K.GE.3 ) THEN 00261 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 00262 ELSE IF( K.EQ.2 ) THEN 00263 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 00264 ELSE 00265 TEMP = Y( K ) 00266 END IF 00267 AK = A( K ) 00268 PERT = SIGN( TOL, AK ) 00269 70 CONTINUE 00270 ABSAK = ABS( AK ) 00271 IF( ABSAK.LT.ONE ) THEN 00272 IF( ABSAK.LT.SFMIN ) THEN 00273 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00274 $ THEN 00275 AK = AK + PERT 00276 PERT = 2*PERT 00277 GO TO 70 00278 ELSE 00279 TEMP = TEMP*BIGNUM 00280 AK = AK*BIGNUM 00281 END IF 00282 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00283 AK = AK + PERT 00284 PERT = 2*PERT 00285 GO TO 70 00286 END IF 00287 END IF 00288 Y( K ) = TEMP / AK 00289 80 CONTINUE 00290 END IF 00291 * 00292 DO 90 K = N, 2, -1 00293 IF( IN( K-1 ).EQ.0 ) THEN 00294 Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K ) 00295 ELSE 00296 TEMP = Y( K-1 ) 00297 Y( K-1 ) = Y( K ) 00298 Y( K ) = TEMP - C( K-1 )*Y( K ) 00299 END IF 00300 90 CONTINUE 00301 END IF 00302 * 00303 * End of DLAGTS 00304 * 00305 END
1.7.3 
