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.