LIN_SOL_GEN

   more...

Solves a general system of linear equations Ax = b. Using optional arguments, any of several related computations can be performed. These extra tasks include computing the LU factorization of A using partial pivoting, representing the determinant of A, computing the inverse matrix A-1, and solving ATx = b or Ax = b given the LU factorization of A.

Required Arguments

A — Array of size n × n containing the matrix. (Input [/Output])
If the packaged option lin_sol_gen_save_LU is used then the LU factorization of A is saved in A. For solving efficiency, the diagonal reciprocals of the matrix U are saved in the diagonal entries of A.

B — Array of size n × nb containing the right-hand side matrix. (Input [/Output])
If the packaged option lin_sol_gen_save_LU 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)

NRHS = nb (Input)
Uses 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 that contains the individual row interchanges. To construct the permuted order so that no pivoting is required, define an integer array ip(n). Initialize ip(i) = i, i = 1, n and then execute the loop, after calling lin_sol_gen,

k=pivots(i)
interchange ip(i) and ip(k), i=1,n

The matrix defined by the array assignment that permutes the rows, A(1:n, 1:n) = A(ip(1:n), 1:n), requires no pivoting for maintaining numerical stability. Now, the optional argument “iopt=” and the packaged option number ?_lin_sol_gen_no_pivoting can be safely used for increased efficiency during the LU factorization of A.

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 this expression is given by 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.

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:

iopt(IO) = ?_options(?_lin_sol_gen_set_small, Small)
Replaces a diagonal term of the matrix U if it is smaller in magnitude than the value Small using the same sign or complex direction as the diagonal. The system is declared singular. A solution is approximated based on this replacement if no overflow results.
Default: the smallest number that can be reciprocated safely

iopt(IO) = ?_options(?_lin_sol_gen_save_LU, ?_dummy)
Saves the LU factorization of A. Requires the optional argument “pivots=” if the routine will be used later for solving 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 U are saved in the diagonal entries of A.

iopt(IO) = ?_options(?_lin_sol_gen_solve_A, ?_dummy)
Uses the LU factorization of A computed and saved to solve Ax = b.

iopt(IO) = ?_options(?_lin_sol_gen_solve_ADJ, ?_dummy)
Uses the LU factorization of A computed and saved to solve ATx = b.

iopt(IO) = ?_options(?_lin_sol_gen_no_pivoting, ?_dummy)
Does no row pivoting. The array pivots (:), if present, are output as pivots (i) = i, for i = 1, …, n.

iopt(IO) = ?_options(?_lin_sol_gen_scan_for_NaN, ?_dummy)
Examines each input array entry to find the first value such that

isNaN(a(i,j)) .or. isNan(b(i,j)) == .true.

See the isNaN() function, Chapter 10.
Default: Does not scan for NaNs.

iopt(IO) = ?_options(?_lin_sol_gen_no_sing_mess, ?_dummy)
Do not output an error message when the matrix A is singular.

iopt(IO) = ?_options(?_lin_sol_gen_A_is_sparse, ?_dummy)
Uses an indirect updating loop for the LU factorization that is efficient for sparse matrices where all matrix entries are stored.

FORTRAN 90 Interface

Generic: CALL LIN_SOL_GEN (A, B, X [, …])

Specific: The specific interface names are S_LIN_SOL_GEN, D_LIN_SOL_GEN, C_LIN_SOL_GEN, and Z_LIN_SOL_GEN.

Description

Routine LIN_SOL_GEN solves a system of linear algebraic equations with a nonsingular coefficient matrix A. It first computes the LU factorization of A with partial pivoting such that LU = A. The matrix U is upper triangular, while the following is true:

The factors Pi and Li are defined by the partial pivoting. Each Pi is an interchange of row i with row j ≥ i. Thus, Pi is defined by that value of j. Every

is an elementary elimination matrix. The vector is zero in entries 1, …, i. This vector is stored as column i in the strictly lower-triangular part of the working array containing the decomposition information. The reciprocals of the diagonals of the matrix U are saved in the diagonal of the working array. The solution of the linear system Ax = b is found by solving two simpler systems,

