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" /** io.h requires hips.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;
  /*INTS numflag;*/
  COEF *xx=NULL, *rhs;
  INTS proc_id, n;
  INTL *ia, nnz;
  INTS *ja;
  COEF *a;
  INTS domsize, nproc;
  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 parameter from the file "Inputs"                                  */
  /* The parameter are set up inside this functions                         */
  /* this function contains calls HIPS_SetOptionsINT and HIPS_SetOptionREAL */
  /**************************************************************************/
  ierr = HIPS_ReadOptionsFromFile(id, NULL, &sym_pattern, &sym_matrix, matrixfile, rhsfile);
  HIPS_ExitOnError(ierr);

  if(argc >= 2)
    {
      /** parameter domsize is an argument of testHIPS.ex **/
      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      */
  /**********************************/
  if(proc_id == 0)
    {
      CSRread(matrixfile, &n, &nnz, &ia, &ja, &a, &sym_matrix);
      dumpcsr(stderr, a, ja, ia, n);
    }

  /*** Everyone need to know if the matrix is symmetric ***/
  MPI_Bcast(&sym_pattern, EE(1), COMM_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&sym_matrix, EE(1), COMM_INT, 0, MPI_COMM_WORLD);

  HIPS_SetOptionINT (id, HIPS_SYMMETRIC, sym_matrix);

  
  if(sym_pattern != 0 && sym_matrix != 1)
    /** Disable the graph symmetrization to save memory **/
    HIPS_SetOptionINT(id, HIPS_GRAPH_SYM, 0);

  if (proc_id == 0)
    fprintf(stdout, "Matrix : dim=%ld nnz=%ld\n", (long)n, (long)nnz);

  /** C : numbering starts from 0 **/
  if(proc_id == 0)
    CSR_Fnum2Cnum(ja, ia, n);

  HIPS_SetOptionINT(id, HIPS_FORTRAN_NUMBERING, 0);

  
  /***************************************************/
  /*                                                 */
  /* ENTER THE MATRIX GRAPH : CENTRALIZED INTERFACE  */
  /*                                                 */
  /***************************************************/
  ierr = HIPS_GraphGlobalCSR(id, n, ia, ja, 0);
  HIPS_ExitOnError(ierr);

  /**********************************************************/
  /*                                                        */
  /* ENTER THE MATRIX COEFFICIENTS : CENTRALIZED INTERFACE  */
  /*                                                        */
  /**********************************************************/
  ierr = HIPS_MatrixGlobalCSR(id, n, ia, ja, a, 0, HIPS_ASSEMBLY_OVW, sym_matrix);
  HIPS_ExitOnError(ierr);

  if(proc_id == 0)
    {
      free(ia);
      free(ja);
      free(a);
    }
  


  /*fprintf(stdout, "Transpose the matrix \n");
    HIPS_TransposeMatrix(id);*/
  




  /****************************************/
  /* Set the right hand side (from proc 0)*/
  /****************************************/
  if(proc_id == 0)
    {
      rhs = (COEF *)malloc(sizeof(COEF)*n);
      if(strcmp(rhsfile, "0") != 0)
	VECread(rhsfile, n, rhs);
      else
	for(i=0;i<n;i++)
	  rhs[i] = 1; 
    }
  else
    rhs = NULL;
  
 
  /****************************************************/
  /* Set the global rhs                               */
  /* rhs is only significant on the master processor  */
  /****************************************************/ 
  ierr = HIPS_SetGlobalRHS(id, rhs, 0, 0);
  HIPS_ExitOnError(ierr);
  



  /****************************************************/
  /* Get the global solution on processor 0           */
  /* Original ordering                                */
  /****************************************************/ 
  if(proc_id == 0)
    xx = (COEF *)malloc(sizeof(COEF)*n);

  ierr = HIPS_GetGlobalSolution(id, xx, 0);
  HIPS_ExitOnError(ierr);

  {
    int i;
    for(i=0;i<n;i++)
      fprintf(stdout, "x[%d] = %g \n", i, xx[i]);
  }
 

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

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


  if(proc_id == 0)
    {
      free(xx);
      free(rhs);
    }

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