CompressSyst.F

#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c*    $Id: CompressSyst.F,v 1.4 1996/04/24 19:14:00 turner Exp $
c*
c*    Removes redundant unknowns in a linear system.
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,ia(_JT_nrows_)
c*      x_large(i) = x_small(ix(i))
c*    enddo
c*
c*    <PARAMETER LIST>
c*
c*     In/Out:
c*      a - matrix
c*      ia - integer vector containing info about how "a" is stored
c*        NOTE: see description of ia below
c*      ja - column map for matrix
c*      x - original vector of unknowns
c*      b - source vector
c*
c*     Output:
c*      ix - scatter vector
c*      status - return status
c*        -1  ==>  invalid argument(s)
c*         0  ==>  success
c*
#include "iadesc.inc"
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     JT_CompressSystem_ELL   for A stored in ELLPACK-ITPACK format
c*     JT_CompressSystem_Full  for A stored in full conventional format
c*
#include "copyright.inc"
c**********************************************************************
      subroutine JT_CompressSystem (a, ia, ja, x, b, ix, status)
      implicit none
c
c ... In/Out:
      integer ia(_JT_no_of_storage_parameters_), ja(*)
      real a(*)
      real x(*)
c
c ... Output:
      integer status
      integer ix(*)
      real b(*)
c
      if (ia(_JT_storage_) .eq. _JT_storage_full_) then     ! Full conventional format.
       call JT_CompressSystem_Full (ia(_JT_idim_), ia(_JT_nrows_),
     &      a, x, b, ix, status)
      elseif (ia(_JT_storage_) .eq. _JT_storage_ELL_) then  ! ELLPACK-ITPACK format.
       call JT_CompressSystem_ELL (ia(_JT_idim_), ia(_JT_maxnz_),
     &      ia(_JT_nrows_), a, ja, x, b, ix, status)
      else
       status = -1
      endif
c
      return
      end