![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAGTS 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAGTS + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlagts.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlagts.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlagts.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, JOB, N 00025 * DOUBLE PRECISION TOL 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IN( * ) 00029 * DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DLAGTS may be used to solve one of the systems of equations 00039 *> 00040 *> (T - lambda*I)*x = y or (T - lambda*I)**T*x = y, 00041 *> 00042 *> where T is an n by n tridiagonal matrix, for x, following the 00043 *> factorization of (T - lambda*I) as 00044 *> 00045 *> (T - lambda*I) = P*L*U , 00046 *> 00047 *> by routine DLAGTF. The choice of equation to be solved is 00048 *> controlled by the argument JOB, and in each case there is an option 00049 *> to perturb zero or very small diagonal elements of U, this option 00050 *> being intended for use in applications such as inverse iteration. 00051 *> \endverbatim 00052 * 00053 * Arguments: 00054 * ========== 00055 * 00056 *> \param[in] JOB 00057 *> \verbatim 00058 *> JOB is INTEGER 00059 *> Specifies the job to be performed by DLAGTS as follows: 00060 *> = 1: The equations (T - lambda*I)x = y are to be solved, 00061 *> but diagonal elements of U are not to be perturbed. 00062 *> = -1: The equations (T - lambda*I)x = y are to be solved 00063 *> and, if overflow would otherwise occur, the diagonal 00064 *> elements of U are to be perturbed. See argument TOL 00065 *> below. 00066 *> = 2: The equations (T - lambda*I)**Tx = y are to be solved, 00067 *> but diagonal elements of U are not to be perturbed. 00068 *> = -2: The equations (T - lambda*I)**Tx = y are to be solved 00069 *> and, if overflow would otherwise occur, the diagonal 00070 *> elements of U are to be perturbed. See argument TOL 00071 *> below. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] N 00075 *> \verbatim 00076 *> N is INTEGER 00077 *> The order of the matrix T. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] A 00081 *> \verbatim 00082 *> A is DOUBLE PRECISION array, dimension (N) 00083 *> On entry, A must contain the diagonal elements of U as 00084 *> returned from DLAGTF. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] B 00088 *> \verbatim 00089 *> B is DOUBLE PRECISION array, dimension (N-1) 00090 *> On entry, B must contain the first super-diagonal elements of 00091 *> U as returned from DLAGTF. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] C 00095 *> \verbatim 00096 *> C is DOUBLE PRECISION array, dimension (N-1) 00097 *> On entry, C must contain the sub-diagonal elements of L as 00098 *> returned from DLAGTF. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] D 00102 *> \verbatim 00103 *> D is DOUBLE PRECISION array, dimension (N-2) 00104 *> On entry, D must contain the second super-diagonal elements 00105 *> of U as returned from DLAGTF. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] IN 00109 *> \verbatim 00110 *> IN is INTEGER array, dimension (N) 00111 *> On entry, IN must contain details of the matrix P as returned 00112 *> from DLAGTF. 00113 *> \endverbatim 00114 *> 00115 *> \param[in,out] Y 00116 *> \verbatim 00117 *> Y is DOUBLE PRECISION array, dimension (N) 00118 *> On entry, the right hand side vector y. 00119 *> On exit, Y is overwritten by the solution vector x. 00120 *> \endverbatim 00121 *> 00122 *> \param[in,out] TOL 00123 *> \verbatim 00124 *> TOL is DOUBLE PRECISION 00125 *> On entry, with JOB .lt. 0, TOL should be the minimum 00126 *> perturbation to be made to very small diagonal elements of U. 00127 *> TOL should normally be chosen as about eps*norm(U), where eps 00128 *> is the relative machine precision, but if TOL is supplied as 00129 *> non-positive, then it is reset to eps*max( abs( u(i,j) ) ). 00130 *> If JOB .gt. 0 then TOL is not referenced. 00131 *> 00132 *> On exit, TOL is changed as described above, only if TOL is 00133 *> non-positive on entry. Otherwise TOL is unchanged. 00134 *> \endverbatim 00135 *> 00136 *> \param[out] INFO 00137 *> \verbatim 00138 *> INFO is INTEGER 00139 *> = 0 : successful exit 00140 *> .lt. 0: if INFO = -i, the i-th argument had an illegal value 00141 *> .gt. 0: overflow would occur when computing the INFO(th) 00142 *> element of the solution vector x. This can only occur 00143 *> when JOB is supplied as positive and either means 00144 *> that a diagonal element of U is very small, or that 00145 *> the elements of the right-hand side vector y are very 00146 *> large. 00147 *> \endverbatim 00148 * 00149 * Authors: 00150 * ======== 00151 * 00152 *> \author Univ. of Tennessee 00153 *> \author Univ. of California Berkeley 00154 *> \author Univ. of Colorado Denver 00155 *> \author NAG Ltd. 00156 * 00157 *> \date November 2011 00158 * 00159 *> \ingroup auxOTHERauxiliary 00160 * 00161 * ===================================================================== 00162 SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO ) 00163 * 00164 * -- LAPACK auxiliary routine (version 3.4.0) -- 00165 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00167 * November 2011 00168 * 00169 * .. Scalar Arguments .. 00170 INTEGER INFO, JOB, N 00171 DOUBLE PRECISION TOL 00172 * .. 00173 * .. Array Arguments .. 00174 INTEGER IN( * ) 00175 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * ) 00176 * .. 00177 * 00178 * ===================================================================== 00179 * 00180 * .. Parameters .. 00181 DOUBLE PRECISION ONE, ZERO 00182 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00183 * .. 00184 * .. Local Scalars .. 00185 INTEGER K 00186 DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP 00187 * .. 00188 * .. Intrinsic Functions .. 00189 INTRINSIC ABS, MAX, SIGN 00190 * .. 00191 * .. External Functions .. 00192 DOUBLE PRECISION DLAMCH 00193 EXTERNAL DLAMCH 00194 * .. 00195 * .. External Subroutines .. 00196 EXTERNAL XERBLA 00197 * .. 00198 * .. Executable Statements .. 00199 * 00200 INFO = 0 00201 IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN 00202 INFO = -1 00203 ELSE IF( N.LT.0 ) THEN 00204 INFO = -2 00205 END IF 00206 IF( INFO.NE.0 ) THEN 00207 CALL XERBLA( 'DLAGTS', -INFO ) 00208 RETURN 00209 END IF 00210 * 00211 IF( N.EQ.0 ) 00212 $ RETURN 00213 * 00214 EPS = DLAMCH( 'Epsilon' ) 00215 SFMIN = DLAMCH( 'Safe minimum' ) 00216 BIGNUM = ONE / SFMIN 00217 * 00218 IF( JOB.LT.0 ) THEN 00219 IF( TOL.LE.ZERO ) THEN 00220 TOL = ABS( A( 1 ) ) 00221 IF( N.GT.1 ) 00222 $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) ) 00223 DO 10 K = 3, N 00224 TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ), 00225 $ ABS( D( K-2 ) ) ) 00226 10 CONTINUE 00227 TOL = TOL*EPS 00228 IF( TOL.EQ.ZERO ) 00229 $ TOL = EPS 00230 END IF 00231 END IF 00232 * 00233 IF( ABS( JOB ).EQ.1 ) THEN 00234 DO 20 K = 2, N 00235 IF( IN( K-1 ).EQ.0 ) THEN 00236 Y( K ) = Y( K ) - C( K-1 )*Y( K-1 ) 00237 ELSE 00238 TEMP = Y( K-1 ) 00239 Y( K-1 ) = Y( K ) 00240 Y( K ) = TEMP - C( K-1 )*Y( K ) 00241 END IF 00242 20 CONTINUE 00243 IF( JOB.EQ.1 ) THEN 00244 DO 30 K = N, 1, -1 00245 IF( K.LE.N-2 ) THEN 00246 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 00247 ELSE IF( K.EQ.N-1 ) THEN 00248 TEMP = Y( K ) - B( K )*Y( K+1 ) 00249 ELSE 00250 TEMP = Y( K ) 00251 END IF 00252 AK = A( K ) 00253 ABSAK = ABS( AK ) 00254 IF( ABSAK.LT.ONE ) THEN 00255 IF( ABSAK.LT.SFMIN ) THEN 00256 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00257 $ THEN 00258 INFO = K 00259 RETURN 00260 ELSE 00261 TEMP = TEMP*BIGNUM 00262 AK = AK*BIGNUM 00263 END IF 00264 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00265 INFO = K 00266 RETURN 00267 END IF 00268 END IF 00269 Y( K ) = TEMP / AK 00270 30 CONTINUE 00271 ELSE 00272 DO 50 K = N, 1, -1 00273 IF( K.LE.N-2 ) THEN 00274 TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 ) 00275 ELSE IF( K.EQ.N-1 ) THEN 00276 TEMP = Y( K ) - B( K )*Y( K+1 ) 00277 ELSE 00278 TEMP = Y( K ) 00279 END IF 00280 AK = A( K ) 00281 PERT = SIGN( TOL, AK ) 00282 40 CONTINUE 00283 ABSAK = ABS( AK ) 00284 IF( ABSAK.LT.ONE ) THEN 00285 IF( ABSAK.LT.SFMIN ) THEN 00286 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00287 $ THEN 00288 AK = AK + PERT 00289 PERT = 2*PERT 00290 GO TO 40 00291 ELSE 00292 TEMP = TEMP*BIGNUM 00293 AK = AK*BIGNUM 00294 END IF 00295 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00296 AK = AK + PERT 00297 PERT = 2*PERT 00298 GO TO 40 00299 END IF 00300 END IF 00301 Y( K ) = TEMP / AK 00302 50 CONTINUE 00303 END IF 00304 ELSE 00305 * 00306 * Come to here if JOB = 2 or -2 00307 * 00308 IF( JOB.EQ.2 ) THEN 00309 DO 60 K = 1, N 00310 IF( K.GE.3 ) THEN 00311 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 00312 ELSE IF( K.EQ.2 ) THEN 00313 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 00314 ELSE 00315 TEMP = Y( K ) 00316 END IF 00317 AK = A( K ) 00318 ABSAK = ABS( AK ) 00319 IF( ABSAK.LT.ONE ) THEN 00320 IF( ABSAK.LT.SFMIN ) THEN 00321 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00322 $ THEN 00323 INFO = K 00324 RETURN 00325 ELSE 00326 TEMP = TEMP*BIGNUM 00327 AK = AK*BIGNUM 00328 END IF 00329 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00330 INFO = K 00331 RETURN 00332 END IF 00333 END IF 00334 Y( K ) = TEMP / AK 00335 60 CONTINUE 00336 ELSE 00337 DO 80 K = 1, N 00338 IF( K.GE.3 ) THEN 00339 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 ) 00340 ELSE IF( K.EQ.2 ) THEN 00341 TEMP = Y( K ) - B( K-1 )*Y( K-1 ) 00342 ELSE 00343 TEMP = Y( K ) 00344 END IF 00345 AK = A( K ) 00346 PERT = SIGN( TOL, AK ) 00347 70 CONTINUE 00348 ABSAK = ABS( AK ) 00349 IF( ABSAK.LT.ONE ) THEN 00350 IF( ABSAK.LT.SFMIN ) THEN 00351 IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK ) 00352 $ THEN 00353 AK = AK + PERT 00354 PERT = 2*PERT 00355 GO TO 70 00356 ELSE 00357 TEMP = TEMP*BIGNUM 00358 AK = AK*BIGNUM 00359 END IF 00360 ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN 00361 AK = AK + PERT 00362 PERT = 2*PERT 00363 GO TO 70 00364 END IF 00365 END IF 00366 Y( K ) = TEMP / AK 00367 80 CONTINUE 00368 END IF 00369 * 00370 DO 90 K = N, 2, -1 00371 IF( IN( K-1 ).EQ.0 ) THEN 00372 Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K ) 00373 ELSE 00374 TEMP = Y( K-1 ) 00375 Y( K-1 ) = Y( K ) 00376 Y( K ) = TEMP - C( K-1 )*Y( K ) 00377 END IF 00378 90 CONTINUE 00379 END IF 00380 * 00381 * End of DLAGTS 00382 * 00383 END