3. Object-Based Design

Although even JTPACK77 implements a form of object-based design, it is necessarily crude due to the lack of syntactic support for such programming in F77 (for more on how this is accomplished in JTPACK77, see [10]). F90 provides a number of syntactic advances over F77, including:

• derived types, similar to C structures, which provide user-defined types which may consist of entities of various types,
• generic procedures, which provide polymorphism,
• array-valued functions,
• overloading of intrinsics, and
• modules, which allow encapsulation and access control via public and private attributes.

Modules in particular are a powerful addition to the language, and although they can be used in a number of ways, they are often used to group entities that are related in some manner. Two extreme options for a package such as JTPACK90 would be:

• Bundle routines functionally. That is, all routines that perform matrix-vector multiplication for the various sparse storage types could be grouped in a single module. That module would then be used by any routine needing to perform matrix-vector multiplication.
• Use modules to create ``classes''. That is, bundle a type definition, along with all the routines necessary to perform operations using that type.

In JTPACK90 we have chosen primarily the latter approach. Although all routines in JTPACK90 are encapsulated in modules, there are two primary types of modules.

• ``Class'' modules each contain a derived-type definition for a sparse storage type, along with all the routines to operate on that type. An edited version of the JTPACK90 module that defines the class for the ELLPACK-ITPACK (ELL) storage format [9] is shown below.

  module JT_ELL_module
    implicit none
    type JT_ELL_matrix
       real, dimension(:,:), pointer :: values
       integer, dimension(:,:), pointer :: map
    end type JT_ELL_matrix
    interface MatMul
       module procedure Ax
    end interface
    private
    public :: JT_ELL_matrix, MatMul
  contains
    function Ax(a,x)
      type(JT_ELL_matrix), intent(in) :: a
      real, intent(in), dimension(:) :: x
      real, dimension(SIZE(x)) :: Ax
      integer :: j
  
      Ax = zero
      do j=1,SIZE(a%values, dim=2)
         where (a%map(:,j) /= 0) Ax = Ax + a%values(:,j)*x(a%map(:,j))
      end do
  
      return
    end function Ax
  end module JT_ELL_module
  

  The definition of the JT_ELL_matrix type appears in the specification portion of the module, and consists of two rank-2 arrays, one real (for the values of the matrix) and one integer (for the column indices).

  Only the function for performing matrix-vector multiplication using ELL storage is shown. The real module also contains routines for assigning, loading, dumping, writing, extracting the diagonal from, computing the norm of, computing an incomplete factorization of, etc., an ELL object.

  Note that the F90 intrinsic for performing matrix multiplication, MatMul, is overloaded by defining a new routine to be used when the arguments match those of the module procedure Ax. Note also that Ax is an array-valued function, meaning that it returns an array rather than a scalar.

• ``Solver'' modules each contain all the routines that implement a particular algorithm, such as CG, GMRES, etc. A modified version of the JTPACK90 GMRES module is shown below.

  module JT_GMRES_module
    implicit none
    interface JT_GMRES
       module procedure GMRES_Full
       module procedure GMRES_ELL
    end interface
    private
    public :: JT_GMRES
  contains
    subroutine GMRES_Full (status, b, control, x, cpu, &
                           rnormt, errt, rnorm, err, a, ap)
      use JT_Full_module
      real, intent(in), dimension(:,:) :: a
      real, intent(inout), dimension(:,:) :: ap
  #include "GMRES-guts.F90"
    end subroutine GMRES_Full
    subroutine GMRES_ELL (status, b, control, x, cpu, &
                          rnormt, errt, rnorm, err, a, ap)
      use JT_ELL_module
      type(JT_ELL_matrix), intent(in) :: a
      type(JT_ELL_matrix), intent(inout) :: ap
  #include "GMRES-guts.F90"
    end subroutine GMRES_ELL
  end module JT_GMRES_module
  

  Two routines are shown, one for standard, full-storage coefficients and one for ELL coefficients. Note that the only difference between the two routines is the use statement and the declarations of the arrays. Everything else is common between the routines, and is hence pulled out and stored in a common file.

  Though in an important sense this shows the significant advances in abstraction allowed by F90, it also illustrates a shortfall, for these cpp acrobatics would not be necessary if F90 had something akin to templates in C++.

Note that this is a hybrid approach in that some of the routines dealing with matrices of a particular storage type reside in the solver modules rather than in the class modules. This was a conscious decision, since it makes adding a new solver much easier at the expense of requiring slightly more effort when adding a new storage type. With our approach, to add a new solver one must simply create a new solver module. Adding a new storage type requires creating a new class module as well as adding a new routine (though really just a ``template'' with the correct declarations and an include statement) to each of the solver modules.