GeneratePreconditioner.F

#include "parmdefines.h"
#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c*    $Id: GeneratePreconditioner.F,v 1.26 1996/04/24 19:34:33 turner Exp $
c*
c*    Computes an approximation to a or an approximation to inv(a)
c*    as appropriate, to be used as a preconditioner.
c*
c*    <PARAMETER LIST>
c*
c*     Input:
c*      iparm - array of integer parameters
c*      rparm - array of floating point parameters
c*      a - coeff. matrix
c*      ia - integer vector containing info about how "a" is stored
c*        NOTE: see description of ia below
c*      ja - column index array for a
c*
c*     Output:
c*      ap - preconditioner
c*      iap - integer vector containing info about how "ap" is stored
c*        NOTE: see description of ia below
c*      jap - column index array for ap
c*      status - return status
c*        -3  ==>  internal error
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_B_eq_A
c*     JT_IC
c*     JT_ILU
c*     JT_diagB_eq_diagA
c*     JT_diagB_eq_invdiagA
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_GeneratePreconditioner (iparm, rparm,
     &           a, ia, ja, ap, iap, jap, status)
      implicit none
c
c ... Input:
      integer ia(_JT_no_of_storage_parameters_), ja(*)
      integer iparm(_JT_no_of_iparms_)
      real rparm(_JT_no_of_rparms_)
      real a(*)
c
c ... Output:
      integer status
      integer iap(_JT_no_of_storage_parameters_), jap(*)
      real ap(*)
c
c ... Check for invalid value of iparm(_JT_pre_).
      if (iparm(_JT_pre_).lt.-3 .or. iparm(_JT_pre_).gt.4) then
       status = -1
       return
      endif
c
c ... Construct preconditioner as appropriate.
      if (iparm(_JT_pre_) .eq. _JT_pre_IC_) then
c
c .... Incomplete Cholesky factorization.
       call JT_B_eq_A (a, ia, ja, ap, iap, jap, status)
       call JT_IC (rparm(_JT_omega_), ap, iap, jap, status)
      elseif (iparm(_JT_pre_) .eq. _JT_pre_ILU_) then
c
c .... Incomplete LU factorization (store unit lower triangular matrix
c      and upper triangular matrix in ap).
       call JT_B_eq_A (a, ia, ja, ap, iap, jap, status)
       call JT_ILU (rparm(_JT_omega_), ap, iap, jap, status)
      elseif (iparm(_JT_pre_) .eq. -2) then
c
c .... Approximate a by diag(a).
       call JT_diagB_eq_diagA (a, ia, ja, ap, iap, jap, status)
      elseif (iparm(_JT_pre_) .eq. -3) then
c
c .... Approximate the inverse of a by inv(diag(a)).
       call JT_diagB_eq_invdiagA (a, ia, ja, ap, iap, jap, status)
      endif
      if (status .lt. 0) goto 9999
c
      if (iparm(_JT_out_) .ge. 4) then
       call JT_WriteMatrix (iparm(_JT_luout_), ap, iap, jap,
     &     'JT_GeneratePreconditioner: PRECONDITIONER:', status)
      endif
c
 9999 continue
c
c ... Check for internal error.
      if (status .lt. 0) status = -3
c
      return
      end