LIN_SOL_SELF
Solves a system of linear equations Ax = b, where A is a self-adjoint matrix. Using optional arguments, any of several related computations can be performed. These extra tasks include computing and saving the factorization of A using symmetric pivoting, representing the determinant of A, computing the inverse matrix A-1, or computing the solution of Ax = b given the factorization of A. An optional argument is provided indicating that A is positive definite so that the Cholesky decomposition can be used.
Required Arguments
A — Array of size n × n containing the self-adjoint matrix. (Input [/Output])
If the packaged option lin_sol_self_save_factors is used then the factorization of A is saved in A. For solving efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A when the Cholesky method is used.
If the packaged option lin_sol_self_save_factors is used then the factorization of A is saved in A. For solving efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A when the Cholesky method is used.
B — Array of size n × nb containing the right-hand side matrix. (Input [/Output])
If the packaged option lin_sol_self_save_factors is used then input B is used as work storage and is not saved.
If the packaged option lin_sol_self_save_factors is used then input B is used as work storage and is not saved.
X — Array of size n × nb containing the solution matrix. (Output)
Optional Arguments
NROWS = n (Input)
Uses array A(1:n, 1:n) for the input matrix.
Default: n = size(A, 1)
Uses array A(1:n, 1:n) for the input matrix.
Default: n = size(A, 1)
NRHS = nb (Input)
Uses the array b(1:n, 1:nb) for the input right-hand side matrix.
Default: nb = size(b, 2)
Note that b must be a rank-2 array.
Uses the array b(1:n, 1:nb) for the input right-hand side matrix.
Default: nb = size(b, 2)
Note that b must be a rank-2 array.
pivots = pivots(:) (Output [/Input])
Integer array of size n + 1 that contains the individual row interchanges in the first n locations. Applied in order, these yield the permutation matrix P. Location n + 1 contains the number of the first diagonal term no larger than Small, which is defined on the next page of this chapter.
Integer array of size n + 1 that contains the individual row interchanges in the first n locations. Applied in order, these yield the permutation matrix P. Location n + 1 contains the number of the first diagonal term no larger than Small, which is defined on the next page of this chapter.
det = det(1:2) (Output)
Array of size 2 of the same type and kind as A for representing the determinant of the input matrix. The determinant is represented by two numbers. The first is the base with the sign or complex angle of the result. The second is the exponent. When det(2) is within exponent range, the value of the determinant is given by the expression abs(det(1))**det(2) * (det(1))/abs(det(1)). If the matrix is not singular, abs(det(1)) = radix(det); otherwise, det(1) = 0, and det(2) = ‑huge(abs(det(1))).
Array of size 2 of the same type and kind as A for representing the determinant of the input matrix. The determinant is represented by two numbers. The first is the base with the sign or complex angle of the result. The second is the exponent. When det(2) is within exponent range, the value of the determinant is given by the expression abs(det(1))**det(2) * (det(1))/abs(det(1)). If the matrix is not singular, abs(det(1)) = radix(det); otherwise, det(1) = 0, and det(2) = ‑huge(abs(det(1))).
ainv = ainv(:,:) (Output)
Array of the same type and kind as A(1:n, 1:n). It contains the inverse matrix, A-1 when the input matrix is nonsingular.
Array of the same type and kind as A(1:n, 1:n). It contains the inverse matrix, A-1 when the input matrix is nonsingular.
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:
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_self | ||
|---|---|---|
Option Prefix = ? | Option Name | Option Value |
s_, d_, c_, z_ | lin_sol_self_set_small | 1 |
s_, d_, c_, z_ | lin_sol_self_save_factors | 2 |
s_, d_, c_, z_ | lin_sol_self_no_pivoting | 3 |
s_, d_, c_, z_ | lin_sol_self_use_Cholesky | 4 |
s_, d_, c_, z_ | lin_sol_self_solve_A | 5 |
s_, d_, c_, z_ | lin_sol_self_scan_for_NaN | 6 |
s_, d_, c_, z_ | lin_sol_self_no_sing_mess | 7 |
iopt(IO) = ?_options(?_lin_sol_self_set_small, Small)
When Aasen’s method is used, the tridiagonal system Tu = v is solved using LU factorization with partial pivoting. If a diagonal term of the matrix U is smaller in magnitude than the value Small, it is replaced by Small. The system is declared singular. When the Cholesky method is used, the upper-triangular matrix R, (see Description), is obtained. If a diagonal term of the matrix R is smaller in magnitude than the value Small, it is replaced by Small. A solution is approximated based on this replacement in either case.
Default: the smallest number that can be reciprocated safely
When Aasen’s method is used, the tridiagonal system Tu = v is solved using LU factorization with partial pivoting. If a diagonal term of the matrix U is smaller in magnitude than the value Small, it is replaced by Small. The system is declared singular. When the Cholesky method is used, the upper-triangular matrix R, (see Description), is obtained. If a diagonal term of the matrix R is smaller in magnitude than the value Small, it is replaced by Small. A solution is approximated based on this replacement in either case.
Default: the smallest number that can be reciprocated safely
iopt(IO) = ?_options(?_lin_sol_self_save_factors, ?_dummy)
Saves the factorization of A. Requires the optional argument “pivots=” if the routine will be used for solving further systems with the same matrix. This is the only case where the input arrays A and b are not saved. For solving efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A when the Cholesky method is used.
Saves the factorization of A. Requires the optional argument “pivots=” if the routine will be used for solving further systems with the same matrix. This is the only case where the input arrays A and b are not saved. For solving efficiency, the diagonal reciprocals of the matrix R are saved in the diagonal entries of A when the Cholesky method is used.
iopt(IO) = ?_options(?_lin_sol_self_no_pivoting, ?_dummy)
Does no row pivoting. The array pivots(:), if present, satisfies pivots(i) = i + 1 for i = 1, …, n – 1 when using Aasen’s method. When using the Cholesky method, pivots(i) = i for i = 1, …, n.
Does no row pivoting. The array pivots(:), if present, satisfies pivots(i) = i + 1 for i = 1, …, n – 1 when using Aasen’s method. When using the Cholesky method, pivots(i) = i for i = 1, …, n.
iopt(IO) = ?_options(?_lin_sol_self_use_Cholesky, ?_dummy)
The Cholesky decomposition PAPT = RTR is used instead of the Aasen method.
The Cholesky decomposition PAPT = RTR is used instead of the Aasen method.
iopt(IO) = ?_options(?_lin_sol_self_solve_A, ?_dummy)
Uses the factorization of A computed and saved to solve Ax = b.
Uses the factorization of A computed and saved to solve Ax = b.
iopt(IO) = ?_options(?_lin_sol_self_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that
Examines each input array entry to find the first value such that
isNaN(a(i,j)) .or. isNan(b(i,j)) == .true.
iopt(IO) = ?_options(?_lin_sol_self_no_sing_mess, ?_dummy)
Do not print an error message when the matrix A is singular.
Do not print an error message when the matrix A is singular.
FORTRAN 90 Interface
Generic: CALL LIN_SOL_SELF (A, B, X [, …])
Specific: The specific interface names are S_LIN_SOL_SELF, D_LIN_SOL_SELF, C_LIN_SOL_SELF, and Z_LIN_SOL_SELF.
Description
Routine LIN_SOL_SELF routine solves a system of linear algebraic equations with a nonsingular coefficient matrix A. By default, the routine computes the factorization of A using Aasen’s method. This decomposition has the form
where P is a permutation matrix, L is a unit lower-triangular matrix, and T is a tridiagonal
self-adjoint matrix. The solution of the linear system Ax = b is found by solving simpler systems,
self-adjoint matrix. The solution of the linear system Ax = b is found by solving simpler systems,
Tv = u
and
More mathematical details for real matrices are found in Golub and Van Loan (1989, Chapter 4).
When the optional Cholesky algorithm is used with a positive definite, self-adjoint matrix, the factorization has the alternate form
where P is a permutation matrix and R is an upper-triangular matrix. The solution of the linear system Ax = b is computed by solving the systems
and
The permutation is chosen so that the diagonal term is maximized at each step of the decomposition. The individual interchanges are optionally available in the argument “pivots”.
Fatal and Terminal Error Messages
See the messages.gls file for error messages for LIN_SOL_SELF. These error messages are numbered 321–336; 341–356; 361–376; 381–396.
Examples
Example 1: Solving a Linear Least-squares System
This example solves a linear least-squares system Cx ≅ d, where Cmxn is a real matrix with m ≥ n. The least-squares solution is computed using the self-adjoint matrix
and the right-hand side
The n × n self-adjoint system Ax = b is solved for x. This solution method is not as satisfactory, in terms of numerical accuracy, as solving the system Cx ≅ d directly by using the routine lin_sol_lsq. Also, see operator_ex05, Chapter 10.
use lin_sol_self_int
use rand_gen_int
implicit none
! This is Example 1 for LIN_SOL_SELF.
integer, parameter :: m=64, n=32
real(kind(1e0)), parameter :: one=1e0
real(kind(1e0)) err
real(kind(1e0)), dimension(n,n) :: A, b, x, res, y(m*n),&
C(m,n), d(m,n)
! Generate two rectangular random matrices.
call rand_gen(y)
C = reshape(y,(/m,n/))
call rand_gen(y)
d = reshape(y,(/m,n/))
! Form the normal equations for the rectangular system.
A = matmul(transpose(C),C)
b = matmul(transpose(C),d)
! Compute the solution for Ax = b.
call lin_sol_self(A, b, x)
! Check the results for small residuals.
res = b - matmul(A,x)
err = maxval(abs(res))/sum(abs(A)+abs(b))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 1 for LIN_SOL_SELF is correct.'
end if
end
Output
Example 1 for LIN_SOL_SELF is correct.
Example 2: System Solving with Cholesky Method
This example solves the same form of the system as Example 1. The optional argument “iopt=” is used to note that the Cholesky algorithm is used since the matrix A is positive definite and self-adjoint. In addition, the sample covariance matrix
is computed, where
the inverse matrix is returned as the “ainv=” optional argument. The scale factor and Γ are computed after returning from the routine. Also, see operator_ex06, Chapter 10.
use lin_sol_self_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 2 for LIN_SOL_SELF.
integer, parameter :: m=64, n=32
real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0
real(kind(1e0)) err
real(kind(1e0)) a(n,n), b(n,1), c(m,n), d(m,1), cov(n,n), x(n,1), &
res(n,1), y(m*n)
type(s_options) :: iopti(1)=s_options(0,zero)
! Generate a random rectangular matrix and a random right hand side.
call rand_gen(y)
c = reshape(y,(/m,n/))
call rand_gen(d(1:n,1))
! Form the normal equations for the rectangular system.
a = matmul(transpose(c),c)
b = matmul(transpose(c),d)
! Use packaged option to use Cholesky decomposition.
iopti(1) = s_options(s_lin_sol_self_Use_Cholesky,zero)
! Compute the solution of Ax=b with optional inverse obtained.
call lin_sol_self(a, b, x, ainv=cov, &
iopt=iopti)
! Compute residuals, x - (inverse)*b, for consistency check.
res = x - matmul(cov,b)
! Scale the inverse to obtain the covariance matrix.
cov = (sum((d-matmul(c,x))**2)/(m-n)) * cov
! Check the results.
err = sum(abs(res))/sum(abs(cov))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 2 for LIN_SOL_SELF is correct.'
end if
end
Output
Example 2 for LIN_SOL_SELF is correct.
Example 3: Using Inverse Iteration for an Eigenvector
This example illustrates the use of the optional argument “iopt=” to reset the value of a Small diagonal term encountered during the factorization. Eigenvalues of the self-adjoint matrix
are computed using the routine lin_eig_self. An eigenvector, corresponding to one of these eigenvalues, λ, is computed using inverse iteration. This solves the near singular system (A – λI)x = b for an eigenvector, x. Following the computation of a normalized eigenvector
the consistency condition
is checked. Since a singular system is expected, suppress the fatal error message that normally prints when the error post-processor routine error_post is called within the routine lin_sol_self. Also, see operator_ex07, Chapter 10.
use lin_sol_self_int
use lin_eig_self_int
use rand_gen_int
use error_option_packet
implicit none
! This is Example 3 for LIN_SOL_SELF.
integer i, tries
integer, parameter :: m=8, n=4, k=2
integer ipivots(n+1)
real(kind(1d0)), parameter :: one=1.0d0, zero=0.0d0
real(kind(1d0)) err
real(kind(1d0)) a(n,n), b(n,1), c(m,n), x(n,1), y(m*n), &
e(n), atemp(n,n)
type(d_options) :: iopti(4)
! Generate a random rectangular matrix.
call rand_gen(y)
c = reshape(y,(/m,n/))
! Generate a random right hand side for use in the inverse
! iteration.
call rand_gen(y(1:n))
b = reshape(y,(/n,1/))
! Compute the positive definite matrix.
a = matmul(transpose(c),c)
! Obtain just the eigenvalues.
call lin_eig_self(a, e)
! Use packaged option to reset the value of a small diagonal.
iopti = d_options(0,zero)
iopti(1) = d_options(d_lin_sol_self_set_small,&
epsilon(one) * abs(e(1)))
! Use packaged option to save the factorization.
iopti(2) = d_options(d_lin_sol_self_save_factors,zero)
! Suppress error messages and stopping due to singularity
! of the matrix, which is expected.
iopti(3) = d_options(d_lin_sol_self_no_sing_mess,zero)
atemp = a
do i=1, n
a(i,i) = a(i,i) - e(k)
end do
! Compute A-eigenvalue*I as the coefficient matrix.
do tries=1, 2
call lin_sol_self(a, b, x, &
pivots=ipivots, iopt=iopti)
! When code is re-entered, the already computed factorization
! is used.
iopti(4) = d_options(d_lin_sol_self_solve_A,zero)
! Reset right-hand side nearly in the direction of the eigenvector.
b = x/sqrt(sum(x**2))
end do
! Normalize the eigenvector.
x = x/sqrt(sum(x**2))
! Check the results.
err = dot_product(x(1:n,1),matmul(atemp(1:n,1:n),x(1:n,1))) - &
e(k)
! If any result is not accurate, quit with no summary printing.
if (abs(err) <= sqrt(epsilon(one))*e(1)) then
write (*,*) 'Example 3 for LIN_SOL_SELF is correct.'
end if
end
Output
Example 3 for LIN_SOL_SELF is correct.
Example 4: Accurate Least-squares Solution with Iterative Refinement
This example illustrates the accurate solution of the self-adjoint linear system
computed using iterative refinement. This solution method is appropriate for least-squares problems when an accurate solution is required. The solution and residuals are accumulated in double precision, while the decomposition is computed in single precision. Also, see operator_ex08, supplied with the product examples.
use lin_sol_self_int
use rand_gen_int
implicit none
! This is Example 4 for LIN_SOL_SELF.
integer i
integer, parameter :: m=8, n=4
real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0
real(kind(1d0)), parameter :: d_zero=0.0d0
integer ipivots((n+m)+1)
real(kind(1e0)) a(m,n), b(m,1), w(m*n), f(n+m,n+m), &
g(n+m,1), h(n+m,1)
real(kind(1e0)) change_new, change_old
real(kind(1d0)) c(m,1), d(m,n), y(n+m,1)
type(s_options) :: iopti(2)=s_options(0,zero)
! Generate a random matrix.
call rand_gen(w)
a = reshape(w, (/m,n/))
! Generate a random right hand side.
call rand_gen(b(1:m,1))
! Save double precision copies of the matrix and right hand side.
d = a
c = b
! Fill in augmented system for accurately solving the least-squares
! problem.
f = zero
do i=1, m
f(i,i) = one
end do
f(1:m,m+1:) = a
f(m+1:,1:m) = transpose(a)
! Start solution at zero.
y = d_zero
change_old = huge(one)
! Use packaged option to save the factorization.
iopti(1) = s_options(s_lin_sol_self_save_factors,zero)
iterative_refinement: do
g(1:m,1) = c(1:m,1) - y(1:m,1) - matmul(d,y(m+1:m+n,1))
g(m+1:m+n,1) = - matmul(transpose(d),y(1:m,1))
call lin_sol_self(f, g, h, &
pivots=ipivots, iopt=iopti)
y = h + y
change_new = sum(abs(h))
! Exit when changes are no longer decreasing.
if (change_new >= change_old) &
exit iterative_refinement
change_old = change_new
! Use option to re-enter code with factorization saved; solve only.
iopti(2) = s_options(s_lin_sol_self_solve_A,zero)
end do iterative_refinement
write (*,*) 'Example 4 for LIN_SOL_SELF is correct.'
end
Output
Example 4 for LIN_SOL_SELF is correct.
Published date: 03/19/2020
Last modified date: 03/19/2020
