Creates the identity matrix.
Identity matrix of size N x N and type real . (Output)
N — Size of output identity matrix. (Input)
EYE (N)
Creates a rank-2 square array whose diagonals are all the value one. The off-diagonals all have value zero.
Dense Matrix Example (operator_ex07.f90)
use linear_operators
implicit none
! This is the equivalent of Example 3 (using operators) for LIN_SOL_SELF.
integer tries
integer, parameter :: m=8, n=4, k=2
integer ipivots(n+1)
real(kind(1d0)) :: one=1.0d0, err
real(kind(1d0)) a(n,n), b(n,1), c(m,n), x(n,1), &
e(n), ATEMP(n,n)
type(d_options) :: iopti(4)
! Generate a random rectangular matrix.
C = rand(C)
! Generate a random right hand side for use in the inverse
! iteration.
b = rand(b)
! Compute the positive definite matrix.
A = C .tx. C; A = (A+.t.A)/2
! Obtain just the eigenvalues.
E = EIG(A)
! Use packaged option to reset the value of a small diagonal.
iopti(4) = 0
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_lin_sol_self_save_factors
! Suppress error messages and stopping due to singularity
! of the matrix, which is expected.
iopti(3) = d_lin_sol_self_no_sing_mess
ATEMP = A
! Compute A-eigenvalue*I as the coefficient matrix.
! Use eigenvalue number k.
A = A - e(k)*EYE(n)
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_lin_sol_self_solve_A
! Reset right-hand side in the direction of the eigenvector.
B = UNIT(x)
end do
! Normalize the eigenvector.
x = UNIT(x)
! Check the results.
b=ATEMP .x. x
err = dot_product(x(1:n,1), b(1:n,1)) - e(k)
! If any result is not accurate, quit with no printing.
if (abs(err) <= sqrt(epsilon(one))*E(1)) then
write (*,*) 'Example 3 for LIN_SOL_SELF (operators) is correct.'
end if
end
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |