HIPS
Hierarchical Iterative Parallel Solver

Documentation

Examples


PROGRAM main

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

  INTEGER , PARAMETER :: BUFLEN = 200

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

  !*  *!

  !*  *!
  INTEGER                  :: id, idnbr, i
  REAL(KIND=8) , DIMENSION(:), ALLOCATABLE :: xx, rhs
  INTEGER                  :: proc_id, n, ln
  INTEGER          , DIMENSION(:), ALLOCATABLE :: ia
  INTEGER                   :: nnz
  INTEGER         , DIMENSION(:), ALLOCATABLE :: ja
  REAL(KIND=8) , DIMENSION(:), ALLOCATABLE :: a
  INTEGER                  :: domsize, nproc

  INTEGER                  :: job
  INTEGER                                  :: ierr
  INTEGER                                  :: ierror
  !*  *!

  !** 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
      !** parameter domsize is an argument of testHIPS.ex **!
     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      *!
  !**********************************!
  IF (proc_id == 0) THEN
     !*** 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 ***!
     PRINT *, "N = ",n, "NNZ =", nnz
     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))
  END IF

  !*** Everyone need to know if the matrix is symmetric ***!
  CALL MPI_Bcast(sym_pattern, 1, MPI_INTEGER4, 0, MPI_COMM_WORLD,ierr)
  CALL MPI_Bcast(sym_matrix, 1, MPI_INTEGER4, 0, MPI_COMM_WORLD ,ierr)

  CALL HIPS_SETOPTIONINT (id, HIPS_SYMMETRIC, sym_matrix, 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

  IF (proc_id == 0) THEN
     PRINT *,"Matrix : dim=", n , " nnz=",nnz
  END IF

  !***************************************************!
  !*                                                 *!
  !* ENTER THE MATRIX GRAPH : CENTRALIZED INTERFACE  *!
  !*                                                 *!
  !***************************************************!
  CALL HIPS_GRAPHGLOBALCSR(id, n, ia, ja, 0, ierr)
  CALL HIPS_EXITONERROR(ierr)

  !***********************************
  !*                                 *
  !* ENTER THE MATRIX COEFFICIENTS   *
  !*                                 *
  !***********************************
  CALL HIPS_MATRIXGLOBALCSR(id, n, ia, ja, a,  0, HIPS_ASSEMBLY_OVW, sym_matrix, ierr)
  CALL HIPS_EXITONERROR(ierr)

  IF (proc_id == 0) THEN
     DEALLOCATE(ia)
     DEALLOCATE(ja)
     DEALLOCATE(a)
  END IF

     
  !****************************************!
  !* Set the right hand side (from proc 0)*!
  !****************************************!
  
  IF (proc_id == 0) THEN
     ALLOCATE(rhs(n))
     IF (TRIM(rhsfile) /= "0") THEN
        CALL VEC_read(n, rhs, TRIM(rhsfile) )
     ELSE
        rhs = 1 
     END IF
  END IF

  !****************************************************!
  !* Set the global rhs                               *!
  !* rhs is only significant on the master processor  *!
  !****************************************************! 
  CALL HIPS_SETGLOBALRHS(id, rhs, 0, 0, ierr)
  CALL HIPS_EXITONERROR(ierr);

  !****************************************************!
  !* Get the global solution on processor 0           *!
  !* Original ordering                                *!
  !****************************************************! 
  IF (proc_id == 0) THEN
     ALLOCATE( xx(n) )
  END IF

  CALL HIPS_GETGLOBALSOLUTION(id, xx, 0, 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)
  CALL HIPS_EXITONERROR(ierr)

  IF (proc_id == 0) THEN
      DEALLOCATE(xx)
      DEALLOCATE(rhs)
   END IF

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