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  */
  /***************************************************/
   /**** 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)
    pend   = i*(proc_id+1);
  else
    pend   = n;

  fprintf(stdout, "Processor %ld enter rows [%ld %ld] of the graph \n", (long)proc_id, (long)pbegin, (long)pend-1);


  /***************************************************/
  /*                                                 */
  /*            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);
    }

  /***************************************************/
  /*                                                 */
  /*            GET THE LOCAL UNKNOWN LIST           */
  /*                                                 */
  /***************************************************/
  ierr = HIPS_GetLocalUnknownNbr(id, &ln);
  HIPS_ExitOnError(ierr);

  unknownlist = (INTS *)malloc(sizeof(INTS)*ln);

  ierr = HIPS_GetLocalUnknownList(id, unknownlist);
  HIPS_ExitOnError(ierr);


  /***************************************************/
  /*                                                 */
  /*          ENTER THE MATRIX COEFFICIENT           */
  /*                                                 */
  /***************************************************/
 /*  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 coefficients are just ignored. */
  {
    INTS i, m;
    INTL k, nnz;


    nnz = 0;
    for(m=0;m<ln;m++)
      nnz += ia[unknownlist[m]+1] - ia[unknownlist[m]];

    ierr = HIPS_AssemblyBegin(id, nnz, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_RESPECT, sym_matrix);
    HIPS_ExitOnError(ierr);
    
    for(m=0;m<ln;m++)
      {
	i = unknownlist[m];
	for(k=ia[i];k<ia[i+1];k++)
	  {
	    ierr = HIPS_AssemblySetValue(id, i, ja[k], a[k]);
	    HIPS_ExitOnError(ierr);
	  }
      }

    ierr = HIPS_AssemblyEnd(id);
    HIPS_ExitOnError(ierr);


    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];*/
	/*x[i] = 1;*/
	rhsloc[i] = 1;
      /*ierr = HIPS_MatrixVectorProduct(id, x, rhsloc);*/
      HIPS_ExitOnError(ierr);
    }

  ierr = HIPS_SetLocalRHS(id, rhsloc, HIPS_ASSEMBLY_OVW, HIPS_ASSEMBLY_OVW);
  HIPS_ExitOnError(ierr);


  
  /***************************************************/
  /*                                                 */
  /*               SOLVE                             */
  /*                                                 */
  /***************************************************/

  /****************************************************/
  /* Get the local solution                           */
  /****************************************************/ 
  ierr = HIPS_GetLocalSolution(id, x);
  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;
}