LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlagts()

subroutine dlagts ( integer job,
integer n,
double precision, dimension( * ) a,
double precision, dimension( * ) b,
double precision, dimension( * ) c,
double precision, dimension( * ) d,
integer, dimension( * ) in,
double precision, dimension( * ) y,
double precision tol,
integer info )

DLAGTS solves the system of equations (T-λI)x = y or (T-λI)^Tx = y, where T is a general tridiagonal matrix and λ a scalar, using the LU factorization computed by slagtf.

Download DLAGTS + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLAGTS may be used to solve one of the systems of equations
!>
!>    (T - lambda*I)*x = y   or   (T - lambda*I)**T*x = y,
!>
!> where T is an n by n tridiagonal matrix, for x, following the
!> factorization of (T - lambda*I) as
!>
!>    (T - lambda*I) = P*L*U ,
!>
!> by routine DLAGTF. The choice of equation to be solved is
!> controlled by the argument JOB, and in each case there is an option
!> to perturb zero or very small diagonal elements of U, this option
!> being intended for use in applications such as inverse iteration.
!>
Parameters
[in]JOB
!>          JOB is INTEGER
!>          Specifies the job to be performed by DLAGTS as follows:
!>          =  1: The equations  (T - lambda*I)x = y  are to be solved,
!>                but diagonal elements of U are not to be perturbed.
!>          = -1: The equations  (T - lambda*I)x = y  are to be solved
!>                and, if overflow would otherwise occur, the diagonal
!>                elements of U are to be perturbed. See argument TOL
!>                below.
!>          =  2: The equations  (T - lambda*I)**Tx = y  are to be solved,
!>                but diagonal elements of U are not to be perturbed.
!>          = -2: The equations  (T - lambda*I)**Tx = y  are to be solved
!>                and, if overflow would otherwise occur, the diagonal
!>                elements of U are to be perturbed. See argument TOL
!>                below.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix T.
!>
[in]A
!>          A is DOUBLE PRECISION array, dimension (N)
!>          On entry, A must contain the diagonal elements of U as
!>          returned from DLAGTF.
!>
[in]B
!>          B is DOUBLE PRECISION array, dimension (N-1)
!>          On entry, B must contain the first super-diagonal elements of
!>          U as returned from DLAGTF.
!>
[in]C
!>          C is DOUBLE PRECISION array, dimension (N-1)
!>          On entry, C must contain the sub-diagonal elements of L as
!>          returned from DLAGTF.
!>
[in]D
!>          D is DOUBLE PRECISION array, dimension (N-2)
!>          On entry, D must contain the second super-diagonal elements
!>          of U as returned from DLAGTF.
!>
[in]IN
!>          IN is INTEGER array, dimension (N)
!>          On entry, IN must contain details of the matrix P as returned
!>          from DLAGTF.
!>
[in,out]Y
!>          Y is DOUBLE PRECISION array, dimension (N)
!>          On entry, the right hand side vector y.
!>          On exit, Y is overwritten by the solution vector x.
!>
[in,out]TOL
!>          TOL is DOUBLE PRECISION
!>          On entry, with  JOB < 0, TOL should be the minimum
!>          perturbation to be made to very small diagonal elements of U.
!>          TOL should normally be chosen as about eps*norm(U), where eps
!>          is the relative machine precision, but if TOL is supplied as
!>          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
!>          If  JOB > 0  then TOL is not referenced.
!>
!>          On exit, TOL is changed as described above, only if TOL is
!>          non-positive on entry. Otherwise TOL is unchanged.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  overflow would occur when computing the INFO(th)
!>                element of the solution vector x. This can only occur
!>                when JOB is supplied as positive and either means
!>                that a diagonal element of U is very small, or that
!>                the elements of the right-hand side vector y are very
!>                large.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 160 of file dlagts.f.

