CompressSyst_spec.F
c**********************************************************************
#include "author.inc"
c* $Id: CompressSyst_spec.F,v 1.1 1995/11/14 02:20:04 turner Exp $
c*
c* Removes redundant unknowns in a linear system stored in full
c* conventional format.
c*
c* For example, consider a simple 3x3 system:
c*
c* [ a11 a12 a13 ] [ x1 ] [ b1 ]
c* [ a21 a22 a23 ] [ x2 ] = [ b2 ]
c* [ a31 a32 a33 ] [ x3 ] [ b3 ]
c*
c* Now suppose for some reason you wanted to force x2 to have the
c* same result as x3. That implies:
c*
c* [ a11 a12 a13 ] [ x1 ] [ b1 ]
c* [ 0 1 -1 ] [ x2 ] = [ 0 ]
c* [ a31 a32 a33 ] [ x3 ] [ b3 ]
c*
c* If the system is large, solving this system would be significantly
c* more costly than solving the equivalent smaller system. In this
c* simple case, the smaller system is:
c*
c* [ a11 a12+a13 ] [ x1 ] = [ b1 ]
c* [ a31 a32+a33 ] [ x3 ] [ b3 ]
c*
c* This routine replaces the original large system with the
c* equivalent small system.
c*
c* An integer vector ix mapping the unknowns from the small system
c* into the larger system is also returned. That is, after solving
c* the small system, the following loop would place them into a
c* vector the length of the original (large) system:
c*
c* do i=1,nrows_a
c* x_large(i) = x_small(ix(i))
c* enddo
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of a
c*
c* In/Out:
c* n - number of unknowns
c* a - matrix
c* x - unknown vector
c* b - source vector
c*
c* Output:
c* ix - scatter vector
c* status - return status
c*
#include "copyright.inc"
c**********************************************************************
subroutine JT_CompressSystem_Full (idim, n, a, x, b, ix, status)
implicit none
c
c ... Input:
integer idim
c
c ... In/Out:
integer n
real a(idim,n), x(n), b(n)
c
c ... Output:
integer status
integer ix(n)
c
c ... Local:
integer i, ii, j, m, retain
real zero, one
parameter (zero=0.0d0, one=1.0d0)
c
c ... Initialize scatter vector.
do i=1,n
ix(i) = i
enddo
c
c ... Initialize m. This will keep track of the size of the system
c as it is reduced. Note that most loops in the main loop of
c the code go to m rather than n.
m = n
c
c ... Loop through each row (equation).
do i=1,n
c
c .... Check to see if equation has form indicating it is redundant.
c
c First, the source must be zero.
if (b(i) .ne. zero) goto 10
c
c Next, it must have a unit diagonal.
if (a(i,i) .ne. one) goto 10
c
c Next, it must have one and only one coefficient equal to -1.
c If a -1 is found, the column containing it indicates the
c current index (prior to deleting the duplicate unknown) of the
c unknown to be retained.
retain = 0
do j=1,m
if (a(i,j) .eq. -one) then
if (retain .ne. 0) then
goto 10
else
retain = j
endif
endif
enddo
if (retain .eq. 0) goto 10
c
c Finally, it must have all other coefficients equal to zero.
do j=1,m
if (j.ne.i .and. j.ne.retain .and. a(i,j).ne.zero) goto 10
enddo
c
c .... If here, then ith unknown is redundant and must be eliminated.
c
c .... First modify the coefficients of the remaining rows.
do ii=1,m
if (ii .ne. i) then
a(ii,retain) = a(ii,retain) + a(ii,i)
endif
enddo
c
c .... Want to modify the element of the scatter vector that points to
c the unknown being deleted, but the correspondence between
c elements of ix and rows in the system is lost after any unknowns
c are deleted. So just make any element of ix that points to
c the unknown being deleted point to the unknown being retained.
do ii=1,n
if (ix(ii) .eq. i) ix(ii) = retain
enddo
c
c .... Delete ith row.
do j=1,m
do ii=i+1,m
a(ii-1,j) = a(ii,j)
enddo
enddo
do ii=i+1,m
x(ii-1) = x(ii)
b(ii-1) = b(ii)
enddo
c
c .... Also need to decrement elements of ix that point to unknowns with
c index greater than or equal to the one being deleted. This could
c possibly be combined with the modification of the scatter vector
c above.
do ii=1,n
if (ix(ii) .ge. i) ix(ii) = ix(ii) - 1
enddo
c
c .... Delete ith column.
do j=i+1,m
do ii=1,m-1
a(ii,j-1) = a(ii,j)
enddo
enddo
c
c .... Decrement m.
m = m - 1
c
10 continue
c
c .... Check to see if out of rows.
if (i .ge. m) goto 20
c
enddo
20 continue
c
c ... Set size of compressed system.
n = m
c
status = 0
return
end
c**********************************************************************
#include "author.inc"
c* $Id: CompressSyst_spec.F,v 1.1 1995/11/14 02:20:04 turner Exp $
c*
c* Removes redundant unknowns in a linear system stored in ELLPACK-
c* ITPACK format.
c*
c* For example, consider a simple 3x3 system in full format:
c*
c* [ a11 a12 a13 ] [ x1 ] [ b1 ]
c* [ a21 a22 a23 ] [ x2 ] = [ b2 ]
c* [ a31 a32 a33 ] [ x3 ] [ b3 ]
c*
c* Now suppose for some reason you wanted to force x2 to have the
c* same result as x3. That implies:
c*
c* [ a11 a12 a13 ] [ x1 ] [ b1 ]
c* [ 0 1 -1 ] [ x2 ] = [ 0 ]
c* [ a31 a32 a33 ] [ x3 ] [ b3 ]
c*
c* If the system is large, solving this system would be significantly
c* more costly than solving the equivalent smaller system. In this
c* simple case, the smaller system is:
c*
c* [ a11 a12+a13 ] [ x1 ] = [ b1 ]
c* [ a31 a32+a33 ] [ x3 ] [ b3 ]
c*
c* This routine replaces the original large system with the
c* equivalent small system.
c*
c* NOTE: This routine is not completely general. Specifically,
c* it is possible that ja would need to be expanded beyond the
c* initial number of columns. This routine
c*
c* An integer vector ix mapping the unknowns from the small system
c* into the larger system is also returned. That is, after solving
c* the small system, the following loop would place them into a
c* vector the length of the original (large) system:
c*
c* do i=1,nrows_a
c* x_large(i) = x_small(ix(i))
c* enddo
c*
c* <PARAMETER LIST>
c*
c* Input:
c* idim - leading dimension of a
c* maxnz - maximum number of non-zero elements in any row of a
c*
c* In/Out:
c* n - number of unknowns
c* a - matrix
c* x - unknown vector
c* b - source vector
c*
c* Output:
c* ix - scatter vector
c* status - return status
c*
#include "copyright.inc"
c**********************************************************************
subroutine JT_CompressSystem_ELL (idim, maxnz, n, a, ja, x, b,
& ix, status)
implicit none
c
c ... Input:
integer idim, maxnz
c
c ... In/Out:
integer n
integer ja(idim,maxnz)
real a(idim,maxnz), x(n), b(n)
c
c ... Output:
integer status
integer ix(n)
c
c ... Local:
integer i, ii, j, keep, delete, m, retain
real zero, one
parameter (zero=0.0d0, one=1.0d0)
c
c ... Initialize scatter vector.
do i=1,n
ix(i) = i
enddo
c
c ... Initialize m. This will keep track of the size of the system
c as it is reduced. Note that most loops in the main loop of
c the code go to m rather than n.
m = n
c
c ... Loop through each row (equation).
do i=1,n
c
c .... Check to see if equation has form indicating it is redundant.
c
c First, the source must be zero.
if (b(i) .ne. zero) goto 10
c
c Next, it must have a unit diagonal.
if (a(i,1) .ne. one) goto 10
c
c Next, it must have one and only one coefficient equal to -1.
c If a -1 is found, the column containing it indicates the
c current index (prior to deleting the duplicate unknown) of
c the unknown to be retained.
retain = 0
do j=1,maxnz
if (a(i,j) .eq. -one) then
if (retain .ne. 0) then
goto 10
else
retain = ja(i,j)
endif
endif
enddo
if (retain .eq. 0) goto 10
c
c Finally, it must have all other coefficients equal to zero.
do j=1,maxnz
if (ja(i,j).ne.i .and. ja(i,j).ne.retain .and.
& a(i,j).ne.zero) goto 10
enddo
c
c .... If here, then ith unknown is redundant and must be eliminated.
c
c .... First modify the coefficients of the remaining rows.
do ii=1,m
if (ii .ne. i) then
keep = 0
delete = 0
do j=1,maxnz
if (ja(ii,j) .eq. retain) then
keep = j
elseif (ja(ii,j) .eq. i) then
delete = j
endif
enddo
c
if (delete .ne. 0) then
c
c Index of unknown to be deleted was found, which means that the
c coefficient for the one to be kept must be modified.
if (keep .ne. 0) then
c
c Index of unknown to be kept was found, so modify that element
c and remove the element that was multiplying the unknown being
c deleted.
a(ii,keep) = a(ii,keep) + a(ii,delete)
a(ii,delete) = zero
ja(ii,delete) = 0
else
c
c Index of unknown to be kept was not found, so just change the
c column map entry for the unknown to be deleted to point to
c the retained unknown. That element of the coeff. already
c contains the correct value.
ja(ii,delete) = retain
endif
endif
endif
enddo
c
c .... Want to modify the element of the scatter vector that points to
c the unknown being deleted, but the correspondence between
c elements of ix and rows in the system is lost after any unknowns
c are deleted. So just make any element of ix that points to
c the unknown being deleted point to the unknown being retained.
do ii=1,n
if (ix(ii) .eq. i) ix(ii) = retain
enddo
c
c .... Delete ith row.
do j=1,maxnz
do ii=i+1,m
a(ii-1,j) = a(ii,j)
ja(ii-1,j) = ja(ii,j)
enddo
enddo
do ii=i+1,m
x(ii-1) = x(ii)
b(ii-1) = b(ii)
enddo
c
c .... Also need to decrement elements of ix that point to unknowns with
c index greater than or equal to the one being deleted. This could
c possibly be combined with the modification of the scatter vector
c above.
do ii=1,n
if (ix(ii) .ge. i) ix(ii) = ix(ii) - 1
enddo
c
c .... Delete ith column.
do j=1,maxnz
do ii=1,m-1
if (ja(ii,j) .gt. i) then
ja(ii,j) = ja(ii,j) - 1
endif
enddo
enddo
c
c .... Decrement m.
m = m - 1
c
10 continue
c
c .... Check to see if out of rows.
if (i .ge. m) goto 20
c
enddo
20 continue
c
c ... Set size of compressed system.
n = m
c
status = 0
return
end