Each matrix Aj is tridiagonal with the same dimension, n. The default solution method is based on LU factorization computed using cyclic reduction or, optionally, Gaussian elimination with partial pivoting.
Required Arguments
C — Array of size 2n×k containing the upper diagonals of the matrices Aj. Each upper diagonal is entered in array locations c(1:n – 1, j). The data C(n, 1:k) are not used. (Input [/Output]) The input data is overwritten. See note below.
D — Array of size 2n×k containing the diagonals of the matrices Aj. Each diagonal is entered in array locations D(1:n, j). (Input [/Output]) The input data is overwritten. See note below.
B — Array of size 2n×k containing the lower diagonals of the matrices Aj. Each lower diagonal is entered in array locations B(2:n, j). The data B(1, 1:k) are not used. (Input [/Output]) The input data is overwritten. See note below.
Y — Array of size 2n×k containing the right-hand sides, yj. Each right-hand side is entered in array locations Y(1:n, j). The computed solution xj is returned in locations Y(1:n, j). (Input [/Output])
NOTE: The required arguments have the Input data overwritten. If these quantities are used later, they must be saved in user-defined arrays. The routine uses each array's locations (n + 1:2 * n, 1:k) for scratch storage and intermediate data in the LU factorization. The default values for problem dimensions are n = (size (D, 1))/2 and k = size (D, 2).
Optional Arguments
NCOLS= n (Input) Uses arrays C(1:n – 1, 1:k), D(1:n, 1:k), and B(2:n, 1:k) as the upper, main and lower diagonals for the input tridiagonal matrices. The right-hand sides and solutions are in array Y(1:n, 1:k). Note that each of these arrays are rank-2. Default: n = (size(D, 1))/2
NPROB= k (Input) The number of systems solved. Default: k = size(D, 2)
iopt= iopt(:) (Input) Derived type array with the same precision as the input matrix. Used for passing optional data to the routine. The options are as follows:
Packaged Options forLIN_SOL_TRI
Option Prefix = ?
Option Name
Option Value
s_, d_, c_, z_
lin_sol_tri_set_small
1
s_, d_, c_, z_
lin_sol_tri_set_jolt
2
s_, d_, c_, z_
lin_sol_tri_scan_for_NaN
3
s_, d_, c_, z_
lin_sol_tri_factor_only
4
s_, d_, c_, z_
lin_sol_tri_solve_only
5
s_, d_, c_, z_
lin_sol_tri_use_Gauss_elim
6
iopt(IO)= ?_options(?_lin_sol_tri_set_small, Small) Whenever a reciprocation is performed on a quantity smaller than Small, it is replaced by that value plus 2 ×jolt. Default: 0.25 ×epsilon()
iopt(IO)= ?_options(?_lin_sol_tri_scan_for_NaN, ?_dummy) Examines each input array entry to find the first value such that
isNaN(C(i,j)) .or.
isNaN(D(i,j)) .or.
isNaN(B(i,j)) .or.
isNaN(Y(i,j)) == .true.
See the isNaN() function, Chapter 10. Default: Does not scan for NaNs.
iopt(IO)= ?_options(?_lin_sol_tri_factor_only, ?_dummy) Obtain the LU factorization of the matrices Aj. Does not solve for a solution. Default: Factor the matrices and solve the systems.
iopt(IO)= ?_options(?_lin_sol_tri_solve_only, ?_dummy) Solve the systems Ajxj = yj using the previously computed LU factorization. Default: Factor the matrices and solve the systems.
iopt(IO)= ?_options(?_lin_sol_tri_use_Gauss_elim, ?_dummy) The accuracy, numerical stability or efficiency of the cyclic reduction algorithm may be inferior to the use of LU factorization with partial pivoting. Default: Use cyclic reduction to compute the factorization.
FORTRAN 90 Interface
Generic: CALLLIN_SOL_TRI (C, D, B, Y[, …])
Specific: The specific interface names are S_LIN_SOL_TRI, D_LIN_SOL_TRI, C_LIN_SOL_TRI, and Z_LIN_SOL_TRI.
Description
Routine lin_sol_tri solves k systems of tridiagonal linear algebraic equations, each problem of dimension n×n. No relation between k and n is required. See Kershaw, pages 86–88 in Rodrigue (1982) for further details. To deal with poorly conditioned or singular systems, a specific regularizing term is added to each reciprocated value. This technique keeps the factorization process efficient and avoids exceptions from overflow or division by zero. Each occurrence of an array reciprocal a-1 is replaced by the expression (a + t)-1, where the array temporary t has the value 0 whenever the corresponding entry satisfies |a| > Small. Alternately, t has the value 2 ×jolt. (Every small denominator gives rise to a finite “jolt”.) Since this tridiagonal solver is used in the routines lin_svd and lin_eig_self for inverse iteration, regularization is required. Users can reset the values of Small and jolt for their own needs. Using the default values for these parameters, it is generally necessary to scale the tridiagonal matrix so that the maximum magnitude has value approximately one. This is normally not an issue when the systems are nonsingular.
The routine is designed to use cyclic reduction as the default method for computing the LU factorization. Using an optional parameter, standard elimination and partial pivoting will be used to compute the factorization. Partial pivoting is numerically stable but is likely to be less efficient than cyclic reduction.
Fatal, Terminal, and Warning Error Messages
See the messages.gls file for error messages for LIN_SOL_TRI. These error messages are numbered 1081–1086; 1101–1106; 1121–1126; 1141–1146.
Examples
Example 1: Solution of Multiple Tridiagonal Systems
The upper, main and lower diagonals of n systems of size n×n are generated randomly. A scalar is added to the main diagonal so that the systems are positive definite. A random vector xj is generated and right-hand sides yj = Ajyj are computed. The routine is used to compute the solution, using the Aj and yj. The results should compare closely with the xj used to generate the right-hand sides. Also, see operator_ex17, supplied with the product examples.
! Set Ax = y. The vector x generates y. Note the use
! of EOSHIFT and array operations to compute the matrix
! product, n distinct ones as one array operation.
y(1:n,1:n)=d(1:n,1:n)*x(1:n,1:n) + &
c(1:n,1:n)*EOSHIFT(x(1:n,1:n),SHIFT=+1,DIM=1) + &
b(1:n,1:n)*EOSHIFT(x(1:n,1:n),SHIFT=-1,DIM=1)
! Compute the solution returned in y. (The input values of c,
! d, b, and y are overwritten by lin_sol_tri.) Check for any
! error messages.
call lin_sol_tri (c, d, b, y)
! Check the size of the residuals, y-x. They should be small,
! relative to the size of values in x.
res = x(1:n,1:n) - y(1:n,1:n)
err = sum(abs(res)) / sum(abs(x(1:n,1:n)))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 1 for LIN_SOL_TRI is correct.'
end if
end
Output
Example 1 for LIN_SOL_TRI is correct.
Example 2: Iterative Refinement and Use of Partial Pivoting
This program unit shows usage that typically gives acceptable accuracy for a large class of problems. Our goal is to use the efficient cyclic reduction algorithm when possible, and keep on using it unless it will fail. In exceptional cases our program switches to the LU factorization with partial pivoting. This use of both factorization and solution methods enhances reliability and maintains efficiency on the average. Also, see operator_ex18, supplied with the product examples.
write (*,*) 'Example 2 for LIN_SOL_TRI is correct.'
end if
end
Output
Example 2 for LIN_SOL_TRI is correct.
Example 3: Eigenvectors of Tridiagonal Matrices
The eigenvalues λ1, .…λn of a tridiagonal real, self-adjoint matrix are computed. Note that the computation is performed using the IMSL MATH/LIBRARY FORTRAN 77 interface to routine EVASB. The user may write this interface based on documentation of the arguments (IMSL 2003, p. 480), or use the module Numerical_Libraries as we have done here. The eigenvectors corresponding to k < n of the eigenvalues are required. These vectors are computed using inverse iteration for all the eigenvalues at one step. See Golub and Van Loan (1989, Chapter 7). The eigenvectors are then orthogonalized. Also, see operator_ex19, supplied with the product examples.
write (*,*) 'Example 3 for LIN_SOL_TRI is correct.'
end if
end
Output
Example 3 for LIN_SOL_TRI is correct.
Example 4: Tridiagonal Matrix Solving within Diffusion Equations
The normalized partial differential equation
is solved for values of 0 ≤x≤π and t > 0. A boundary value problem consists of choosing the value
such that the equation
is satisfied. Arbitrary values
and
are used for illustration of the solution process. The one-parameter equation
The variables are changed to
that v(0, t) = 0. The function v(x, t) satisfies the differential equation. The one-parameter equation solved is therefore
To solve this equation for , use the standard technique of the variational equation,
Thus
Since the initial data for
the variational equation initial condition is
w(x, 0) = –1
This model problem illustrates the method of lines and Galerkin principle implemented with the differential-algebraic solver, D2SPG (IMSL 2003, pp. 889–911). We use the integrator in “reverse communication” mode for evaluating the required functions, derivatives, and solving linear algebraic equations. See Example 4 of routine DASPG for a problem that uses reverse communication. Next see Example 4 of routine IVPAG for the development of the piecewise-linear Galerkin discretization method to solve the differential equation. This present example extends parts of both previous examples and illustrates Fortran 90 constructs. It further illustrates how a user can deal with a defect of an integrator that normally functions using only dense linear algebra factorization methods for solving the corrector equations. See the comments in Brenan et al. (1989, esp. p. 137). Also, see operator_ex20, supplied with the product examples.
