ApplyPreconditioner.F

#include "parmdefines.h"
#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c*    $Id: ApplyPreconditioner.F,v 1.30 1996/04/24 19:34:29 turner Exp $
c*
c*    Applies preconditioner as appropriate.
c*
c*    <PARAMETER LIST>
c*
c*     Input:
c*      iparm - array of integer parameters
c*      rparm - array of floating point parameters
c*      a - coefficient matrix
c*      ia - integer vector containing info about how "a" is stored
c*        NOTE: see description of ia below
c*      ja - integer array that maps columns in a to true columns
c*      ap - preconditioning matrix
c*      iap - integer vector containing info about how "ap" is stored
c*        NOTE: see description of ia below
c*      jap - column map for preconditioning matrix
c*      x - x-vector (source)
c*
c*     Output:
c*      y - y-vector (solution)
c*      status - return status
c*        -3  ==>  internal error
c*        -2  ==>  memory allocation failure
c*        -1  ==>  invalid argument(s)
c*         0  ==>  success
c*
#include "parmdesc.inc"
c*
#include "iadesc.inc"
c*
c*     Variables in iap are analogous to those in ia.
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     JT_FillVectorFloat
c*     JT_FlushUnit
c*     JT_Jacobi
c*     JT_SolveLTriang
c*     JT_SolveUTriang
c*     JT_SOR
c*     JT_WriteVectorFloat
c*     JT_y_eq_Ax
c*     JT_y_eq_diagAx
c*     JT_y_eq_invdiagAx
c*     JT_y_eq_x
c*
c*    <UNDOCUMENTED FEATURE>
c*
c*     Two simple implementations of Jacobi preconditioning have been
c*     implemented, primarily as templates for development of other
c*     preconditioners.
c*
c*     These preconditioners can be accessed by setting iparm(_JT_pre_) to -2 or 
c*     -3:
c*          -2 ==> Jacobi approximation of the coefficient
c*          -3 ==> Jacobi approximation of the inverse of coefficient
c*
#include "copyright.inc"
c**********************************************************************
      subroutine JT_ApplyPreconditioner (iparm, rparm,
     &           a, ia, ja, ap, iap, jap, x, y, status)
      implicit none
c
c ... Input:
      integer ia(_JT_no_of_storage_parameters_), ja(*)
      integer iap(_JT_no_of_storage_parameters_), jap(*)
      integer iparm(_JT_no_of_iparms_)
      real rparm(_JT_no_of_rparms_)
      real a(*), ap(*)
#ifdef strict_f77
      real x(*)
#else
      real x(ia(_JT_nrows_))
#endif
c
c ... Output:
      integer status
#ifdef strict_f77
      real y(*)
#else
      real y(ia(_JT_nrows_))
#endif
c
c ... Local:
      integer i
      real zero
c
      parameter (zero=0.0d0)
c
      if (iparm(_JT_pre_) .eq. _JT_pre_none_) then
c
c .... No preconditioning.
       call JT_y_eq_x (ia(_JT_nrows_), x, y, status)
      elseif (iparm(_JT_pre_) .eq. _JT_pre_Jacobi_) then
c
c .... m-step Jacobi preconditioning with relaxation.
       call JT_FillVectorFloat (ia(_JT_nrows_), zero, y, status)
       do i=1,iparm(_JT_steps_)
        call JT_Jacobi (rparm(_JT_omega_), x, a, ia, ja, y, status)
        if (status .lt. 0) goto 9999
       enddo
      elseif (iparm(_JT_pre_) .eq. _JT_pre_SSOR_) then
c
c .... m-step SSOR preconditioning.
       call JT_FillVectorFloat (ia(_JT_nrows_), zero, y, status)
       do i=1,iparm(_JT_steps_)
        call JT_SOR (1, rparm(_JT_omega_), x, a, ia, ja, y, status)
        if (status .lt. 0) goto 9999
       enddo
      elseif (iparm(_JT_pre_) .eq. _JT_pre_IC_) then
c
c .... Ap is an incomplete factorization, stored as a lower
c      triangular matrix.
       call JT_y_eq_x (ia(_JT_nrows_), x, y, status)
       call JT_SolveLTriang (0, 0, ap, iap, jap, y, status)
       if (status .lt. 0) goto 9999
       if (iparm(_JT_out_) .ge. 6) then
        call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), y, 
     &      'JT_ApplyPreconditioner: AFTER FORWARD SWEEP (y):', status)
#ifndef dont_flush
        call JT_FlushUnit (iparm(_JT_luout_), status)
#endif
       endif
       call JT_SolveUTriang (0, 1, ap, iap, jap, y, status)
      elseif (iparm(_JT_pre_) .eq. _JT_pre_ILU_) then
c
c .... Ap is an incomplete factorization, stored as a unit
c      lower triangular matrix and an upper triangular matrix.
       call JT_y_eq_x (ia(_JT_nrows_), x, y, status)
       call JT_SolveLTriang (1, 0, ap, iap, jap, y, status)
       if (status .lt. 0) goto 9999
       if (iparm(_JT_out_) .ge. 6) then
        call JT_WriteVectorFloat (iparm(_JT_luout_), ia(_JT_nrows_), y, 
     &      'JT_ApplyPreconditioner: AFTER FORWARD SWEEP (y):', status)
#ifndef dont_flush
        call JT_FlushUnit (iparm(_JT_luout_), status)
#endif
       endif
       call JT_SolveUTriang (0, 0, ap, iap, jap, y, status)
      elseif (iparm(_JT_pre_) .eq. -2) then
c
c .... Array ap contains diag(a).
       call JT_y_eq_invdiagAx (x, ap, iap, jap, y, status)
      elseif (iparm(_JT_pre_) .eq. -3) then
c
c .... Array ap contains inv(diag(a)).
       call JT_y_eq_diagAx (x, ap, iap, jap, y, status)
      else
c
c .... Invalid value of iparm(_JT_pre_).
       status = -1
       return
      endif
c
 9999 continue
c
c ... Check for internal error.
      if (status .lt. 0) then
       if (status .ne. -2) status = -3
      endif
c
      return
      end