Matrix Optional Data Changes
To reset tolerances for determining singularity and to allow for other data changes, non-allocated “hidden” variables are defined within the modules. These variables can be allocated first, then assigned values which result in the use of different tolerances or greater efficiency in the executable program. The non-allocated variables, whose scope is limited to the module, are hidden from the casual user. Default values or rules are applied if these arrays are not allocated. In more detail, the inverse matrix operator “.i.” applied to a square matrix first uses the LU factorization code LIN_SOL_GEN and row pivoting. The default value for a small diagonal term is defined to be:
sqrt(epsilon(A))*sum(abs(A))/(n*n+1)
If the system is singular, a generalized matrix inverse is computed with the QR factorization code LIN_SOL_LSQ using this same tolerance. Both row and column pivoting are used. If the system is singular, an error message will be printed and a Fortran 90 STOP is executed. Users may want to change this rule. This is illustrated by continuing and not printing the error message. The following is a additional source to accomplish this, for all following invocations of the operator “.i.”:
allocate(s_inv_options(1))
s_inv_options (1) = skip_error_processing
B = .i. A
There are additional self-documenting integer parameters, packaged in the module linear_operators, that allow users other choices, such as changing the value of the tolerance, as noted above. Included is the ability to have the option apply for just the next invocation of the operator. Options are available that allow optional data to be passed to supporting Fortran 90 subroutines. This is illustrated in the following example:
Operator_ex36.f90
use linear_operators
implicit none
! This is the equivalent of Example 4 for LIN_GEIG_GEN (using operators).
integer, parameter :: n=32
real(kind(1d0)), parameter :: one=1d0, zero=0d0
real(kind(1d0)) a(n,n), b(n,n), bta(n), err
complex(kind(1d0)) alpha(n), v(n,n)
! Generate random matrices for both A and B.
A = rand(A); B = rand(B)
! Set the option, a larger tolerance than default for lin_sol_lsq.
allocate(d_eig_options(6))
d_eig_options(1) = options_for_lin_geig_gen
d_eig_options(2) = 4
d_eig_options(3) = d_lin_geig_gen_for_lin_sol_lsq
d_eig_options(4) = 2
d_eig_options(5) = d_options(d_lin_sol_lsq_set_small,&
sqrt(epsilon(one))*norm(B,1))
d_eig_options(6) = d_lin_sol_lsq_no_sing_mess
! Compute the generalized eigenvalues.
alpha = EIG(A, B=B, D=bta, W=V)
! Check the residuals.
err = norm((A .x. V .x. diag(bta)) - (B .x. V .x. diag(alpha)),1)/&
(norm(A,1)*norm(bta,1)+norm(B,1)*norm(alpha,1))
if (err <= sqrt(epsilon(one))) then
write (*,*) 'Example 4 for LIN_GEIG_GEN (operators) is correct.'
end if
! Clean up the allocated array. This is good housekeeping.
deallocate(d_eig_options)
end
Note that in this example one first allocates the array by which the user will pass the new options for EIG to use. This array is named d_eig_options in accordance with the name of the unallocated option array specified in the documentation for EIG. A size of 6 is specified because a total of six options must be passed to EIG to accomplish the resetting of the singular value tolerance and to turn off the printing of the error message when the matrix is singular. The first entry of d_eig_options specifies which of the options for EIG will be set. The next entry designates the number of entries which follows that apply to “options_for_lin_geig_gen”. The third entry specifies the option value of LIN_GEIG_GEN to be set, d_lin_geig_gen_for_lin_sol_lsq. The fourth entry specifies the number of entries that follow which apply to LIN_SOL_LSQ. Finally, the fifth and sixth entries set the two LIN_SOL_LSQ options that we desire.
