program main !==================================================================== ! Evaluate the condition number of matrix A ! Method: straightforward evaluation !==================================================================== implicit none integer, parameter :: n=3 double precision a(n,n), CondA, Mcond integer i,j ! matrix A ! data (a(1,i), i=1,3) / 3.0, 2.0, 4.0 / ! data (a(2,i), i=1,3) / 2.0, -3.0, 1.0 / ! data (a(3,i), i=1,3) / 1.0, 1.0, 2.0 / data (a(1,i), i=1,3) / 3.0, 2.0, 4.0 / data (a(2,i), i=1,3) / 3.000001, 2.0, 4.000002 / data (a(3,i), i=1,3) / 1.0, 1.0, 2.0 / ! print a header and the original matrix write (*,200) do i=1,n write (*,201) (a(i,j),j=1,n) end do CondA = MCond(A,n) ! print the result write(*,202) CondA 200 format (' Evaluate the condition number of matrix A ',/,/, & ' Matrix A') 201 format (6f12.6) 202 format (' The condition number C(A) = ',1pe12.4) end program main double precision function MCond(A,n) !============================================================ ! Evaluate the condition number of matrix A ! Method: C(A) = norm(A*A^{-1}) ! Alex G. December 2009 !----------------------------------------------------------- ! input ... ! a(n,n) - array of coefficients for matrix A ! n - dimension ! output ... ! Matnorm- the condition number ! comments ... calls !=========================================================== implicit none integer n double precision a(n,n), c(n,n) double precision Mnorm integer i,j MCond = Mnorm(A,n) call inverse(a,c,n) MCond = MCond*Mnorm(C,n) end function MCond double precision function Mnorm(a,n) !=========================================================== ! evaluate the Euclidean norm of matrix A ! Alex G. December 2009 !----------------------------------------------------------- ! input ... ! a(n,n) - array of coefficients for matrix A ! n - dimension ! output ... ! Mnorm - the Euclidean norm !=========================================================== implicit none integer n double precision a(n,n) double precision Mnorn integer i, j Mnorm = 0.0 do i=1,n do j=1,n Mnorm = Mnorm + a(i,j)**2 end do end do Mnorm = sqrt(Mnorm) end function Mnorm subroutine inverse(a,c,n) !============================================================ ! Inverse matrix ! Method: Based on Doolittle LU factorization for Ax=b ! Alex G. December 2009 !----------------------------------------------------------- ! input ... ! a(n,n) - array of coefficients for matrix A ! n - dimension ! output ... ! c(n,n) - inverse matrix of A ! comments ... ! the original matrix a(n,n) will be destroyed ! during the calculation !=========================================================== implicit none integer n double precision a(n,n), c(n,n) double precision L(n,n), U(n,n), b(n), d(n), x(n) double precision coeff integer i, j, k ! step 0: initialization for matrices L and U and b ! Fortran 90/95 aloows such operations on matrices L=0.0 U=0.0 b=0.0 ! step 1: forward elimination do k=1, n-1 do i=k+1,n coeff=a(i,k)/a(k,k) L(i,k) = coeff do j=k+1,n a(i,j) = a(i,j)-coeff*a(k,j) end do end do end do ! Step 2: prepare L and U matrices ! L matrix is a matrix of the elimination coefficient ! + the diagonal elements are 1.0 do i=1,n L(i,i) = 1.0 end do ! U matrix is the upper triangular part of A do j=1,n do i=1,j U(i,j) = a(i,j) end do end do ! Step 3: compute columns of the inverse matrix C do k=1,n b(k)=1.0 d(1) = b(1) ! Step 3a: Solve Ld=b using the forward substitution do i=2,n d(i)=b(i) do j=1,i-1 d(i) = d(i) - L(i,j)*d(j) end do end do ! Step 3b: Solve Ux=d using the back substitution x(n)=d(n)/U(n,n) do i = n-1,1,-1 x(i) = d(i) do j=n,i+1,-1 x(i)=x(i)-U(i,j)*x(j) end do x(i) = x(i)/u(i,i) end do ! Step 3c: fill the solutions x(n) into column k of C do i=1,n c(i,k) = x(i) end do b(k)=0.0 end do end subroutine inverse