y = L-1b and x = U-1y

More mathematical details are found in Golub and Van Loan (1989, Chapter 3).

Fatal and Terminal Error Messages

See the messages.gls file for error messages for LIN_SOL_GEN. The messages are numbered 161–175; 181–195; 201–215; 221–235.

Examples

Example 1: Solving a Linear System of Equations

This example solves a linear system of equations. This is the simplest use of lin_sol_gen. The equations are generated using a matrix of random numbers, and a solution is obtained corresponding to a random right-hand side matrix. Also, see operator_ex01, supplied with the product examples, for this example using the operator notation.

      use lin_sol_gen_int 
      use rand_gen_int 
      use error_option_packet 
 
      implicit none 
 
! This is Example 1 for LIN_SOL_GEN.  
 
      integer, parameter :: n=32 
      real(kind(1e0)), parameter :: one=1e0 
      real(kind(1e0)) err 
      real(kind(1e0)) A(n,n), b(n,n), x(n,n), res(n,n), y(n**2) 
 
! Generate a random matrix. 
      call rand_gen(y) 
 
 
! Generate random right-hand sides. 
      call rand_gen(y) 
      b = reshape(y,(/n,n/)) 
 
! Compute the solution matrix of Ax=b. 
      call lin_sol_gen(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_GEN is correct.' 
      end if 
 
      end

Output

Example 1 for LIN_SOL_GEN is correct.

Example 2: Matrix Inversion and Determinant

This example computes the inverse and determinant of A, a random matrix. Tests are made on the conditions

and

Also, see operator_ex02.

use lin_sol_gen_int

use rand_gen_int

implicit none

! This is Example 2 for LIN_SOL_GEN.

integer i

integer, parameter :: n=32

real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0

real(kind(1e0)) err

real(kind(1e0)) A(n,n), b(n,0), inv(n,n), x(n,0), res(n,n), &

y(n**2), determinant(2), inv_determinant(2)

! Generate a random matrix.

call rand_gen(y)

A = reshape(y,(/n,n/))

! Compute the matrix inverse and its determinant.

call lin_sol_gen(A, b, x, nrhs=0, &

ainv=inv, det=determinant)

! Compute the determinant for the inverse matrix.

call lin_sol_gen(inv, b, x, nrhs=0, &

det=inv_determinant)

! Check residuals, A times inverse = Identity.

res = matmul(A,inv)

do i=1, n

res(i,i) = res(i,i) - one

end do

err = sum(abs(res)) / sum(abs(a))

if (err <= sqrt(epsilon(one))) then

if (determinant(1) == inv_determinant(1) .and. &

(abs(determinant(2)+inv_determinant(2)) &

<= abs(determinant(2))*sqrt(epsilon(one)))) then

write (*,*) 'Example 2 for LIN_SOL_GEN is correct.'

end if

end if

end

Output

Example 2 for LIN_SOL_GEN is correct.

Example 3: Solving a System with Iterative Refinement

This example computes a factorization of a random matrix using single-precision arithmetic. The double-precision solution is corrected using iterative refinement. The corrections are added to the developing solution until they are no longer decreasing in size. The initialization of the derived type array iopti(1:2) = s_option(0,0.0e0) leaves the integer part of the second element of iopti(:) at the value zero. This stops the internal processing of options inside lin_sol_gen. It results in the LU factorization being saved after exit. The next time the routine is entered the integer entry of the second element of iopt(:) results in a solve step only. Since the LU factorization is saved in arrays A(:,:) and ipivots(:), at the final step, solve only steps can occur in subsequent entries to lin_sol_gen. Also, see operator_ex03, Chapter 10.

use lin_sol_gen_int

use rand_gen_int

implicit none

! This is Example 3 for LIN_SOL_GEN.

integer, parameter :: n=32

real(kind(1e0)), parameter :: one=1.0e0, zero=0.0e0

real(kind(1d0)), parameter :: d_zero=0.0d0

integer ipivots(n)

real(kind(1e0)) a(n,n), b(n,1), x(n,1), w(n**2)

real(kind(1e0)) change_new, change_old

real(kind(1d0)) c(n,1), d(n,n), y(n,1)

type(s_options) :: iopti(2)=s_options(0,zero)

! Generate a random matrix.

call rand_gen(w)

a = reshape(w, (/n,n/))

! Generate a random right hand side.

call rand_gen(b(1:n,1))

! Save double precision copies of the matrix and right hand side.

d = a

c = b

! 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_gen_save_LU,zero)

