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: }