SolveUTriang_spec.F
c**********************************************************************
#include "author.inc"
c* $Id: SolveUTriang_spec.F,v 1.1 1995/11/14 02:19:45 turner Exp $
c*
c* Solves upper triangular system in full conventional format.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idiag - indicates whether diagonal is unity
c* 0 => general case
c* non-zero => unit diagonal
c* itrans - indicates whether to use transpose
c* 0 => no
c* non-zero => use transpose
c* idim - leading dimension of a
c* n - number of unknowns
c* a - matrix
c*
c* In/Out:
c* x - source vector on input, solution on output
c*
c* Output:
c* status - return status
c* -4 ==> breakdown
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_y_eq_y_plus_sx
c*
#include "copyright.inc"
c**********************************************************************
subroutine JT_SolveUTriang_Full (idiag, itrans, idim, n, a,
& x, status)
implicit none
c
c ... Input:
integer idiag, itrans, idim, n
real a(idim,n)
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
#ifdef use_lapack
integer info
#else
integer i, j
real zero
c
parameter (zero=0.0d0)
#endif
c
c ... Set return status.
status = 0
c
c ... Check arguments.
if (idim.le.0 .or. n.le.0 .or. n.gt.idim) then
status = -1
return
endif
c
#ifdef use_lapack
if (itrans .eq. 0) then
c
c .... Use coeff. as is.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal.
call STRTRS ('U', 'N', 'N', n, 1, a, idim, x, n, info)
else
c
c ..... Unit diagonal.
call STRTRS ('U', 'N', 'U', n, 1, a, idim, x, n, info)
endif
else
c
c .... Use transpose.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal.
call STRTRS ('L', 'T', 'N', n, 1, a, idim, x, n, info)
else
c
c ..... Unit diagonal.
call STRTRS ('L', 'T', 'U', n, 1, a, idim, x, n, info)
endif
endif
if (info .lt. 0) then
status = -1
elseif (info .gt. 0) then
status = -4
endif
#else
c
c ... Check diagonal elements.
if (idiag .eq. 0) then
do i=1,n
if (a(i,i) .eq. zero) then
status = -4
return
endif
enddo
endif
c
if (itrans .eq. 0) then
c
c .... Use coeff. as is.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal (used in ILU preconditioning).
do j=n,2,-1
x(j) = x(j) / a(j,j)
call JT_y_eq_y_plus_sx (j-1, -x(j), a(1,j), x(1), status)
enddo
x(1) = x(1) / a(1,1)
else
c
c ..... Unit diagonal (UNUSED, UNTESTED).
do j=n,2,-1
call JT_y_eq_y_plus_sx (j-1, -x(j), a(1,j), x(1), status)
enddo
endif
else
c
c .... Use transpose.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal (used in IC preconditioning).
x(n) = x(n) / a(n,n)
do i=n-1,1,-1
do j=i+1,n
x(i) = x(i) - a(j,i)*x(j)
enddo
x(i) = x(i) / a(i,i)
enddo
else
c
c ..... Unit diagonal (UNUSED, UNTESTED).
do i=n-1,1,-1
do j=i+1,n
x(i) = x(i) - a(j,i)*x(j)
enddo
enddo
endif
endif
#endif
c
return
end
c**********************************************************************
#include "author.inc"
c* $Id: SolveUTriang_spec.F,v 1.1 1995/11/14 02:19:45 turner Exp $
c*
c* Solves upper triangular system stored in ELL format.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idiag - indicates whether diagonal is unity
c* 0 => general case
c* non-zero => unit diagonal
c* itrans - indicates whether to use transpose
c* 0 => no
c* non-zero => use transpose
c* idim - leading dimension of a
c* n - number of unknowns
c* maxnz - maximum number of non-zero elements in any row of a
c* a - matrix
c* ja - column map for a
c*
c* In/Out:
c* x - source vector on input, solution 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_SolveUTriang_ELL (idiag, itrans, idim, n,
& maxnz, a, ja, x, status)
implicit none
c
c ... Input:
integer idiag, itrans, idim, maxnz, n
integer ja(idim,maxnz)
real a(idim,maxnz)
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j, ji
real zero
c
parameter (zero=0.0d0)
c
c ... Set return status.
status = 0
c
c ... Check arguments.
if (idim.le.0 .or. n.le.0 .or. maxnz.le.0 .or. n.gt.idim) then
status = -1
return
endif
c
c ... Check diagonal elements.
if (idiag .eq. 0) then
do i=1,n
if (a(i,1) .eq. zero) then
status = -4
return
endif
enddo
endif
c
c ... Back-substitution.
if (itrans .eq. 0) then
c
c .... Use coeff. as is.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal (used in ILU preconditioning).
x(n) = x(n)/a(n,1)
do i=n-1,1,-1
do j=2,maxnz
if (ja(i,j) .eq. 0) then
goto 10
elseif (ja(i,j) .gt. i) then
x(i) = x(i) - a(i,j)*x(ja(i,j))
endif
enddo
10 continue
x(i) = x(i)/a(i,1)
enddo
else
c
c ..... Unit diagonal (UNUSED, UNTESTED).
do i=n-1,1,-1
do j=2,maxnz
if (ja(i,j) .eq. 0) then
goto 15
elseif (ja(i,j) .gt. i) then
x(i) = x(i) - a(i,j)*x(ja(i,j))
endif
enddo
15 continue
enddo
endif
else
c
c .... Use transpose.
if (idiag .eq. 0) then
c
c ..... Non-unit diagonal (used in IC preconditioning).
x(n) = x(n)/a(n,1)
do i=n-1,1,-1
do j=i+1,n
c
c ....... Find ji such that ja(j,ji) = i. If one is found, we also
c know that a(j,ji) is non-zero.
do ji=2,maxnz
if (ja(j,ji) .eq. 0) then
c
c ......... Out of elements, so a(j,ji) = 0, so bail out and get
c another j.
goto 20
elseif (ja(j,ji) .eq. i) then
c
c ......... Found ji.
goto 19
endif
enddo
c
c ....... Didn't find ji.
goto 20
c
c ....... Modify i-th element of x.
19 continue
x(i) = x(i) - a(j,ji)*x(j)
c
c ....... Get another j.
20 continue
enddo
c
c ...... Finished with i-th element of x.
x(i) = x(i) / a(i,1)
enddo
else
c
c ..... Unit diagonal (UNUSED, UNTESTED).
do i=n-1,1,-1
do j=i+1,n
c
c ....... Find ji such that ja(j,ji) = i. If one is found, we also
c know that a(j,ji) is non-zero.
do ji=2,maxnz
if (ja(j,ji) .eq. 0) then
c
c ......... Out of elements, so a(j,ji) = 0, so bail out and get
c another j.
goto 25
elseif (ja(j,ji) .eq. i) then
c
c ......... Found ji.
goto 24
endif
enddo
c
c ....... Didn't find ji.
goto 25
c
c ....... Modify i-th element of x.
24 continue
x(i) = x(i) - a(j,ji)*x(j)
25 continue
enddo
enddo
endif
endif
c
return
end