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 :: unknownlist
  REAL(KIND=8) , DIMENSION(:), ALLOCATABLE :: x, rhsloc, rhs
  INTEGER                  :: proc_id, n, ln, p
  INTEGER          , DIMENSION(:), ALLOCATABLE :: ia
  INTEGER                   :: nnz
  INTEGER         , DIMENSION(:), ALLOCATABLE :: ja, mapptr, mapp
  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)

  
  !  We compute a trivial unknown partition for the example 
  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

  ln = pend-pbegin
  ALLOCATE( unknownlist(ln) )
  DO i = pbegin, pend-1
    unknownlist(i-pbegin+1) = i
  END DO		  

  !***************************************************!
  !* ENTER THE GRAPH : PARALLEL INTERFACE            *!
  !* Each processor deals with a part of the 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

  !***************************************************!
  !*            ENTER A USER PARTITION               *!
  !***************************************************!
  !** Only the master processor (0 by default) needs to enter the partition **!
  IF ( proc_id == 0 .AND. argc == 0) THEN 
      !** In this example we enter a non overlaped partition **!
      !** But you can use any overlaped partition **!
    
      ALLOCATE( mapptr(nproc+1) )
      ALLOCATE( mapp(n) )

      i = n/nproc
      DO p = 1, nproc
	mapptr(p) = i*(p-1) + 1 
      END DO	  
      mapptr(nproc+1) = n+1

      DO i = 1, n
	mapp(i) = i
      END DO	   
      CALL HIPS_SETPARTITION(id, nproc, mapptr, mapp, ierr)
      CALL HIPS_ExitOnError(ierr)

      DEALLOCATE(mapptr)
      DEALLOCATE(mapp)
   END IF

 

  job = 1 !**  **!

  nnz = 0
  DO i = pbegin, pend-1
     nnz = nnz + ia( i+1 ) - ia( i )
  END DO
 
  !***********************************
  !*                                 *
  !* ENTER THE MATRIX COEFFICIENTS   *
  !*                                 *
  !***********************************

  !The processor enter the rows pbegin to pend-1 (naive partition) 
  ! In the HIPS_ASSEMBLY_FOOL mode any processor can enter any coefficient :
  !    nevertheless it is better to enter as much as possible coefficients
  !    that are in the local matrix of the processor 
  PRINT * , "Assembly nnz = ",nnz
  CALL HIPS_ASSEMBLYBEGIN(id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL, sym_matrix, ierr)
  CALL HIPS_EXITONERROR(ierr)				
  DO i = pbegin, pend-1
     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( unknownlist )
     DEALLOCATE( rhs )
  ELSE
     rhsloc = 1.
  END IF

  !****************************************************!
  !* Set the rhs corresponding to the user partition  *!
  !****************************************************! 
  CALL HIPS_SETRHS(id, ln, unknownlist, rhsloc, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL, ierr)
  CALL HIPS_EXITONERROR(ierr)				

  !****************************************************!
  !* Get the solution corresponding to user partition *!
  !****************************************************! 
  CALL HIPS_GETSOLUTION(id, ln, unknownlist, x, HIPS_ASSEMBLY_FOOL, 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(unknownlist)
  DEALLOCATE(x)
  DEALLOCATE(rhsloc)

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