Actual source code: metisnd.c


  2: #include <petscmat.h>
  3: #include <petsc/private/matorderimpl.h>
  4: #include <metis.h>

  6: /*
  7:     MatGetOrdering_METISND - Find the nested dissection ordering of a given matrix.
  8: */
  9: PETSC_INTERN PetscErrorCode MatGetOrdering_METISND(Mat mat,MatOrderingType type,IS *row,IS *col)
 10: {
 12:   PetscInt       i,j,iptr,ival,nrow,*xadj,*adjncy,*perm,*iperm;
 13:   const PetscInt *ia,*ja;
 14:   int            status;
 15:   Mat            B = NULL;
 16:   idx_t          options[METIS_NOPTIONS];
 17:   PetscBool      done;

 19:   MatGetRowIJ(mat,0,PETSC_TRUE,PETSC_TRUE,&nrow,&ia,&ja,&done);
 20:   if (!done) {
 21:     MatConvert(mat,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
 22:     MatGetRowIJ(B,0,PETSC_TRUE,PETSC_TRUE,&nrow,&ia,&ja,&done);
 23:   }
 24:   METIS_SetDefaultOptions(options);
 25:   options[METIS_OPTION_NUMBERING] = 0;
 26:   PetscOptionsBegin(PetscObjectComm((PetscObject)mat),((PetscObject)mat)->prefix,"METISND Options","Mat");

 28:   ival = (PetscInt)options[METIS_OPTION_NSEPS];
 29:   PetscOptionsInt("-mat_ordering_metisnd_nseps","number of different separators per level","None",ival,&ival,NULL);
 30:   options[METIS_OPTION_NSEPS] = (idx_t)ival;

 32:   ival = (PetscInt)options[METIS_OPTION_NITER];
 33:   PetscOptionsInt("-mat_ordering_metisnd_niter","number of refinement iterations","None",ival,&ival,NULL);
 34:   options[METIS_OPTION_NITER] = (idx_t)ival;

 36:   ival = (PetscInt)options[METIS_OPTION_UFACTOR];
 37:   PetscOptionsInt("-mat_ordering_metisnd_ufactor","maximum allowed imbalance","None",ival,&ival,NULL);
 38:   options[METIS_OPTION_UFACTOR] = (idx_t)ival;

 40:   ival = (PetscInt)options[METIS_OPTION_PFACTOR];
 41:   PetscOptionsInt("-mat_ordering_metisnd_pfactor","minimum degree of vertices that will be ordered last","None",ival,&ival,NULL);
 42:   options[METIS_OPTION_PFACTOR] = (idx_t)ival;

 44:   PetscOptionsEnd();

 46:   PetscMalloc4(nrow+1,&xadj,ia[nrow],&adjncy,nrow,&perm,nrow,&iperm);
 47:   /* The adjacency list of a vertex should not contain the vertex itself.
 48:   */
 49:   iptr = 0;
 50:   xadj[iptr] = 0;
 51:   for (j=0; j<nrow; j++) {
 52:     for (i=ia[j]; i<ia[j+1]; i++) {
 53:       if (ja[i] != j)
 54:         adjncy[iptr++] = ja[i];
 55:     }
 56:     xadj[j+1] = iptr;
 57:   }

 59:   status = METIS_NodeND(&nrow,(idx_t*)xadj,(idx_t*)adjncy,NULL,options,(idx_t*)perm,(idx_t*)iperm);
 60:   switch (status) {
 61:   case METIS_OK: break;
 62:   case METIS_ERROR:
 63:     SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_LIB,"METIS returned with an unspecified error");
 64:   case METIS_ERROR_INPUT:
 65:     SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_LIB,"METIS received an invalid input");
 66:   case METIS_ERROR_MEMORY:
 67:     SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_MEM,"METIS could not compute ordering");
 68:   default:
 69:     SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_LIB,"Unexpected return value");
 70:   }

 72:   if (B) {
 73:     MatRestoreRowIJ(B,0,PETSC_TRUE,PETSC_TRUE,NULL,&ia,&ja,&done);
 74:     MatDestroy(&B);
 75:   } else {
 76:     MatRestoreRowIJ(mat,0,PETSC_TRUE,PETSC_TRUE,NULL,&ia,&ja,&done);
 77:   }

 79:   ISCreateGeneral(PETSC_COMM_SELF,nrow,perm,PETSC_COPY_VALUES,row);
 80:   ISCreateGeneral(PETSC_COMM_SELF,nrow,perm,PETSC_COPY_VALUES,col);
 81:   PetscFree4(xadj,adjncy,perm,iperm);
 82:   return 0;
 83: }