Skip to content

Latest commit

 

History

History
185 lines (137 loc) · 6.86 KB

inverse_iteration.md

File metadata and controls

185 lines (137 loc) · 6.86 KB

Find Eigenvalues of a Square Linear System Using Inverse Iteration (LU Decomposition)

This is a part of the student software manual project for Math 5610: Computational Linear Algebra and Solution of Systems of Equations.

Routine Name: inverse_iteration

Author: Christian Bolander

Language: Fortran. This code can be compiled using the GNU Fortran compiler by

$ gfortran -c inverse_iteration.f90

and can be added to a program using

$ gfortran program.f90 inverse_iteration.o

Description/Purpose: This subroutine uses the inverse iteration method with shifting to iteratively find any eigenvalue and the corresponding eigenvector given an appropriate shift (a shift of zero produces the smallest eigenvalue). The idea is that if the eigenvalues of A are , then the eigenvalues of A - I are - and the eigenvalues of B = (A - I) are

Meaning that a value can be found in the inverse problem that satisfies a given eigenvalue in the original problem.

Input:

A : REAL - a coefficient matrix of size n x n

n : INTEGER - the rank of A and the length of the eigenvectors

v0 : REAL - the initial guess of the eigenvetor

alpha : REAL - a shift given to the A matrix

tol : REAL - the exit tolerance for the algorithm

maxiter : INTEGER - the maximum number of iterations before exit

printit : INTEGER - flag for printing final convergence information

Output:

v : REAL - the final approximation of the eigenvector corresponding to the largest eigenvalue

lam0 : REAL - the final approximation of the largest eigenvalue

Usage/Example:

This routine can be implemented in a program as follows

INTEGER :: n, maxiter, i, printit
REAL*8, ALLOCATABLE :: A(:, :), v0(:), v(:)
REAL*8 :: tol, lam, alpha

n = 3
ALLOCATE(A(n, n), v0(n), v(n))
lam = 0.D0
tol = 10.D-16
maxiter = 10000
printit = 1
v0 = 1.D0
alpha = 0.0D0
A = RESHAPE((/1.D0, 2.D0, 0.D0, &
			& 2.D0, 1.D0, 2.D0, &
            & 0.D0, 2.D0, 1.D0/), (/n, n/), ORDER=(/2, 1/))
CALL inverse_iteration(A, n, v0, alpha, tol, maxiter, printit, v, lam)
WRITE(*,*) lam, v

The outputs from the above code:

   8.8817841970012523E-016          28  ! Convergence error at exit and the exit iterations
   -1.8284271247461898                   ! Approximation of largest eigenvalue
   0.70710678052474663   -1.0000000000000000   0.70710679885817773  ! Normalized eigenvector

Additionally, using this routine on an 20 x 20 Hilbert matrix can be done with the following code:

INTEGER :: n, maxiter, i, printit, j
REAL*8, ALLOCATABLE :: A(:, :), v0(:), v(:)
REAL*8 :: tol, lam, alpha

n = 20
ALLOCATE(A(n, n), v0(n), v(n))
lam = 0.D0
tol = 10.D-16
maxiter = 10000
printit = 1
v0 = 1.D0
alpha = 0.D0
DO i = 1, n
	DO j = 1, n
		A(i, j) = 1.0D0/(REAL(i) + REAL(j) - 1.0D0)
	END DO
END DO
CALL inverse_iteration(A, n, v0, tol, alpha, maxiter, printit, v, lam)
WRITE(*,*) lam, v

and the output is

   8.4534686233368191E-020          31   ! Convergence error and number of iterations
   -1.8190197007123243E-018                    ! Eigenvalue
   ! Eigenvector
   7.1434857804492041E-009  -1.1494603201052072E-006   4.5632987553809310E-005  -7.7840548626109121E-004   7.0518472705676013E-003  -3.7471861751792644E-002  0.12163121197985835      -0.23992887217693432       0.27453697897244583      -0.20088876891016452       0.30726509218971626      -0.74604545697557068        1.0000000000000000      -0.62497068848836279        2.3580846900107453E-002  0.35574384095812017      -0.54925444346690810       0.52524084608927935      -0.27334781169413414        5.7591158369571113E-002

Implementation/Code: The code for inverse_iteration can be seen below.

SUBROUTINE inverse_iteration(A, n, v0, alpha, tol, maxiter, printit, v, lam0)
	IMPLICIT NONE
	
	! Implements the inverse iteration method for finding any eigenvalue
	! in the system (provided an accurate guess of `alpha` is given).
	
	! Takes as an input the matrix `A` of rank `n` that contains the
	! system to be analyzed, an initial guess for the eigenvector, `v0`,
	! the shift, `alpha` that is to be used to find the eigenvalue,
	! a tolerance, `tol`, for exiting the iterative solver as well as a
	! maximum number of iterations, `maxiter`. Finally, a flag to print
	! the final number of iterations and convergence error, `printit` is
	! an input. The output of the subroutine is the final eigenvector,
	! `v1` produced from the algorithm (scaled by the element with the
	! maximum absolute value) as well as the final eigenvalue, `lam0`.
	! If `alpha` is zero, the minimum eigenvalue will be found.
	INTEGER, INTENT(IN) :: n, maxiter, printit
	REAL*8, INTENT(IN) :: A(n, n), v0(n), alpha, tol
	REAL*8, INTENT(OUT) :: v(n), lam0
	
	REAL*8 :: v1(n), lam1, error, norm, ALU(n, n)
	INTEGER :: i, k
	
	! Initializes variables
	v = 0.D0
	lam1 = 0.D0
	error = 10.D0*tol
	norm = 0.D0
	k = 0
	lam0 = 0.D0
	v1 = v0
	ALU = A
	
	! Shift the `A` matrix.
	DO i = 1, n
		ALU(i, i) = ALU(i, i) - alpha
	END DO
	
	! Perform an LU factorization on the shifted matrix for efficiency.
	CALL lu_factor(ALU, n)
	
	! Iterate until the error or number of iterations reaches the given
	! limits
	DO WHILE (error > tol .AND. k < maxiter)
		! Solve the LU equation with the previous eigenvalue guess
		CALL lu_solve(ALU, n, v1, v)
		
		! Normalize the newest guess
		CALL l2_vec_norm(v, n, norm)
		v1 = v/norm
		
		! Calculate the new guess for the eigenvalue
		CALL mat_prod(A, v1, n, n, 1, v)
		CALL vec_dot_prod(v1, v, n, lam1)
		
		! Calculate convergence error and increment values
		error = DABS(lam1 - lam0)
		lam0 = lam1
		k = k + 1
	END DO
	
	! Normalizes the eigenvector according to the maximum absolute value
	! of the elements.
	v = v1/MAXVAL(DABS(v1))
	
	! Prints the error and number of iterations when exiting.
	IF (printit == 1) THEN
		WRITE(*,*) error, k
	END IF
END SUBROUTINE