y_minus_Ax.F

#include "iadefines.h"
c**********************************************************************
#include "author.inc"
c*    $Id: y_minus_Ax.F,v 1.20 1996/04/24 19:15:14 turner Exp $
c*
c*    Computes y=y-Ax.
c*
c*    <PARAMETER LIST>
c*
c*     Input:
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 - x-vector
c*
c*     In/Out:
c*      y - y-vector
c*
c*     Output:
c*      status - return status
c*        -3  ==>  internal error
c*        -2  ==>  memory allocation failure
c*        -1  ==>  invalid argument(s)
c*         0  ==>  success
c*
#include "iadesc.inc"
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     User_y_eq_Ax  user-provided routine to compute y=Ax (actual name
c*                   passed in parameter list)
c*     JT_y_eq_y_minus_Ax_Full  for A stored in full conventional format
c*     JT_y_eq_y_minus_Ax_ELL   for A stored in ELLPACK-ITPACK format
c*     JT_y_eq_y_minus_Ax_COO   for A stored in coordinate format
c*     JT_y_eq_y_minus_Ax_RSS   for A stored in row sequential storage
c*     JT_y_eq_y_minus_Ax_CSS   for A stored in column sequential storage
c*     JT_y_eq_y_minus_x
c*
#include "copyright.inc"
c**********************************************************************
      subroutine JT_y_eq_y_minus_Ax (a, ia, ja, x, y, status)
      implicit none
c
c ... Input:
      integer ia(_JT_no_of_storage_parameters_), ja(*)
      real a(*)
#ifdef strict_f77
      real x(*)
#else
      real x(ia(_JT_nrows_))
#endif
c
c ... In/Out:
#ifdef strict_f77
      real y(*)
#else
      real y(ia(_JT_nrows_))
#endif
c
c ... Output:
      integer status
c
      if (ia(_JT_storage_) .eq. _JT_storage_full_) then
       call JT_y_eq_y_minus_Ax_Full (ia(_JT_idim_), ia(_JT_nrows_), a, x, y, status)
      elseif (ia(_JT_storage_) .eq. _JT_storage_ELL_) then
       call JT_y_eq_y_minus_Ax_ELL
     &      (ia(_JT_idim_), ia(_JT_nrows_), ia(_JT_maxnz_), a, ja, x, y, status)
       if (status .eq. -1) status = -3
      elseif (ia(_JT_storage_) .eq. _JT_storage_COO_) then
       call JT_y_eq_y_minus_Ax_COO
     &      (ia(_JT_nelem_), ia(_JT_nrows_), a, ja, x, y, status)
       if (status .eq. -1) status = -3
      elseif (ia(_JT_storage_) .eq. _JT_storage_RSS_) then
       call JT_y_eq_y_minus_Ax_RSS
     &      (ia(_JT_nelem_), ia(_JT_nrows_), a, ja, x, y, status)
       if (status .eq. -1) status = -3
      elseif (ia(_JT_storage_) .eq. _JT_storage_CSS_) then
       call JT_y_eq_y_minus_Ax_CSS
     &      (ia(_JT_nelem_), ia(_JT_nrows_), a, ja, x, y, status)
       if (status .eq. -1) status = -3
      else  ! invalid value for ia(_JT_storage_)
       status = -1
      endif
c
      return
      end