161*
162* -- LAPACK auxiliary routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166* .. Scalar Arguments ..
167 INTEGER INFO, JOB, N
168 DOUBLE PRECISION TOL
169* ..
170* .. Array Arguments ..
171 INTEGER IN( * )
172 DOUBLE PRECISION A( * ), B( * ), C( * ), D( * ), Y( * )
173* ..
174*
175* =====================================================================
176*
177* .. Parameters ..
178 DOUBLE PRECISION ONE, ZERO
179 parameter( one = 1.0d+0, zero = 0.0d+0 )
180* ..
181* .. Local Scalars ..
182 INTEGER K
183 DOUBLE PRECISION ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
184* ..
185* .. Intrinsic Functions ..
186 INTRINSIC abs, max, sign
187* ..
188* .. External Functions ..
189 DOUBLE PRECISION DLAMCH
190 EXTERNAL dlamch
191* ..
192* .. External Subroutines ..
193 EXTERNAL xerbla
194* ..
195* .. Executable Statements ..
196*
197 info = 0
198 IF( ( abs( job ).GT.2 ) .OR. ( job.EQ.0 ) ) THEN
199 info = -1
200 ELSE IF( n.LT.0 ) THEN
201 info = -2
202 END IF
203 IF( info.NE.0 ) THEN
204 CALL xerbla( 'DLAGTS', -info )
205 RETURN
206 END IF
207*
208 IF( n.EQ.0 )
209 $ RETURN
210*
211 eps = dlamch( 'Epsilon' )
212 sfmin = dlamch( 'Safe minimum' )
213 bignum = one / sfmin
214*
215 IF( job.LT.0 ) THEN
216 IF( tol.LE.zero ) THEN
217 tol = abs( a( 1 ) )
218 IF( n.GT.1 )
219 $ tol = max( tol, abs( a( 2 ) ), abs( b( 1 ) ) )
220 DO 10 k = 3, n
221 tol = max( tol, abs( a( k ) ), abs( b( k-1 ) ),
222 $ abs( d( k-2 ) ) )
223 10 CONTINUE
224 tol = tol*eps
225 IF( tol.EQ.zero )
226 $ tol = eps
227 END IF
228 END IF
229*
230 IF( abs( job ).EQ.1 ) THEN
231 DO 20 k = 2, n
232 IF( in( k-1 ).EQ.0 ) THEN
233 y( k ) = y( k ) - c( k-1 )*y( k-1 )
234 ELSE
235 temp = y( k-1 )
236 y( k-1 ) = y( k )
237 y( k ) = temp - c( k-1 )*y( k )
238 END IF
239 20 CONTINUE
240 IF( job.EQ.1 ) THEN
241 DO 30 k = n, 1, -1
242 IF( k.LE.n-2 ) THEN
243 temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
244 ELSE IF( k.EQ.n-1 ) THEN
245 temp = y( k ) - b( k )*y( k+1 )
246 ELSE
247 temp = y( k )
248 END IF
249 ak = a( k )
250 absak = abs( ak )
251 IF( absak.LT.one ) THEN
252 IF( absak.LT.sfmin ) THEN
253 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
254 $ THEN
255 info = k
256 RETURN
257 ELSE
258 temp = temp*bignum
259 ak = ak*bignum
260 END IF
261 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
262 info = k
263 RETURN
264 END IF
265 END IF
266 y( k ) = temp / ak
267 30 CONTINUE
268 ELSE
269 DO 50 k = n, 1, -1
270 IF( k.LE.n-2 ) THEN
271 temp = y( k ) - b( k )*y( k+1 ) - d( k )*y( k+2 )
272 ELSE IF( k.EQ.n-1 ) THEN
273 temp = y( k ) - b( k )*y( k+1 )
274 ELSE
275 temp = y( k )
276 END IF
277 ak = a( k )
278 pert = sign( tol, ak )
279 40 CONTINUE
280 absak = abs( ak )
281 IF( absak.LT.one ) THEN
282 IF( absak.LT.sfmin ) THEN
283 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
284 $ THEN
285 ak = ak + pert
286 pert = 2*pert
287 GO TO 40
288 ELSE
289 temp = temp*bignum
290 ak = ak*bignum
291 END IF
292 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
293 ak = ak + pert
294 pert = 2*pert
295 GO TO 40
296 END IF
297 END IF
298 y( k ) = temp / ak
299 50 CONTINUE
300 END IF
301 ELSE
302*
303* Come to here if JOB = 2 or -2
304*
305 IF( job.EQ.2 ) THEN
306 DO 60 k = 1, n
307 IF( k.GE.3 ) THEN
308 temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
309 ELSE IF( k.EQ.2 ) THEN
310 temp = y( k ) - b( k-1 )*y( k-1 )
311 ELSE
312 temp = y( k )
313 END IF
314 ak = a( k )
315 absak = abs( ak )
316 IF( absak.LT.one ) THEN
317 IF( absak.LT.sfmin ) THEN
318 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
319 $ THEN
320 info = k
321 RETURN
322 ELSE
323 temp = temp*bignum
324 ak = ak*bignum
325 END IF
326 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
327 info = k
328 RETURN
329 END IF
330 END IF
331 y( k ) = temp / ak
332 60 CONTINUE
333 ELSE
334 DO 80 k = 1, n
335 IF( k.GE.3 ) THEN
336 temp = y( k ) - b( k-1 )*y( k-1 ) - d( k-2 )*y( k-2 )
337 ELSE IF( k.EQ.2 ) THEN
338 temp = y( k ) - b( k-1 )*y( k-1 )
339 ELSE
340 temp = y( k )
341 END IF
342 ak = a( k )
343 pert = sign( tol, ak )
344 70 CONTINUE
345 absak = abs( ak )
346 IF( absak.LT.one ) THEN
347 IF( absak.LT.sfmin ) THEN
348 IF( absak.EQ.zero .OR. abs( temp )*sfmin.GT.absak )
349 $ THEN
350 ak = ak + pert
351 pert = 2*pert
352 GO TO 70
353 ELSE
354 temp = temp*bignum
355 ak = ak*bignum
356 END IF
357 ELSE IF( abs( temp ).GT.absak*bignum ) THEN
358 ak = ak + pert
359 pert = 2*pert
360 GO TO 70
361 END IF
362 END IF
363 y( k ) = temp / ak
364 80 CONTINUE
365 END IF
366*
367 DO 90 k = n, 2, -1
368 IF( in( k-1 ).EQ.0 ) THEN
369 y( k-1 ) = y( k-1 ) - c( k-1 )*y( k )
370 ELSE
371 temp = y( k-1 )
372 y( k-1 ) = y( k )
373 y( k ) = temp - c( k-1 )*y( k )
374 END IF
375 90 CONTINUE
376 END IF
377*
378* End of DLAGTS
379*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
Here is the call graph for this function:
Here is the caller graph for this function: