MILU_spec.F
c****************************************************************
#include "author.inc"
c* $Id: MILU_spec.F,v 1.1 1995/11/14 02:19:55 turner Exp $
c*
c* Overwrite a matrix in full-storage format by a level zero
c* modified incomplete LU-factorization (fill elements are
c* lumped into diagonal rather than simply being dropped).
c*
c* NOTE: No pivoting is performed.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of a
c* n - actual size of a to be used
c*
c* In/Out:
c* a - matrix to be factored
c*
c* Output:
c* status - return status
c* -4 ==> breakdown
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c****************************************************************
subroutine JT_MILU_Full (idim, n, a, status)
implicit none
c
c ... Input:
integer idim, n
c
c ... In/Out:
real a(idim,n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j, k
real amult, one, zero
c
parameter (zero=0.0d0, one=1.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (n.le.0 .or. idim.le.0 .or. n.gt.idim) then
status = -1
return
endif
c
c ... Test for zero main diagonal element.
do k=1,n
if (a(k,k) .eq. zero) then
status = -4
return
endif
enddo
c
do k=1,n-1
amult = one/a(k,k)
do i=k+1,n
if (a(i,k) .ne. zero) then
a(i,k) = amult*a(i,k)
do j=k+1,n
if (a(k,j) .ne. zero) then
if (a(i,j) .ne. zero) then
a(i,j) = a(i,j) - a(i,k)*a(k,j)
else
a(i,i) = a(i,i) - a(i,k)*a(k,j)
endif
endif
enddo
endif
enddo
enddo
c
return
end
c****************************************************************
#include "author.inc"
c* $Id: MILU_spec.F,v 1.1 1995/11/14 02:19:55 turner Exp $
c*
c* Overwrite a matrix in ELL storage format by a level zero
c* modified incomplete LU-factorization (fill elements are
c* lumped into diagonal rather than simply being dropped).
c*
c* Note that the main difference between this routine and
c* JT_ILU_ELL is that the j loop goes from i+1 to n rather
c* than from 1 to maxnz. This increases cost significantly
c* for large systems. Is there a better way?
c*
c* NOTE: No pivoting is performed.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of array a
c* n - actual size of a to be used
c* maxnz - maximum number of non-zero elements in any row
c* ja - column map for array a
c*
c* In/Out:
c* a - matrix to be factored
c*
c* Output:
c* status - return status
c* -4 ==> breakdown
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c****************************************************************
subroutine JT_MILU_ELL (idim, n, maxnz, ja, a, status)
implicit none
c
c ... Input:
integer idim, maxnz, n
integer ja(idim,maxnz)
c
c ... In/Out:
real a(idim,maxnz)
c
c ... Output:
integer status
c
c ... Local:
integer i, ij, j, k, kj, kk
real one, zero, amult
c
parameter (zero=0.0d0, one=1.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (n.le.0 .or. idim.le.0 .or. maxnz.le.0 .or. n.gt.idim) then
status = -1
return
endif
c
c ... Test for zero main diagonal element.
do k=1,n
if (a(k,1) .eq. zero) then
status = -4
return
endif
enddo
c
do k=1,n-1
amult = one/a(k,1)
c
c .... Loop over rows k+1 to n.
do i=k+1,n
c
c ..... Find kk such that ja(i,kk) = k. Note that if one is found,
c we also know that a(i,kk) is non-zero.
do kk=2,maxnz
if (ja(i,kk) .eq. 0) then
c
c ....... Out of elements, so a(i,kk) = 0, so bail out and get
c another i.
goto 35
elseif (ja(i,kk) .eq. k) then
c
c ....... Found kk.
goto 15
endif
enddo
c
c ..... Didn't find kk, so a(i,kk) = 0, so bail out and get another i.
goto 35
c
15 continue
a(i,kk) = amult*a(i,kk)
c
c ..... Loop over columns k+1 to n.
do j=k+1,n
c
c ...... Find kj such that ja(k,kj) = j. Note that if one is found,
c we also know that a(k,kj) is non-zero.
do kj=2,maxnz
if (ja(k,kj) .eq. 0) then
c
c ........ Out of elements, so a(k,kj) = 0, so bail out and get
c another j.
goto 30
elseif (ja(k,kj) .eq. j) then
c
c ........ Found kj.
goto 20
endif
enddo
c
c ...... Didn't find kj, so a(k,kj) = 0, so bail out and get another j.
goto 30
c
c ...... Find ij such that ja(i,ij) = j. Whether one is found or
c not will determine whether to modify the diagonal or the
c off-diagonal.
20 continue
do ij=2,maxnz
if (ja(i,ij) .eq. 0) then
c
c ........ Out of elements, so a(i,ij) = 0, so bail out and modify
c diagonal.
goto 25
elseif (ja(i,ij) .eq. j) then
c
c ........ Found ij, so modify element a(i,ij).
a(i,ij) = a(i,ij) - a(i,kk)*a(k,kj)
goto 30
endif
enddo
c
c ...... Didn't find ij, so a(i,ij) = 0, so modify diagonal.
25 continue
a(i,1) = a(i,1) - a(i,kk)*a(k,kj)
30 continue
enddo ! end of loop over columns (j)
35 continue
enddo ! end of loop over rows (i)
enddo ! end of main loop (k)
c
return
end