LIN_SOL_TRI
Solves multiple systems of linear equations
Ajxj = yj, j = 1, …, k
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 for LIN_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_set_jolt, jolt)
Default: epsilon(), machine precision
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: CALL LIN_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 = Aj yj 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.
use lin_sol_tri_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 1 for LIN_SOL_TRI.
integer i
integer, parameter :: n=128
real(kind(1d0)), parameter :: one=1d0, zero=0d0
real(kind(1d0)) err
real(kind(1d0)), dimension(2*n,n) :: d, b, c, res(n,n), &
t(n), x, y
! Generate the upper, main, and lower diagonals of the
! n matrices A_i. For each system a random vector x is used
! to construct the right-hand side, Ax = y. The lower part
! of each array remains zero as a result.
c = zero; d=zero; b=zero; x=zero
do i = 1, n
call rand_gen (c(1:n,i))
call rand_gen (d(1:n,i))
call rand_gen (b(1:n,i))
call rand_gen (x(1:n,i))
end do
! Add scalars to the main diagonal of each system so that
! all systems are positive definite.
t = sum(c+d+b,DIM=1)
d(1:n,1:n) = d(1:n,1:n) + spread(t,DIM=1,NCOPIES=n)
! 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
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.
use lin_sol_tri_int
use rand_gen_int
implicit none
! This is Example 2 for LIN_SOL_TRI.
integer i, nopt
integer, parameter :: n=128
real(kind(1e0)), parameter :: s_one=1e0, s_zero=0e0
real(kind(1d0)), parameter :: d_one=1d0, d_zero=0d0
real(kind(1e0)), dimension(2*n,n) :: d, b, c, res(n,n), &
x, y
real(kind(1e0)) change_new, change_old, err
type(s_options) :: iopt(2) = s_options(0,s_zero)
real(kind(1d0)), dimension(n,n) :: d_save, b_save, c_save, &
x_save, y_save, x_sol
logical solve_only
c = s_zero; d=s_zero; b=s_zero; x=s_zero
! Generate the upper, main, and lower diagonals of the
! matrices A. A random vector x is used to construct the
! right-hand sides: y=A*x.
do i = 1, n
call rand_gen (c(1:n,i))
call rand_gen (d(1:n,i))
call rand_gen (b(1:n,i))
call rand_gen (x(1:n,i))
end do
! Save double precision copies of the diagonals and the
! right-hand side.
c_save = c(1:n,1:n); d_save = d(1:n,1:n)
b_save = b(1:n,1:n); x_save = x(1:n,1:n)
y_save(1:n,1:n) = d(1:n,1:n)*x_save + &
c(1:n,1:n)*EOSHIFT(x_save,SHIFT=+1,DIM=1) + &
b(1:n,1:n)*EOSHIFT(x_save,SHIFT=-1,DIM=1)
! Iterative refinement loop.
factorization_choice: do nopt=0, 1
! Set the logical to flag the first time through.
solve_only = .false.
x_sol = d_zero
change_old = huge(s_one)
iterative_refinement: do
! This flag causes a copy of data to be moved to work arrays
! and a factorization and solve step to be performed.
if (.not. solve_only) then
c(1:n,1:n)=c_save; d(1:n,1:n)=d_save
b(1:n,1:n)=b_save
end if
! Compute current residuals, y - A*x, using current x.
y(1:n,1:n) = -y_save + &
d_save*x_sol + &
c_save*EOSHIFT(x_sol,SHIFT=+1,DIM=1) + &
b_save*EOSHIFT(x_sol,SHIFT=-1,DIM=1)
call lin_sol_tri (c, d, b, y, iopt=iopt)
x_sol = x_sol + y(1:n,1:n)
change_new = sum(abs(y(1:n,1:n)))
! If size of change is not decreasing, stop the iteration.
if (change_new >= change_old) exit iterative_refinement
change_old = change_new
iopt(nopt+1) = s_options(s_lin_sol_tri_solve_only,s_zero)
solve_only = .true.
end do iterative_refinement
! Use Gaussian Elimination if Cyclic Reduction did not get an
! accurate solution.
! It is an exceptional event when Gaussian Elimination is required.
if (sum(abs(x_sol - x_save)) / sum(abs(x_save)) &
<= sqrt(epsilon(d_one))) exit factorization_choice
iopt = s_options(0,s_zero)
iopt(nopt+1) = s_options(s_lin_sol_tri_use_Gauss_elim,s_zero)
end do factorization_choice
! Check on accuracy of solution.
res = x(1:n,1:n)- x_save
err = sum(abs(res)) / sum(abs(x_save))
if (err <= sqrt(epsilon(d_one))) then
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.
use lin_sol_tri_int
use rand_gen_int
use Numerical_Libraries
implicit none
! This is Example 3 for LIN_SOL_TRI.
integer i, j, nopt
integer, parameter :: n=128, k=n/4, ncoda=1, lda=2
real(kind(1e0)), parameter :: s_one=1e0, s_zero=0e0
real(kind(1e0)) A(lda,n), EVAL(k)
type(s_options) :: iopt(2)=s_options(0,s_zero)
real(kind(1e0)) d(n), b(n), d_t(2*n,k), c_t(2*n,k), perf_ratio, &
b_t(2*n,k), y_t(2*n,k), eval_t(k), res(n,k), temp
logical small
! This flag is used to get the k largest eigenvalues.
small = .false.
! Generate the main diagonal and the co-diagonal of the
! tridiagonal matrix.
call rand_gen (b)
call rand_gen (d)
A(1,1:)=b; A(2,1:)=d
! Use Numerical Libraries routine for the calculation of k
! largest eigenvalues.
CALL EVASB (N, K, A, LDA, NCODA, SMALL, EVAL)
EVAL_T = EVAL
! Use DNFL tridiagonal solver for inverse iteration
! calculation of eigenvectors.
factorization_choice: do nopt=0,1
! Create k tridiagonal problems, one for each inverse
! iteration system.
b_t(1:n,1:k) = spread(b,DIM=2,NCOPIES=k)
c_t(1:n,1:k) = EOSHIFT(b_t(1:n,1:k),SHIFT=1,DIM=1)
d_t(1:n,1:k) = spread(d,DIM=2,NCOPIES=k) - &
spread(EVAL_T,DIM=1,NCOPIES=n)
! Start the right-hand side at random values, scaled downward
! to account for the expected 'blowup' in the solution.
do i=1, k
call rand_gen (y_t(1:n,i))
end do
! Do two iterations for the eigenvectors.
do i=1, 2
y_t(1:n,1:k) = y_t(1:n,1:k)*epsilon(s_one)
call lin_sol_tri(c_t, d_t, b_t, y_t, &
iopt=iopt)
iopt(nopt+1) = s_options(s_lin_sol_tri_solve_only,s_zero)
end do
! Orthogonalize the eigenvectors. (This is the most
! intensive part of the computing.)
do j=1,k-1 ! Forward sweep of HMGS orthogonalization.
temp=s_one/sqrt(sum(y_t(1:n,j)**2))
y_t(1:n,j)=y_t(1:n,j)*temp
y_t(1:n,j+1:k)=y_t(1:n,j+1:k)+ &
spread(-matmul(y_t(1:n,j),y_t(1:n,j+1:k)), &
DIM=1,NCOPIES=n)* spread(y_t(1:n,j),DIM=2,NCOPIES=k-j)
end do
temp=s_one/sqrt(sum(y_t(1:n,k)**2))
y_t(1:n,k)=y_t(1:n,k)*temp
do j=k-1,1,-1 ! Backward sweep of HMGS.
y_t(1:n,j+1:k)=y_t(1:n,j+1:k)+ &
spread(-matmul(y_t(1:n,j),y_t(1:n,j+1:k)), &
DIM=1,NCOPIES=n)* spread(y_t(1:n,j),DIM=2,NCOPIES=k-j)
end do
! See if the performance ratio is smaller than the value one.
! If it is not the code will re-solve the systems using Gaussian
! Elimination. This is an exceptional event. It is a necessary
! complication for achieving reliable results.
res(1:n,1:k) = spread(d,DIM=2,NCOPIES=k)*y_t(1:n,1:k) + &
spread(b,DIM=2,NCOPIES=k)* &
EOSHIFT(y_t(1:n,1:k),SHIFT=-1,DIM=1) + &
EOSHIFT(spread(b,DIM=2,NCOPIES=k)*y_t(1:n,1:k),SHIFT=1) &
-y_t(1:n,1:k)*spread(EVAL_T(1:k),DIM=1,NCOPIES=n)
! If the factorization method is Cyclic Reduction and perf_ratio is
! larger than one, re-solve using Gaussian Elimination. If the
! method is already Gaussian Elimination, the loop exits
! and perf_ratio is checked at the end.
perf_ratio = sum(abs(res(1:n,1:k))) / &
sum(abs(EVAL_T(1:k))) / &
epsilon(s_one) / (5*n)
if (perf_ratio <= s_one) exit factorization_choice
iopt(nopt+1) = s_options(s_lin_sol_tri_use_Gauss_elim,s_zero)
end do factorization_choice
if (perf_ratio <= s_one) then
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.
use lin_sol_tri_int
use rand_gen_int
use Numerical_Libraries
implicit none
! This is Example 4 for LIN_SOL_TRI.
integer, parameter :: n=1000, ichap=5, iget=1, iput=2, &
inum=6, irnum=7
real(kind(1e0)), parameter :: zero=0e0, one = 1e0
integer i, ido, in(50), inr(20), iopt(6), ival(7), &
iwk(35+n)
real(kind(1e0)) hx, pi_value, t, u_0, u_1, atol, rtol, sval(2), &
tend, wk(41+11*n), y(n), ypr(n), a_diag(n), &
a_off(n), r_diag(n), r_off(n), t_y(n), t_ypr(n), &
t_g(n), t_diag(2*n,1), t_upper(2*n,1), &
t_lower(2*n,1), t_sol(2*n,1)
type(s_options) :: iopti(2)=s_options(0,zero)
character(2) :: pi(1) = 'pi'
! Define initial data.
t = 0.0e0
u_0 = 1
u_1 = 0.5
tend = one
! Initial values for the variational equation.
y = -one; ypr= zero
pi_value = const(pi)
hx = pi_value/(n+1)
a_diag = 2*hx/3
a_off = hx/6
r_diag = -2/hx
r_off = 1/hx
! Get integer option numbers.
iopt(1) = inum
call iumag ('math', ichap, iget, 1, iopt, in)
! Get floating point option numbers.
iopt(1) = irnum
call iumag ('math', ichap, iget, 1, iopt, inr)
! Set for reverse communication evaluation of the DAE.
iopt(1) = in(26)
ival(1) = 0
! Set for use of explicit partial derivatives.
iopt(2) = in(5)
ival(2) = 1
! Set for reverse communication evaluation of partials.
iopt(3) = in(29)
ival(3) = 0
! Set for reverse communication solution of linear equations.
iopt(4) = in(31)
ival(4) = 0
! Storage for the partial derivative array are not allocated or
! required in the integrator.
iopt(5) = in(34)
ival(5) = 1
! Set the sizes of iwk, wk for internal checking.
iopt(6) = in(35)
ival(6) = 35 + n
ival(7) = 41 + 11*n
! Set integer options:
call iumag ('math', ichap, iput, 6, iopt, ival)
! Reset tolerances for integrator:
atol = 1e-3; rtol= 1e-3
sval(1) = atol; sval(2) = rtol
iopt(1) = inr(5)
! Set floating point options:
call sumag ('math', ichap, iput, 1, iopt, sval)
! Integrate ODE/DAE. Use dummy external names for g(y,y')
! and partials.
ido = 1
Integration_Loop: do
call d2spg (n, t, tend, ido, y, ypr, dgspg, djspg, iwk, wk)
! Find where g(y,y') goes. (It only goes in one place here, but can
! vary where divided differences are used for partial derivatives.)
iopt(1) = in(27)
call iumag ('math', ichap, iget, 1, iopt, ival)
! Direct user response:
select case(ido)
case(1,4)
! This should not occur.
write (*,*) ' Unexpected return with ido = ', ido
stop
case(3)
! Reset options to defaults. (This is good housekeeping but not
! required for this problem.)
in = -in
call iumag ('math', ichap, iput, 50, in, ival)
inr = -inr
call sumag ('math', ichap, iput, 20, inr, sval)
exit Integration_Loop
case(5)
! Evaluate partials of g(y,y').
t_y = y; t_ypr = ypr
t_g = r_diag*t_y + r_off*EOSHIFT(t_y,SHIFT=+1) &
+ EOSHIFT(r_off*t_y,SHIFT=-1) &
- (a_diag*t_ypr + a_off*EOSHIFT(t_ypr,SHIFT=+1) &
+ EOSHIFT(a_off*t_ypr,SHIFT=-1))
! Move data from the assumed size to assumed shape arrays.
do i=1, n
wk(ival(1)+i-1) = t_g(i)
end do
cycle Integration_Loop
case(6)
! Evaluate partials of g(y,y').
! Get value of c_j for partials.
iopt(1) = inr(9)
call sumag ('math', ichap, iget, 1, iopt, sval)
! Subtract c_j from diagonals to compute (partials for y')*c_j.
! The linear system is tridiagonal.
t_diag(1:n,1) = r_diag - sval(1)*a_diag
t_upper(1:n,1) = r_off - sval(1)*a_off
t_lower = EOSHIFT(t_upper,SHIFT=+1,DIM=1)
cycle Integration_Loop
case(7)
! Compute the factorization.
iopti(1) = s_options(s_lin_sol_tri_factor_only,zero)
call lin_sol_tri (t_upper, t_diag, t_lower, &
t_sol, iopt=iopti)
cycle Integration_Loop
case(8)
! Solve the system.
iopti(1) = s_options(s_lin_sol_tri_solve_only,zero)
! Move data from the assumed size to assumed shape arrays.
t_sol(1:n,1)=wk(ival(1):ival(1)+n-1)
call lin_sol_tri (t_upper, t_diag, t_lower, &
t_sol, iopt=iopti)
! Move data from the assumed shape to assumed size arrays.
wk(ival(1):ival(1)+n-1)=t_sol(1:n,1)
cycle Integration_Loop
case(2)
! Correct initial value to reach u_1 at t=tend.
u_0 = u_0 - (u_0*y(n/2) - (u_1-u_0)) / (y(n/2) + 1)
! Finish up internally in the integrator.
ido = 3
cycle Integration_Loop
end select
end do Integration_Loop
write (*,*) 'The equation u_t = u_xx, with u(0,t) = ', u_0
write (*,*) 'reaches the value ',u_1, ' at time = ', tend, '.'
write (*,*) 'Example 4 for LIN_SOL_TRI is correct.'
end
Output
Example 4 for LIN_SOL_TRI is correct.
