psblas3/base/modules
Salvatore Filippone a695b63c3c New linmap internal structure.
..
aux Merge branch 'psblas-3.6-maint' of https://github.com/sfilippone/psblas3 into psblas-3.6-maint
comm New linmap internal structure.
desc New fnd_owner_halo methods. Bump version 3.6.1.
psblas New interface for methods with global reductions.
serial Implement clean_zeros in CSR and CSC
tools New RENUM_BLOCK function to create descriptor with contiguous global numbering.
Makefile Rename C implementation of hashval
README.F2003 Fixed README.
cutil.c psblas3:
error.f90 Updated copyright statement
fakempi.c psblas3:
psb_base_mod.f90 Updated copyright statement
psb_check_mod.f90 Updated copyright statement
psb_const_mod.F90 New fnd_owner_halo methods. Bump version 3.6.1.
psb_error_impl.F90 Simple caching for psb_info. Make barrier in errcomm optional, default false.
psb_error_mod.F90 Updated copyright statement
psb_internals.h psblas3:
psb_penv_mod.F90 psblas3:
psb_realloc_mod.F90 Updated copyright statement
psi_bcast_mod.F90 Updated copyright statement
psi_c_mod.f90 Updated copyright statement
psi_comm_buffers_mod.F90 Updated copyright statement
psi_d_mod.f90 Updated copyright statement
psi_i_mod.f90 Modified internals to have tmp_ovrlap in local numbering.
psi_mod.f90 Updated copyright statement
psi_p2p_mod.F90 Updated copyright statement
psi_penv_mod.F90 Take out HAVE_MOLD.
psi_reduce_mod.F90 Updated copyright statement
psi_s_mod.f90 Updated copyright statement
psi_z_mod.f90 Updated copyright statement

README.F2003

Design principles for this directory.

1. What is a sparse matrix? 
   It is an object which does have some properties (number of rows,
   number of columns, whether it is a triangle, and in that case
   upper/lower, unit/nonunit), a state (null, build, assembled,
   update), a type (real/complex, single/double), and a storage
   format.  
   Thus we have a three-level inheritance chain: 
   i.   The  base object, defining the methods to set/query the various
        properties, and allocate and  free. Some of the property
        getters/setters, allocate and free depend on the storage
        format, so at this level they will just throw an error. 
   ii.  The X_base_object, where X=s,d,c,z  thus defining the
        type. At this level we define the computational interfaces to
	MV and SV, since they require the type of the vectors/scalars
   	involved (should also add NRMI here!!!!), but again they will
   	be empty shells. We also define the interface to CSPUT,
   	required to build the object, and TO_COO,FROM_COO (see
   	below). 
   iii. The X_YYY_object where the real implementation of the
      	MV/SV/NRMI/CSPUT/ALLOCATE/FREE/TO_COO/FROM_COO takes place.  

2. What is a sparse matrix (take 2)? 
   The above structure by itself does not allow a sparse matrix to
   switch among different storage formats during its life. To do this,
   we define all of the above to be INNER objects, encapsulated in an
   OUTER object which is what the rest of the library sees, as
   follows: 
     
   type :: psbn_d_sparse_mat

     class(psbn_d_base_sparse_mat), allocatable  :: a 
    
   end type psbn_d_sparse_mat
   type(psbn_d_sparse_mat) :: a

   In this way we can have an outer object whose type is stable
   both statically (at compile time) and at runtime, while at runtime
   the type of the inner object switches from COO to CSR to whatever as
   needed. All of the methods are simply thrown onto the corresponding
   methods of the (allocatable, polymorphic) component A%A as needed
   (provided the component is allocated, that is).
   This is what is called a STATE design pattern (different from the
   internal state we discussed above).

   As an example, consider the allocate/build/assembly cycle: 
   the outer code would do the following: 
   1.  Allocate(psbn_d_coo_sparse_mat ::  a%a) 

   2. During the build loop a call to A%CSINS() gets translated into
      CALL A%A%CSINS()
    
   3. At assembly time the code would do the following
      subroutine psb_spasb(a,....)
      type(psbn_d_sparse_mat), intent(inout)  :: a 

      class(psbn_d_base_sparse_mat), allocatable :: temp

      select case (TYPE)
      case('CSR')
         allocate(psbn_d_csr_sparse_mat :: temp, stat=info)
      end select				 
      call temp%from_coo(a%a)
      call a%a%free()
      call move_alloc(temp,a%a)
      

   4. Note in the above that to_coo, from_coo are defined so that every
      conceivable  storage representation provides just 2 conversion
      routines, avoiding quadratic explosion. But since all have to
      provide them, the to_coo/from_coo is defined in d_base_mat_mod
      together with d_coo_sparse_mat, which enjoys the "crown prince"
      status with respect to all the other types derived from
      d_base_sparse_mat (its "siblings"). 

   5. How does a user add a new storage format? Very simple. After
      deriving the class and implementing all the necessary methods,
      the user declares in the program a dummy variable of the new
      inner type 

      type(X_YYY_sparse_mat) :: reftype


      then calls 
        call psb_spasb(a,....,mold=reftype)

      In psb_spasb we have
      class(psbn_d_base_sparse_mat), intent(in), optional :: mold

      if (present(mold)) then 
        allocate(temp,source=mold,stat=info)
      end select				 
      call temp%from_coo(a%a)
      call a%a%free()
      call move_alloc(temp,a%a)
  
      AND IT'S DONE! Nothing else in the library requires the explicit
      knowledge of type of  MOLD.