HIPS
Hierarchical Iterative Parallel Solver

Documentation

Examples


PROGRAM main

  IMPLICIT NONE
  INCLUDE "mpif.h"
  INCLUDE "hips.inc"

  INTEGER , PARAMETER :: BUFLEN = 200
  CHARACTER(LEN=BUFLEN) ::  argv
  INTEGER               ::  argc

  !* to read parameters *!
  INTEGER  :: sym_pattern, sym_matrix
  CHARACTER(len=BUFLEN) ::  matrixfile
  CHARACTER(len=BUFLEN) ::  rhsfile
  INTEGER  ::  method

  !*  *!

  !*  *!
  INTEGER                  :: id, idnbr, i, j
  INTEGER , DIMENSION(:), ALLOCATABLE :: nodelist
  REAL(KIND=8) , DIMENSION(:), ALLOCATABLE :: x, rhsloc, rhs
  INTEGER                  :: proc_id, n, ln
  INTEGER          , DIMENSION(:), ALLOCATABLE :: ia
  INTEGER                   :: nnz
  INTEGER         , DIMENSION(:), ALLOCATABLE :: ja
  REAL(KIND=8) , DIMENSION(:), ALLOCATABLE :: a
  INTEGER                  :: interior
  INTEGER                  :: domsize, nproc
  INTEGER                  :: pbegin, pend

  INTEGER                                  :: ierror
  INTEGER                                  :: ierr

  INTEGER :: job, mode, m 
  INTEGER  :: k

  !*  *!

  !** Init MPI environment **!
  CALL MPI_Init(ierror)
  CALL MPI_Comm_rank(MPI_COMM_WORLD, proc_id , ierror)
  CALL MPI_Comm_size(MPI_COMM_WORLD, nproc, ierror)

  !***************************************!
  !* Initialize HIPS for one problem     *!
  !***************************************!
  idnbr = 1 !* total *!
  CALL HIPS_INITIALIZE(idnbr, ierr)
  CALL HIPS_EXITONERROR(ierr)

  id = 0 !** id of the linear system **!

  !**************************************************************************!
  !* Read parameter from the file "Inputs"                                  *!
  !* The parameter are set up inside this functions                         *!
  !* this function contains calls HIPS_SetOptionsINT and HIPS_SetOptionREAL *!
  !**************************************************************************!
  CALL HIPS_READOPTIONSFROMFILE(id, sym_pattern, sym_matrix, "", matrixfile, rhsfile, ierr)
  CALL HIPS_EXITONERROR(ierr)

  argc = command_argument_count()
  !** parameter domsize is an argument of testHIPS.ex **!
  IF (argc >= 1) THEN
     CALL GETARG(1,argv)     
     READ(argv , *) domsize
     CALL HIPS_SetOptionINT(id, HIPS_PARTITION_TYPE, 0, ierr)
     CALL HIPS_SetOptionINT(id, HIPS_DOMSIZE, domsize, ierr)
  ELSE
     CALL HIPS_SetOptionINT(id, HIPS_DOMNBR, nproc, ierr)
  END IF

 !**********************************!
 !* Read the matrix from file      *!
 !**********************************!
 !*** In Fortran we first read the sizes of the matrix ***!
 job = 0
 CALL MATRIX_READ(job, n, nnz, ia, ja, a, sym_matrix, TRIM(matrixfile))
 !*** then we allocate memory ***!
 ALLOCATE( ia(n+1),ja(nnz),a(nnz) )
 job = 1
 !*** eventually we set the coefficients ***!
 CALL MATRIX_READ(job, n, nnz, ia, ja, a, sym_matrix, TRIM(matrixfile))
 CALL HIPS_EXITONERROR(ierr)

 CALL HIPS_SetOptionINT (id, HIPS_SYMMETRIC, sym_matrix, ierr)

  !***************************************************!
  !* ENTER THE GRAPH : PARALLEL INTERFACE            *!
  !* Every processor deals with a part of the graph  *!
  !***************************************************!
   !**** Each processor keep a part of the matrix between the row pbegin and pend *****!
  i = n/nproc
  pbegin = i*proc_id
  if (proc_id /= nproc-1) then
     pend   = i*(proc_id+1)
  else
     pend   = n
  end if

  pbegin = pbegin + 1
  pend   = pend   + 1
  PRINT *, "Processor" , proc_id, "enter rows [", pbegin, pend-1, "] of the graph"

  !*************************************!
  !* Enter the adjacency graph         *!
  !*************************************!
  CALL HIPS_GRAPHBEGIN(id, n, ia(pend)-ia(pbegin), ierr)
  CALL HIPS_EXITONERROR(ierr)

  DO i = pbegin,pend-1
     DO j = ia(i), ia(i+1)-1
        CALL HIPS_GRAPHEDGE(id, i, ja(j), ierr)
     END DO
  END DO
  CALL HIPS_GRAPHEND(id, ierr)
  CALL HIPS_EXITONERROR(ierr)				

  IF (sym_pattern /= 0 .AND. sym_matrix /= 1) THEN
    !** Disable the graph symmetrization to save memory **!
     CALL HIPS_SETOPTIONINT(id, HIPS_GRAPH_SYM, 0, ierr)
  END IF

  !***************************************!
  !* Get the ordered list of local nodes *!
  !***************************************!
  CALL HIPS_GETLOCALNODENBR(id, ln, ierr)

  PRINT *,"Processor", proc_id, ": domain total size =", ln, ", overlap =", ln - interior, "nodes"
  ALLOCATE( nodelist(ln) )

  CALL HIPS_GETLOCALNODELIST(id, nodelist, ierr)

  job = 1 !**  **!
  nnz = 0
  DO m = 1, ln
     nnz = nnz + ia( nodelist(m)+1 ) - ia( nodelist(m) )
  END DO
 
  !***********************************
  !*                                 *
  !* ENTER THE MATRIX COEFFICIENTS   *
  !*                                 *
  !***********************************
  ! In this example, each processor enter a row of the matrix. In HIPS_ASSEMBLY_RESPECT
  ! mode you can enter any LOCAL coefficient you like on a processor (according to HIPS_GetUnknownlist).
  ! That is to say only coefficient A(i,j) with i in Unknownlist and j in Unknownlist are taken into 
  ! account. Other coefficient are just ignored.
  PRINT * , "Assembly"
  CALL HIPS_ASSEMBLYBEGIN(id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_RESPECT, sym_matrix, ierr)
  CALL HIPS_EXITONERROR(ierr)				
  DO m = 1, ln
     i = nodelist(m)
     DO k = ia(i) , ia(i+1) - 1 
        CALL HIPS_ASSEMBLYSETVALUE(id, i, ja(k), a(k), ierr)
     END DO
  END DO
  CALL HIPS_ASSEMBLYEND(id, ierr)
  CALL HIPS_EXITONERROR(ierr)				

  DEALLOCATE(ia)
  DEALLOCATE(ja)
  DEALLOCATE(a)
  

  !****************************************!
  !* Set the local right hand side        *!
  !****************************************!
  ALLOCATE( rhsloc(ln) )
  ALLOCATE( x(ln) )
  if (TRIM(rhsfile) /= "0") THEN
     ALLOCATE( rhs(n) )
     CALL VEC_read(n, rhs , rhsfile)
     rhsloc = rhs( nodelist )
     DEALLOCATE( rhs )
  else
     x = 1.
     CALL HIPS_MATRIXVECTORPRODUCT(id, x, rhsloc, ierr)
     CALL HIPS_EXITONERROR(ierr)				
  END if

  !****************************************************!
  !* Set the local rhs                                *!
  !****************************************************! 
  CALL HIPS_SETLOCALRHS(id, rhsloc, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, ierr)
  CALL HIPS_EXITONERROR(ierr)				

  !****************************************************!
  !* Get the local solution                           *!
  !****************************************************! 
  CALL HIPS_GETLOCALSOLUTION(id, x, ierr)
  CALL HIPS_EXITONERROR(ierr)				

  !************************************************!
  !* Free HIPS internal structure for problem id  *!
  !************************************************!
  CALL HIPS_CLEAN(id, ierr)

  !**********************************!
  !* Free HIPS internal structures  *!
  !**********************************!
  CALL HIPS_FINALIZE(ierr)

  DEALLOCATE(nodelist)
  DEALLOCATE(x)
  DEALLOCATE(rhsloc)

  
  !** End MPI **!
   CALL MPI_Finalize(ierror)
  
END PROGRAM main