HIPS
Hierarchical Iterative Parallel Solver

Documentation

Examples




#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include "math.h"

#include "hips.h"
#include "io.h"

#define BUFLEN 200

int main(int argc, char *argv[])
{  

  /* to read parameters */
  INTS  sym_pattern, sym_matrix;
  char matrixfile[BUFLEN];
  char rhsfile   [BUFLEN];
  /*  */

  /*  */
  INTS id, idnbr, i, j;
  INTS *unknownlist;
  COEF *x, *rhsloc;
  INTS proc_id, n, ln;
  INTL *ia, nnz;
  INTS *ja;
  COEF *a;
  INTS domsize, nproc;
  INTS pbegin, pend;
  INTS ierr;
  /*  */

  /*  */

  /** Init MPI environment **/
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
  MPI_Comm_size(MPI_COMM_WORLD, &nproc);

  /***************************************/
  /* Initialize HIPS for one problem     */
  /***************************************/
  idnbr = 1; /* total */
  ierr = HIPS_Initialize(idnbr);
  HIPS_ExitOnError(ierr);

  id = 0; /** id of the linear system **/
  /******************************************/
  /* Read the parameters from "inputs" file */
  /******************************************/
  ierr = HIPS_ReadOptionsFromFile(id, NULL, &sym_pattern, &sym_matrix, matrixfile, rhsfile);
  HIPS_ExitOnError(ierr);
  
  /** parameter domsize is an argument of testHIPS.ex **/
  if(argc >= 2)
    {
      domsize = atoi(argv[1]);
      HIPS_SetOptionINT(id, HIPS_PARTITION_TYPE, 0);
      HIPS_SetOptionINT(id, HIPS_DOMSIZE, domsize);
    }
  else
    {
      HIPS_SetOptionINT(id, HIPS_DOMNBR, nproc);
    }


  /**********************************/
  /* Read the matrix from file      */
  /**********************************/
  CSRread(matrixfile, &n, &nnz, &ia, &ja, &a, &sym_matrix);
  HIPS_SetOptionINT (id, HIPS_SYMMETRIC, sym_matrix);


  /** C : numbering starts from 0 **/
  CSR_Fnum2Cnum(ja, ia, n);
  HIPS_SetOptionINT(id, HIPS_FORTRAN_NUMBERING, 0);
 
  /***************************************************/
  /* ENTER THE GRAPH : PARALLEL INTERFACE            */
  /* Every processor deals with a part of the graph  */
  /***************************************************/
  
  i = n/nproc;
  pbegin = i*proc_id;
  if(proc_id != nproc-1)
    pend   = i*(proc_id+1);
  else
    pend   = n;
  
  /*** We compute a trivial unknown partition for the example **/
  ln = pend-pbegin;
  unknownlist = (INTS *)malloc(sizeof(INTS)*ln);
  for(i=0;i<ln;i++)
    unknownlist[i] = pbegin+i;
  
  /***************************************************/
  /*                                                 */
  /*            ENTER THE GRAPH                      */
  /*                                                 */
  /***************************************************/
  ierr = HIPS_GraphBegin(id, n, ia[pend]-ia[pbegin]);
  HIPS_ExitOnError(ierr);

  for(i=pbegin;i<pend;i++)
    for(j=ia[i];j<ia[i+1];j++)
      {
	ierr = HIPS_GraphEdge(id, i, ja[j]);
	HIPS_ExitOnError(ierr);
      }

  ierr = HIPS_GraphEnd(id);
  HIPS_ExitOnError(ierr);

  if(sym_pattern != 0 && sym_matrix != 1)
    {
      /** Disable the graph symmetrization (it saves memory).. **/
      /** If you are not sure that the graph is symmetric or not DO NOT  set this parameter to 0 **/
      HIPS_SetOptionINT(id, HIPS_GRAPH_SYM, 0);
    }

  /***************************************************/
  /*                                                 */
  /*            ENTER A USER PARTITION               */
  /*                                                 */
  /***************************************************/
  /** Only the master processor (0 by default) needs to enter the partition **/
  if(proc_id == 0 && argc < 2)
    {
      /** In this example we enter a non overlaped partition **/
      /** But you can use any overlaped partition **/
      INTS *mapptr, *mapp;
      INTS p;
      mapptr = (INTS *)malloc(sizeof(INTS)*(nproc+1));
      mapp = (INTS *)malloc(sizeof(INTS)*n);

      i = n/nproc;
      for(p=0;p<nproc;p++)
	mapptr[p] = i*p;
      mapptr[nproc] = n;

      for(i=0;i<n;i++)
	mapp[i] = i;
      
      /** If a domsize was specified then we cannot use this **/
      ierr = HIPS_SetPartition(id, nproc, mapptr, mapp);
      HIPS_ExitOnError(ierr);

      free(mapptr);
      free(mapp);
    }

  /***************************************************/
  /*                                                 */
  /*          ENTER THE MATRIX COEFFICIENT           */
  /*                                                 */
  /***************************************************/
  /** 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 **/
  nnz = 0;
  for(i=pbegin;i<pend;i++)
    nnz += ia[i+1] - ia[i];

  ierr = HIPS_AssemblyBegin(id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL, sym_matrix);
  HIPS_ExitOnError(ierr);
    
  for(i=pbegin;i<pend;i++)
    for(j=ia[i];j<ia[i+1];j++)
      {
	ierr = HIPS_AssemblySetValue(id, i, ja[j], a[j]);
	HIPS_ExitOnError(ierr);
      }
  ierr = HIPS_AssemblyEnd(id);
  HIPS_ExitOnError(ierr);


  /** Do not need this anymore **/
  free(ia);
  free(ja);
  free(a);


  /***************************************************/
  /*                                                 */
  /*          ENTER THE RIGHT-HAND-SIDE              */
  /*                                                 */
  /***************************************************/
  rhsloc = (COEF *)malloc(sizeof(COEF)*ln);
  x = (COEF *)malloc(sizeof(COEF)*ln);
  if(strcmp(rhsfile, "0") != 0)
    {
      COEF* rhs = (COEF *)malloc(sizeof(COEF)*n);
      VECread(rhsfile, n, rhs);
      for(i=0;i<ln;i++)
	rhsloc[i] = rhs[unknownlist[i]];
      free(rhs);
    }
  else
    {
      for(i=0;i<ln;i++)
	/*x[i] = unknownlist[i];*/
	rhsloc[i] = 1;
      /*ierr = HIPS_MatrixVectorProduct(id, x, rhsloc);
	HIPS_ExitOnError(ierr);*/
    }

  ierr = HIPS_SetRHS(id, ln, unknownlist, rhsloc, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_FOOL);
  HIPS_ExitOnError(ierr);

  
  /***************************************************/
  /*                                                 */
  /*               SOLVE                             */
  /*                                                 */
  /***************************************************/
  /****************************************************/
  /* Get the local solution                           */
  /****************************************************/ 
  ierr = HIPS_GetSolution(id, ln, unknownlist, x, HIPS_ASSEMBLY_FOOL);
  HIPS_ExitOnError(ierr);
  

  /**************************************************/
  /* Free HIPS internal structure for problem "id"  */
  /**************************************************/
  ierr = HIPS_Clean(id);
  HIPS_ExitOnError(ierr);


  /**********************************/
  /* Free HIPS internal structures  */
  /**********************************/
  ierr = HIPS_Finalize();
  HIPS_ExitOnError(ierr);

  free(unknownlist);
  free(x);
  free(rhsloc);

  /** End MPI **/
  MPI_Finalize();
  
  return 0;
}