iterative_refinement: do

b = c - matmul(d,y)

call lin_sol_gen(a, b, x, &

pivots=ipivots, iopt=iopti)

y = x + y

change_new = sum(abs(x))

! 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_gen_solve_A,zero)

end do iterative_refinement

write (*,*) 'Example 3 for LIN_SOL_GEN is correct.'

end

Output

Example 3 for LIN_SOL_GEN is correct.

Example 4: Evaluating the Matrix Exponential

This example computes the solution of the ordinary differential equation problem

with initial values y(0) = y0. For this example, the matrix A is real and constant with respect to t. The unique solution is given by the matrix exponential:

This method of solution uses an eigenvalue-eigenvector decomposition of the matrix

to evaluate the solution with the equivalent formula

where

is computed using the complex arithmetic version of lin_sol_gen. The results for y(t) are real quantities, but the evaluation uses intermediate complex-valued calculations. Note that the computation of the complex matrix X and the diagonal matrix D is performed using the IMSL MATH/LIBRARY FORTRAN 77 interface to routine EVCRG. This is an illustration of intermixing interfaces of FORTRAN 77 and Fortran 90 code. The information is made available to the Fortran 90 compiler by using the FORTRAN 77 interface for EVCRG. Also, see operator_ex04, supplied with the product examples, where the Fortran 90 function EIG() has replaced the call to EVCRG.

use lin_sol_gen_int

use rand_gen_int

use Numerical_Libraries

implicit none

! This is Example 4 for LIN_SOL_GEN.

integer, parameter :: n=32, k=128

real(kind(1e0)), parameter :: one=1.0e0, t_max=1, delta_t=t_max/(k-1)

real(kind(1e0)) err, A(n,n), atemp(n,n), ytemp(n**2)

real(kind(1e0)) t(k), y(n,k), y_prime(n,k)

complex(kind(1e0)) EVAL(n), EVEC(n,n)

complex(kind(1e0)) x(n,n), z_0(n,1), y_0(n,1), d(n)

integer i

! Generate a random matrix in an F90 array.

call rand_gen(ytemp)

atemp = reshape(ytemp,(/n,n/))

! Assign data to an F77 array.

A = atemp

! Use IMSL Numerical Libraries F77 subroutine for the

! eigenvalue-eigenvector calculation.

CALL EVCRG(N, A, N, EVAL, EVEC, N)

! Generate a random initial value for the ODE system.

call rand_gen(ytemp(1:n))

y_0(1:n,1) = ytemp(1:n)

! Assign the eigenvalue-eigenvector data to F90 arrays.

d = EVAL; x = EVEC

! Solve complex data system that transforms the initial values, Xz_0=y_0.

call lin_sol_gen(x, y_0, z_0)

t = (/(i*delta_t,i=0,k-1

! Compute y and y' at the values t(1:k).

y = matmul(x, exp(spread(d,2,k)*spread(t,1,n))* &

spread(z_0(1:n,1),2,k))

y_prime = matmul(x, spread(d,2,k)* &

exp(spread(d,2,k)*spread(t,1,n))* &

spread(z_0(1:n,1),2,k))

! Check results. Is y' - Ay = 0?

err = sum(abs(y_prime-matmul(atemp,y))) / &

(sum(abs(atemp))*sum(abs(y)))

if (err <= sqrt(epsilon(one))) then

write (*,*) 'Example 4 for LIN_SOL_GEN is correct.'

end if

end

Output

Example 4 for LIN_SOL_GEN is c orrect.