ILU_spec.F
c****************************************************************
#include "author.inc"
c* $Id: ILU_spec.F,v 1.1 1995/11/14 02:19:50 turner Exp $
c*
c* Overwrite a matrix in full-storage format by a level zero
c* incomplete LU-factorization.
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_ILU_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 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. 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)
endif
endif
enddo
endif
enddo
enddo
c
return
end
c****************************************************************
#include "author.inc"
c* $Id: ILU_spec.F,v 1.1 1995/11/14 02:19:50 turner Exp $
c*
c* Overwrite a matrix in ELL storage format by a level zero
c* incomplete LU-factorization.
c*
c* NOTE: No pivoting is performed.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of arrays
c* n - number of rows to use
c* maxnz - maximum number of non-zero elements in any row
c* ja - column map for array a
c*
c* In/Out:
c* a - original matrix on input, factored matrix on output
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_ILU_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, k, kj, kk
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. 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 30
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 30
c
15 continue
a(i,kk) = amult*a(i,kk)
c
c ..... Loop over *all* columns (not just 2:maxnz).
do ij=1,maxnz
c
c ...... Only want to consider columns k+1 to n, so if ja(i,ij) < k+1
c get another ij. This also automatically ensures we get only
c the non-zero elements of a, since zero elements have
c ja(i,ij) = 0.
if (ja(i,ij) .ge. k+1) then
c
c ....... Find kj such that ja(k,kj) = ja(i,ij). Again, if one is found,
c we 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 ij.
goto 25
elseif (ja(k,kj) .eq. ja(i,ij)) then
c
c ......... Found kj.
goto 20
endif
enddo
c
c ....... Didn't find one, so a(k,kj) = 0, so bail out and get
c another ij.
goto 25
c
c ....... Finally ready to do the update.
20 continue
a(i,ij) = a(i,ij) - a(i,kk)*a(k,kj)
c
endif
25 continue
enddo ! end of loop over columns (ij)
30 continue
enddo ! end of loop over rows (i)
enddo ! end of main loop (k)
c
return
end