SOR_spec.F
c***********************************************************************
#include "author.inc"
c* $Id: SOR_spec.F,v 1.1 1995/11/14 02:19:59 turner Exp $
c*
c* Performs one iteration of SOR or SSOR for a system in full-storage
c* format.
c*
c* NOTES:
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* imeth - determines whether or not to perform SSOR
c* idim - leading dimension of array a
c* n - number of unknowns
c* omega - relaxation parameter
c* b - source vector
c* a - coefficient matrix
c*
c* In/Out:
c* x - solution vector
c*
c* Output:
c* status - return status
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_SOR_Full (imeth, idim, n, omega, b, a, x, status)
implicit none
c
c ... Input:
integer imeth, idim, n
real a(idim,n), b(n), omega
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j
real sum, zero, two
c
parameter (zero=0.0d0, two=2.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (idim.le.0 .or. n.le.0 .or. n.gt.idim .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
do i=1,n
sum = zero
do j = 1,i-1
sum = sum + a(i,j)*x(j)
enddo
do j = i+1,n
sum = sum + a(i,j)*x(j)
enddo
sum = ( b(i) - sum ) / a(i,i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
if (imeth .eq. 1) then
do i=n,1,-1
sum = zero
do j = 1,i-1
sum = sum + a(i,j)*x(j)
enddo
do j = i+1,n
sum = sum + a(i,j)*x(j)
enddo
sum = ( b(i) - sum ) / a(i,i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
endif
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: SOR_spec.F,v 1.1 1995/11/14 02:19:59 turner Exp $
c*
c* Performs one iteration of SOR or SSOR for a system in ELL format.
c*
c* NOTES:
c* - The main diagonal of the coefficient is assumed to be in
c* the first column.
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* imeth - determines whether or not to perform SSOR
c* idim - leading dimension of arrays a, ja
c* n - number of unknowns
c* maxnz - maximum number of non-zero elements in any row
c* omega - relaxation parameter
c* b - source vector
c* a - coefficient matrix
c* ja - integer map for coefficient
c*
c* In/Out:
c* x - solution vector
c*
c* Output:
c* status - return status
c* -2 ==> memory allocation failure
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
c* <SUBROUTINES REQUIRED>
c*
c* JT_y_eq_x
c*
#include "arrays-SOR_ELL.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_SOR_ELL (imeth, idim, n, maxnz, omega, b, a, ja,
& x, status)
implicit none
c
c ... Input:
integer imeth, idim, maxnz, n
integer ja(idim,maxnz)
real a(idim,maxnz), b(n), omega
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j
real sum, zero, two
#include "declare-SOR_ELL.inc"
c
parameter (zero=0.0d0, two=2.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (idim.le.0 .or. n.le.0 .or. maxnz.le.0 .or. n.gt.idim .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
#include "allocate-SOR_ELL.inc"
c
c ... Copy x into a temporary array that starts at zero and use that
c array in the iteration (since ja(i,j) can be zero).
call JT_y_eq_x (n, x, xtmp(1), status)
xtmp(0) = zero
c
do i=1,n
if (maxnz .eq. 5) then
sum = a(i,2)*xtmp(ja(i,2)) +
& a(i,3)*xtmp(ja(i,3)) +
& a(i,4)*xtmp(ja(i,4)) +
& a(i,5)*xtmp(ja(i,5))
elseif (maxnz .eq. 7) then
sum = a(i,2)*xtmp(ja(i,2)) +
& a(i,3)*xtmp(ja(i,3)) +
& a(i,4)*xtmp(ja(i,4)) +
& a(i,5)*xtmp(ja(i,5)) +
& a(i,6)*xtmp(ja(i,6)) +
& a(i,7)*xtmp(ja(i,7))
else
sum = zero
do j = 2,maxnz
sum = sum + a(i,j)*xtmp(ja(i,j))
enddo
endif
sum = ( b(i) - sum ) / a(i,1)
xtmp(i) = xtmp(i) + omega*( sum - xtmp(i) )
enddo
c
if (imeth .eq. 1) then
do i=n,1,-1
if (maxnz .eq. 5) then
sum = a(i,2)*xtmp(ja(i,2)) +
& a(i,3)*xtmp(ja(i,3)) +
& a(i,4)*xtmp(ja(i,4)) +
& a(i,5)*xtmp(ja(i,5))
elseif (maxnz .eq. 7) then
sum = a(i,2)*xtmp(ja(i,2)) +
& a(i,3)*xtmp(ja(i,3)) +
& a(i,4)*xtmp(ja(i,4)) +
& a(i,5)*xtmp(ja(i,5)) +
& a(i,6)*xtmp(ja(i,6)) +
& a(i,7)*xtmp(ja(i,7))
else
sum = zero
do j = 2,maxnz
sum = sum + a(i,j)*xtmp(ja(i,j))
enddo
endif
sum = ( b(i) - sum ) / a(i,1)
xtmp(i) = xtmp(i) + omega*( sum - xtmp(i) )
enddo
endif
c
c ... Copy new iterate back into x.
call JT_y_eq_x (n, xtmp(1), x, status)
c
9999 continue
#include "deallocate-SOR_ELL.inc"
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: SOR_spec.F,v 1.1 1995/11/14 02:19:59 turner Exp $
c*
c* Performs one iteration of SOR or SSOR for a system in coordinate
c* format.
c*
c* NOTES:
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* imeth - determines whether or not to perform SSOR
c* nelem - number of active elements in array a
c* n - number of unknowns
c* omega - relaxation parameter
c* b - source vector
c* a - coefficient matrix
c* ja - integer map for coefficient
c*
c* In/Out:
c* x - solution vector
c*
c* Output:
c* status - return status
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_SOR_COO (imeth, nelem, n, omega, b, a, ja, x,
& status)
implicit none
c
c ... Input:
integer imeth, nelem, n
integer ja(2*nelem)
real a(nelem), b(n), omega
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j
real sum, zero, two
c
parameter (zero=0.0d0, two=2.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (nelem.le.0 .or. n.le.0 .or. n.gt.nelem .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
do i=1,n
sum = zero
do j=n+1,nelem
if (ja(j) .eq. i) then
sum = sum + a(j)*x(ja(j+nelem))
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
if (imeth .eq. 1) then
do i=n,1,-1
sum = zero
do j=n+1,nelem
if (ja(j) .eq. i) then
sum = sum + a(j)*x(ja(j+nelem))
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
endif
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: SOR_spec.F,v 1.1 1995/11/14 02:19:59 turner Exp $
c*
c* Performs one iteration of SOR or SSOR for a system in RSS format.
c*
c* NOTES:
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* imeth - determines whether or not to perform SSOR
c* nelem - number of active elements in array a
c* n - number of unknowns
c* omega - relaxation parameter
c* b - source vector
c* a - coefficient matrix
c* ja - integer map for coefficient
c*
c* In/Out:
c* x - solution vector
c*
c* Output:
c* status - return status
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_SOR_RSS (imeth, nelem, n, omega, b, a, ja, x,
& status)
implicit none
c
c ... Input:
integer imeth, nelem, n
integer ja(nelem)
real a(nelem), b(n), omega
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j, ii, jj
real sum, zero, two
c
parameter (zero=0.0d0, two=2.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (nelem.le.0 .or. n.le.0 .or. n.gt.nelem .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
do i=1,n
sum = zero
do j=n+1,nelem
ii = (ja(j) - 1)/n + 1
if (ii .eq. i) then
jj = ja(j) - (ii-1)*n
sum = sum + a(j)*x(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
if (imeth .eq. 1) then
do i=n,1,-1
sum = zero
do j=n+1,nelem
ii = (ja(j) - 1)/n + 1
if (ii .eq. i) then
jj = ja(j) - (ii-1)*n
sum = sum + a(j)*x(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
endif
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: SOR_spec.F,v 1.1 1995/11/14 02:19:59 turner Exp $
c*
c* Performs one iteration of SOR or SSOR for a system in CSS format.
c*
c* NOTES:
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
c* imeth - determines whether or not to perform SSOR
c* nelem - number of active elements in array a
c* n - number of unknowns
c* omega - relaxation parameter
c* b - source vector
c* a - coefficient matrix
c* ja - integer map for coefficient
c*
c* In/Out:
c* x - solution vector
c*
c* Output:
c* status - return status
c* -1 ==> invalid argument(s)
c* 0 ==> success
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_SOR_CSS (imeth, nelem, n, omega, b, a, ja, x,
& status)
implicit none
c
c ... Input:
integer imeth, nelem, n
integer ja(nelem)
real a(nelem), b(n), omega
c
c ... In/Out:
real x(n)
c
c ... Output:
integer status
c
c ... Local:
integer i, j, ii, jj
real sum, zero, two
c
parameter (zero=0.0d0, two=2.0d0)
c
c ... Initialize return status.
status = 0
c
c ... Check arguments.
if (nelem.le.0 .or. n.le.0 .or. n.gt.nelem .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
do i=1,n
sum = zero
do j=n+1,nelem
jj = (ja(j) - 1)/n + 1
ii = ja(j) - (jj-1)*n
if (ii .eq. i) then
sum = sum + a(j)*x(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
if (imeth .eq. 1) then
do i=n,1,-1
sum = zero
do j=n+1,nelem
jj = (ja(j) - 1)/n + 1
ii = ja(j) - (jj-1)*n
if (ii .eq. i) then
sum = sum + a(j)*x(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
endif
c
return
end