Annotation of rpl/lapack/lapack/dlagts.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAGTS( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
! 2: *
! 3: * -- LAPACK auxiliary routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: INTEGER INFO, JOB, N
! 10: DOUBLE PRECISION TOL
! 11: * ..
! 12: * .. Array Arguments ..
! 13: INTEGER IN( * )
! 14: DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
! 15: * ..
! 16: *
! 17: * Purpose
! 18: * =======
! 19: *
! 20: * DLAGTS may be used to solve one of the systems of equations
! 21: *
! 22: * (T - lambda*I)*x = y or (T - lambda*I)'*x = y,
! 23: *
! 24: * where T is an n by n tridiagonal matrix, for x, following the
! 25: * factorization of (T - lambda*I) as
! 26: *
! 27: * (T - lambda*I) = P*L*U ,
! 28: *
! 29: * by routine DLAGTF. The choice of equation to be solved is
! 30: * controlled by the argument JOB, and in each case there is an option
! 31: * to perturb zero or very small diagonal elements of U, this option
! 32: * being intended for use in applications such as inverse iteration.
! 33: *
! 34: * Arguments
! 35: * =========
! 36: *
! 37: * JOB (input) INTEGER
! 38: * Specifies the job to be performed by DLAGTS as follows:
! 39: * = 1: The equations (T - lambda*I)x = y are to be solved,
! 40: * but diagonal elements of U are not to be perturbed.
! 41: * = -1: The equations (T - lambda*I)x = y are to be solved
! 42: * and, if overflow would otherwise occur, the diagonal
! 43: * elements of U are to be perturbed. See argument TOL
! 44: * below.
! 45: * = 2: The equations (T - lambda*I)'x = y are to be solved,
! 46: * but diagonal elements of U are not to be perturbed.
! 47: * = -2: The equations (T - lambda*I)'x = y are to be solved
! 48: * and, if overflow would otherwise occur, the diagonal
! 49: * elements of U are to be perturbed. See argument TOL
! 50: * below.
! 51: *
! 52: * N (input) INTEGER
! 53: * The order of the matrix T.
! 54: *
! 55: * A (input) DOUBLE PRECISION array, dimension (N)
! 56: * On entry, A must contain the diagonal elements of U as
! 57: * returned from DLAGTF.
! 58: *
! 59: * B (input) DOUBLE PRECISION array, dimension (N-1)
! 60: * On entry, B must contain the first super-diagonal elements of
! 61: * U as returned from DLAGTF.
! 62: *
! 63: * C (input) DOUBLE PRECISION array, dimension (N-1)
! 64: * On entry, C must contain the sub-diagonal elements of L as
! 65: * returned from DLAGTF.
! 66: *
! 67: * D (input) DOUBLE PRECISION array, dimension (N-2)
! 68: * On entry, D must contain the second super-diagonal elements
! 69: * of U as returned from DLAGTF.
! 70: *
! 71: * IN (input) INTEGER array, dimension (N)
! 72: * On entry, IN must contain details of the matrix P as returned
! 73: * from DLAGTF.
! 74: *
! 75: * Y (input/output) DOUBLE PRECISION array, dimension (N)
! 76: * On entry, the right hand side vector y.
! 77: * On exit, Y is overwritten by the solution vector x.
! 78: *
! 79: * TOL (input/output) DOUBLE PRECISION
! 80: * On entry, with JOB .lt. 0, TOL should be the minimum
! 81: * perturbation to be made to very small diagonal elements of U.
! 82: * TOL should normally be chosen as about eps*norm(U), where eps
! 83: * is the relative machine precision, but if TOL is supplied as
! 84: * non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
! 85: * If JOB .gt. 0 then TOL is not referenced.
! 86: *
! 87: * On exit, TOL is changed as described above, only if TOL is
! 88: * non-positive on entry. Otherwise TOL is unchanged.
! 89: *
! 90: * INFO (output) INTEGER
! 91: * = 0 : successful exit
! 92: * .lt. 0: if INFO = -i, the i-th argument had an illegal value
! 93: * .gt. 0: overflow would occur when computing the INFO(th)
! 94: * element of the solution vector x. This can only occur
! 95: * when JOB is supplied as positive and either means
! 96: * that a diagonal element of U is very small, or that
! 97: * the elements of the right-hand side vector y are very
! 98: * large.
! 99: *
! 100: * =====================================================================
! 101: *
! 102: * .. Parameters ..
! 103: DOUBLE PRECISION ONE, ZERO
! 104: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 105: * ..
! 106: * .. Local Scalars ..
! 107: INTEGER K
! 108: DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
! 109: * ..
! 110: * .. Intrinsic Functions ..
! 111: INTRINSIC ABS, MAX, SIGN
! 112: * ..
! 113: * .. External Functions ..
! 114: DOUBLE PRECISION DLAMCH
! 115: EXTERNAL DLAMCH
! 116: * ..
! 117: * .. External Subroutines ..
! 118: EXTERNAL XERBLA
! 119: * ..
! 120: * .. Executable Statements ..
! 121: *
! 122: INFO = 0
! 123: IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
! 124: INFO = -1
! 125: ELSE IF( N.LT.0 ) THEN
! 126: INFO = -2
! 127: END IF
! 128: IF( INFO.NE.0 ) THEN
! 129: CALL XERBLA( 'DLAGTS', -INFO )
! 130: RETURN
! 131: END IF
! 132: *
! 133: IF( N.EQ.0 )
! 134: $ RETURN
! 135: *
! 136: EPS = DLAMCH( 'Epsilon' )
! 137: SFMIN = DLAMCH( 'Safe minimum' )
! 138: BIGNUM = ONE / SFMIN
! 139: *
! 140: IF( JOB.LT.0 ) THEN
! 141: IF( TOL.LE.ZERO ) THEN
! 142: TOL = ABS( A( 1 ) )
! 143: IF( N.GT.1 )
! 144: $ TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
! 145: DO 10 K = 3, N
! 146: TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
! 147: $ ABS( D( K-2 ) ) )
! 148: 10 CONTINUE
! 149: TOL = TOL*EPS
! 150: IF( TOL.EQ.ZERO )
! 151: $ TOL = EPS
! 152: END IF
! 153: END IF
! 154: *
! 155: IF( ABS( JOB ).EQ.1 ) THEN
! 156: DO 20 K = 2, N
! 157: IF( IN( K-1 ).EQ.0 ) THEN
! 158: Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
! 159: ELSE
! 160: TEMP = Y( K-1 )
! 161: Y( K-1 ) = Y( K )
! 162: Y( K ) = TEMP - C( K-1 )*Y( K )
! 163: END IF
! 164: 20 CONTINUE
! 165: IF( JOB.EQ.1 ) THEN
! 166: DO 30 K = N, 1, -1
! 167: IF( K.LE.N-2 ) THEN
! 168: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
! 169: ELSE IF( K.EQ.N-1 ) THEN
! 170: TEMP = Y( K ) - B( K )*Y( K+1 )
! 171: ELSE
! 172: TEMP = Y( K )
! 173: END IF
! 174: AK = A( K )
! 175: ABSAK = ABS( AK )
! 176: IF( ABSAK.LT.ONE ) THEN
! 177: IF( ABSAK.LT.SFMIN ) THEN
! 178: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
! 179: $ THEN
! 180: INFO = K
! 181: RETURN
! 182: ELSE
! 183: TEMP = TEMP*BIGNUM
! 184: AK = AK*BIGNUM
! 185: END IF
! 186: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
! 187: INFO = K
! 188: RETURN
! 189: END IF
! 190: END IF
! 191: Y( K ) = TEMP / AK
! 192: 30 CONTINUE
! 193: ELSE
! 194: DO 50 K = N, 1, -1
! 195: IF( K.LE.N-2 ) THEN
! 196: TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
! 197: ELSE IF( K.EQ.N-1 ) THEN
! 198: TEMP = Y( K ) - B( K )*Y( K+1 )
! 199: ELSE
! 200: TEMP = Y( K )
! 201: END IF
! 202: AK = A( K )
! 203: PERT = SIGN( TOL, AK )
! 204: 40 CONTINUE
! 205: ABSAK = ABS( AK )
! 206: IF( ABSAK.LT.ONE ) THEN
! 207: IF( ABSAK.LT.SFMIN ) THEN
! 208: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
! 209: $ THEN
! 210: AK = AK + PERT
! 211: PERT = 2*PERT
! 212: GO TO 40
! 213: ELSE
! 214: TEMP = TEMP*BIGNUM
! 215: AK = AK*BIGNUM
! 216: END IF
! 217: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
! 218: AK = AK + PERT
! 219: PERT = 2*PERT
! 220: GO TO 40
! 221: END IF
! 222: END IF
! 223: Y( K ) = TEMP / AK
! 224: 50 CONTINUE
! 225: END IF
! 226: ELSE
! 227: *
! 228: * Come to here if JOB = 2 or -2
! 229: *
! 230: IF( JOB.EQ.2 ) THEN
! 231: DO 60 K = 1, N
! 232: IF( K.GE.3 ) THEN
! 233: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
! 234: ELSE IF( K.EQ.2 ) THEN
! 235: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
! 236: ELSE
! 237: TEMP = Y( K )
! 238: END IF
! 239: AK = A( K )
! 240: ABSAK = ABS( AK )
! 241: IF( ABSAK.LT.ONE ) THEN
! 242: IF( ABSAK.LT.SFMIN ) THEN
! 243: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
! 244: $ THEN
! 245: INFO = K
! 246: RETURN
! 247: ELSE
! 248: TEMP = TEMP*BIGNUM
! 249: AK = AK*BIGNUM
! 250: END IF
! 251: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
! 252: INFO = K
! 253: RETURN
! 254: END IF
! 255: END IF
! 256: Y( K ) = TEMP / AK
! 257: 60 CONTINUE
! 258: ELSE
! 259: DO 80 K = 1, N
! 260: IF( K.GE.3 ) THEN
! 261: TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
! 262: ELSE IF( K.EQ.2 ) THEN
! 263: TEMP = Y( K ) - B( K-1 )*Y( K-1 )
! 264: ELSE
! 265: TEMP = Y( K )
! 266: END IF
! 267: AK = A( K )
! 268: PERT = SIGN( TOL, AK )
! 269: 70 CONTINUE
! 270: ABSAK = ABS( AK )
! 271: IF( ABSAK.LT.ONE ) THEN
! 272: IF( ABSAK.LT.SFMIN ) THEN
! 273: IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
! 274: $ THEN
! 275: AK = AK + PERT
! 276: PERT = 2*PERT
! 277: GO TO 70
! 278: ELSE
! 279: TEMP = TEMP*BIGNUM
! 280: AK = AK*BIGNUM
! 281: END IF
! 282: ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
! 283: AK = AK + PERT
! 284: PERT = 2*PERT
! 285: GO TO 70
! 286: END IF
! 287: END IF
! 288: Y( K ) = TEMP / AK
! 289: 80 CONTINUE
! 290: END IF
! 291: *
! 292: DO 90 K = N, 2, -1
! 293: IF( IN( K-1 ).EQ.0 ) THEN
! 294: Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
! 295: ELSE
! 296: TEMP = Y( K-1 )
! 297: Y( K-1 ) = Y( K )
! 298: Y( K ) = TEMP - C( K-1 )*Y( K )
! 299: END IF
! 300: 90 CONTINUE
! 301: END IF
! 302: *
! 303: * End of DLAGTS
! 304: *
! 305: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlagts.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack