Jacobi_spec.F
c***********************************************************************
#include "author.inc"
c* $Id: Jacobi_spec.F,v 1.1 1995/11/14 02:19:52 turner Exp $
c*
c* Performs one Jacobi iteration for a system in full-storage format.
c*
c* NOTES:
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
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* -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-Jacobi_Full.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Jacobi_Full (idim, n, omega, b, a, x, status)
implicit none
c
c ... Input:
integer 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
#include "declare-Jacobi_Full.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. n.gt.idim .or.
& omega.le.zero .or. omega.ge.two) then
status = -1
return
endif
c
#include "allocate-Jacobi_Full.inc"
c
call JT_y_eq_x (n, x, xtmp, status)
do i=1,n
sum = zero
do j = 1,n
sum = sum + a(i,j)*xtmp(j)
enddo
sum = sum - a(i,i)*xtmp(i)
sum = ( b(i) - sum ) / a(i,i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
9999 continue
#include "deallocate-Jacobi_Full.inc"
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: Jacobi_spec.F,v 1.1 1995/11/14 02:19:52 turner Exp $
c*
c* Performs one Jacobi iteration 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* 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-Jacobi_ELL.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Jacobi_ELL (idim, n, maxnz, omega, b, a, ja,
& x, status)
implicit none
c
c ... Input:
integer 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-Jacobi_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-Jacobi_ELL.inc"
c
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)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
9999 continue
#include "deallocate-Jacobi_ELL.inc"
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: Jacobi_spec.F,v 1.1 1995/11/14 02:19:52 turner Exp $
c*
c* Performs one Jacobi iteration for a system in coordinate format.
c*
c* NOTES:
c* - The main diagonal is assumed to be in the first n elements
c* of the coefficient.
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
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* -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-Jacobi_COO.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Jacobi_COO (nelem, n, omega, b, a, ja, x, status)
implicit none
c
c ... Input:
integer 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
#include "declare-Jacobi_COO.inc"
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
#include "allocate-Jacobi_COO.inc"
c
call JT_y_eq_x (n, x, xtmp, status)
do i=1,n
sum = zero
do j=n+1,nelem
if (ja(j) .eq. i) then
sum = sum + a(j)*xtmp(ja(j+nelem))
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
9999 continue
#include "deallocate-Jacobi_COO.inc"
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: Jacobi_spec.F,v 1.1 1995/11/14 02:19:52 turner Exp $
c*
c* Performs one Jacobi iteration for a system in RSS format.
c*
c* NOTES:
c* - The main diagonal is assumed to be in the first n elements
c* of the coefficient.
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
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* -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-Jacobi_RSS.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Jacobi_RSS (nelem, n, omega, b, a, ja, x, status)
implicit none
c
c ... Input:
integer 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
#include "declare-Jacobi_RSS.inc"
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
#include "allocate-Jacobi_RSS.inc"
c
call JT_y_eq_x (n, x, xtmp, status)
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)*xtmp(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
9999 continue
#include "deallocate-Jacobi_RSS.inc"
c
return
end
c***********************************************************************
#include "author.inc"
c* $Id: Jacobi_spec.F,v 1.1 1995/11/14 02:19:52 turner Exp $
c*
c* Performs one Jacobi iteration for a system in CSS format.
c*
c* NOTES:
c* - The main diagonal is assumed to be in the first n elements
c* of the coefficient.
c* - The main diagonal is assumed to be non-zero.
c*
c* <PARAMETER LIST>
c*
c* Input:
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* -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-Jacobi_CSS.inc"
c*
#include "copyright.inc"
c***********************************************************************
subroutine JT_Jacobi_CSS (nelem, n, omega, b, a, ja, x, status)
implicit none
c
c ... Input:
integer 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
#include "declare-Jacobi_CSS.inc"
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
#include "allocate-Jacobi_CSS.inc"
c
call JT_y_eq_x (n, x, xtmp, status)
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)*xtmp(jj)
endif
enddo
sum = ( b(i) - sum ) / a(i)
x(i) = x(i) + omega*( sum - x(i) )
enddo
c
9999 continue
#include "deallocate-Jacobi_CSS.inc"
c
return
end