c ****************************************************************** c John Kerl c 8/18/93 c ****************************************************************** subroutine d_lusolve(A, x, b, z, L, U, n, ier) implicit none integer n, ier double precision A(n,n), x(n), b(n), L(n,n), U(n,n), z(n) logical debug debug=.true. call d_ludecomp(A, L, U, n, ier) if (ier .eq. 1) then write(0, *) "d_lusolve: Matrix couldn't be factored." ier=1 return end if if (debug) then write (*,*) 'A is:' call d_write_matrix(A, n, n, '', 0, 'g', -1) write (*,*) 'L is:' call d_write_matrix(L, n, n, '', 0, 'lt', -1) write (*,*) 'U is:' call d_write_matrix(U, n, n, '', 0, 'ut', -1) end if call d_backsubplus(l, b, z, n, ier) if (ier .eq. 1) then write(0, *) &"d_lusolve: Could not back-substitute solution (divide by zero)." ier=1 return end if call d_backsubminus(u, z, x, n, ier) if (ier .eq. 1) then write(0, *) &"d_lusolve: Could not back-substitute solution (divide by zero)." ier=1 return end if ier=0 end ! subroutine d_lusolve c ****************************************************************** subroutine d_ludecomp(A, L, U, n, ier) ! Need to pack L,U into one mat implicit none integer n, ier double precision A(n,n), L(n,n), U(n,n) integer i, j, k c Initialize workspace to the identity matrix do i=1,n do j=1,n L(i,j)=0.0 U(i,j)=0.0 end do end do do i=1,n L(i,i)=1.0 end do U(1,1)=A(1,1) if (U(1,1)*L(1,1) .eq. 0.0) then write (0,*) "d_ludecomp: Factorization impossible." ier=1 return end if do j=2,n U(1,j) = A(1,j) / L(1,1) L(j,1) = A(j,1) / U(1,1) end do do i=2, n-1 U(i,i) = A(i,i) do k=1, i-1 U(i,i) = U(i,i) - L(i,k) * U(k,i) end do ! k=1, i-1 if (L(i,i)*U(i,i) .eq. 0.0) then write (0,*) "d_ludecomp: Factorization impossible" ier = 1 return end if do j=i+1, n U(i,j) = A(i,j) do k=1, i-1 U(i,j) = U(i,j) - L(i,k) * U(k,j) end do U(i,j) = U(i,j) / L(i,i) L(j,i) = A(j,i) do k=1, i-1 L(j,i) = L(j,i) - L(j,k) * U(k,i) end do L(j,i) = L(j,i) / U(i,i) end do ! j=i+1, n end do ! i=2, n-1 U(n,n) = A(n,n) do k=1, n-1 U(n,n) = U(n,n) - L(n,k)*U(k,n) end do ier=0 end ! subroutine d_ludecomp