Actual source code: sorder.c

      1: #define PETSCMAT_DLL

      3: /*
      4:      Provides the code that allows PETSc users to register their own
      5:   sequential matrix Ordering routines.
      6: */
     7:  #include src/mat/matimpl.h
     8:  #include petscsys.h

     10: PetscFList      MatOrderingList = 0;
     11: PetscTruth MatOrderingRegisterAllCalled = PETSC_FALSE;

     13: EXTERN PetscErrorCode MatOrdering_Flow_SeqAIJ(Mat,const MatOrderingType,IS *,IS *);

     17: PetscErrorCode MatOrdering_Flow(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
     18: {
     20:   SETERRQ(PETSC_ERR_SUP,"Cannot do default flow ordering for matrix type");
     21: #if !defined(PETSC_USE_DEBUG)
     22:   return(0);
     23: #endif
     24: }

     29: PetscErrorCode PETSCMAT_DLLEXPORT MatOrdering_Natural(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
     30: {
     32:   PetscInt       n,i,*ii;
     33:   PetscTruth     done;
     34:   MPI_Comm       comm;

     37:   PetscObjectGetComm((PetscObject)mat,&comm);
     38:   MatGetRowIJ(mat,0,PETSC_FALSE,&n,PETSC_NULL,PETSC_NULL,&done);
     39:   MatRestoreRowIJ(mat,0,PETSC_FALSE,&n,PETSC_NULL,PETSC_NULL,&done);
     40:   if (done) { /* matrix may be "compressed" in symbolic factorization, due to i-nodes or block storage */
     41:     /*
     42:       We actually create general index sets because this avoids mallocs to
     43:       to obtain the indices in the MatSolve() routines.
     44:       ISCreateStride(PETSC_COMM_SELF,n,0,1,irow);
     45:       ISCreateStride(PETSC_COMM_SELF,n,0,1,icol);
     46:     */
     47:     PetscMalloc(n*sizeof(PetscInt),&ii);
     48:     for (i=0; i<n; i++) ii[i] = i;
     49:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,irow);
     50:     ISCreateGeneral(PETSC_COMM_SELF,n,ii,icol);
     51:     PetscFree(ii);
     52:   } else {
     53:     PetscInt start,end;

     55:     MatGetOwnershipRange(mat,&start,&end);
     56:     ISCreateStride(comm,end-start,start,1,irow);
     57:     ISCreateStride(comm,end-start,start,1,icol);
     58:   }
     59:   ISSetIdentity(*irow);
     60:   ISSetIdentity(*icol);
     61:   return(0);
     62: }

     66: /*
     67:      Orders the rows (and columns) by the lengths of the rows.
     68:    This produces a symmetric Ordering but does not require a
     69:    matrix with symmetric non-zero structure.
     70: */
     73: PetscErrorCode PETSCMAT_DLLEXPORT MatOrdering_RowLength(Mat mat,const MatOrderingType type,IS *irow,IS *icol)
     74: {
     76:   PetscInt       n,*ia,*ja,*permr,*lens,i;
     77:   PetscTruth     done;

     80:   MatGetRowIJ(mat,0,PETSC_FALSE,&n,&ia,&ja,&done);
     81:   if (!done) SETERRQ(PETSC_ERR_SUP,"Cannot get rows for matrix");

     83:   PetscMalloc(2*n*sizeof(PetscInt),&lens);
     84:   permr = lens + n;
     85:   for (i=0; i<n; i++) {
     86:     lens[i]  = ia[i+1] - ia[i];
     87:     permr[i] = i;
     88:   }
     89:   MatRestoreRowIJ(mat,0,PETSC_FALSE,&n,&ia,&ja,&done);

     91:   PetscSortIntWithPermutation(n,lens,permr);

     93:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,irow);
     94:   ISCreateGeneral(PETSC_COMM_SELF,n,permr,icol);
     95:   PetscFree(lens);
     96:   return(0);
     97: }

    102: PetscErrorCode PETSCMAT_DLLEXPORT MatOrderingRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(Mat,const MatOrderingType,IS*,IS*))
    103: {
    105:   char           fullname[PETSC_MAX_PATH_LEN];

    108:   PetscFListConcat(path,name,fullname);
    109:   PetscFListAdd(&MatOrderingList,sname,fullname,(void (*)(void))function);
    110:   return(0);
    111: }

    115: /*@C
    116:    MatOrderingRegisterDestroy - Frees the list of ordering routines.

    118:    Not collective

    120:    Level: developer
    121:
    122: .keywords: matrix, register, destroy

    124: .seealso: MatOrderingRegisterDynamic(), MatOrderingRegisterAll()
    125: @*/
    126: PetscErrorCode PETSCMAT_DLLEXPORT MatOrderingRegisterDestroy(void)
    127: {

    131:   if (MatOrderingList) {
    132:     PetscFListDestroy(&MatOrderingList);
    133:     MatOrderingList = 0;
    134:   }
    135:   return(0);
    136: }

    138:  #include src/mat/impls/aij/mpi/mpiaij.h
    141: /*@C
    142:    MatGetOrdering - Gets a reordering for a matrix to reduce fill or to
    143:    improve numerical stability of LU factorization.

    145:    Collective on Mat

    147:    Input Parameters:
    148: +  mat - the matrix
    149: -  type - type of reordering, one of the following:
    150: $      MATORDERING_NATURAL - Natural
    151: $      MATORDERING_ND - Nested Dissection
    152: $      MATORDERING_1WD - One-way Dissection
    153: $      MATORDERING_RCM - Reverse Cuthill-McKee
    154: $      MATORDERING_QMD - Quotient Minimum Degree

    156:    Output Parameters:
    157: +  rperm - row permutation indices
    158: -  cperm - column permutation indices

    161:    Options Database Key:
    162: . -mat_view_ordering_draw - plots matrix nonzero structure in new ordering

    164:    Level: intermediate
    165:
    166:    Notes:
    167:       This DOES NOT actually reorder the matrix; it merely returns two index sets
    168:    that define a reordering. This is usually not used directly, rather use the
    169:    options PCLUSetMatOrdering() or PCILUSetMatOrdering().

    171:    The user can define additional orderings; see MatOrderingRegisterDynamic().

    173: .keywords: matrix, set, ordering, factorization, direct, ILU, LU,
    174:            fill, reordering, natural, Nested Dissection,
    175:            One-way Dissection, Cholesky, Reverse Cuthill-McKee,
    176:            Quotient Minimum Degree

    178: .seealso:   MatOrderingRegisterDynamic(), PCLUSetMatOrdering(), PCILUSetMatOrdering()
    179: @*/
    180: PetscErrorCode PETSCMAT_DLLEXPORT MatGetOrdering(Mat mat,const MatOrderingType type,IS *rperm,IS *cperm)
    181: {
    182:   PetscErrorCode  ierr;
    183:   PetscInt        mmat,nmat,mis,m;
    184:   PetscErrorCode (*r)(Mat,const MatOrderingType,IS*,IS*);
    185:   PetscTruth     flg,isseqdense,ismpidense;

    191:   if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix");
    192:   if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

    194:   PetscTypeCompare((PetscObject)mat,MATSEQDENSE,&isseqdense);
    195:   PetscTypeCompare((PetscObject)mat,MATMPIDENSE,&ismpidense);
    196:   if (isseqdense || ismpidense) {
    197:     MatGetLocalSize(mat,&m,PETSC_NULL);
    198:     /*
    199:        Dense matrices only give natural ordering
    200:     */
    201:     ISCreateStride(PETSC_COMM_SELF,0,m,1,cperm);
    202:     ISCreateStride(PETSC_COMM_SELF,0,m,1,rperm);
    203:     ISSetIdentity(*cperm);
    204:     ISSetIdentity(*rperm);
    205:     ISSetPermutation(*rperm);
    206:     ISSetPermutation(*cperm);
    207:     return(0);
    208:   }

    210:   if (!mat->M) { /* matrix has zero rows */
    211:     ISCreateStride(PETSC_COMM_SELF,0,0,1,cperm);
    212:     ISCreateStride(PETSC_COMM_SELF,0,0,1,rperm);
    213:     ISSetIdentity(*cperm);
    214:     ISSetIdentity(*rperm);
    215:     ISSetPermutation(*rperm);
    216:     ISSetPermutation(*cperm);
    217:     return(0);
    218:   }

    220:   if (!MatOrderingRegisterAllCalled) {
    221:     MatOrderingRegisterAll(PETSC_NULL);
    222:   }

    224:   PetscLogEventBegin(MAT_GetOrdering,mat,0,0,0);
    225:    PetscFListFind(mat->comm,MatOrderingList,type,(void (**)(void)) &r);
    226:   if (!r) {SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Unknown or unregistered type: %s",type);}

    228:   (*r)(mat,type,rperm,cperm);
    229:   ISSetPermutation(*rperm);
    230:   ISSetPermutation(*cperm);

    232:   /*
    233:       Adjust for inode (reduced matrix ordering) only if row permutation
    234:     is smaller then matrix size
    235:   */
    236:   MatGetLocalSize(mat,&mmat,&nmat);
    237:   ISGetLocalSize(*rperm,&mis);
    238:   if (mmat > mis) {
    239:     MatInodeAdjustForInodes(mat,rperm,cperm);
    240:   }

    242:   PetscLogEventEnd(MAT_GetOrdering,mat,0,0,0);

    244:   PetscOptionsHasName(PETSC_NULL,"-mat_view_ordering_draw",&flg);
    245:   if (flg) {
    246:     Mat tmat;
    247:     PetscOptionsHasName(PETSC_NULL,"-mat_view_contour",&flg);
    248:     if (flg) {
    249:       PetscViewerPushFormat(PETSC_VIEWER_DRAW_(mat->comm),PETSC_VIEWER_DRAW_CONTOUR);
    250:     }
    251:     MatPermute(mat,*rperm,*cperm,&tmat);
    252:     MatView(tmat,PETSC_VIEWER_DRAW_(mat->comm));
    253:     if (flg) {
    254:       PetscViewerPopFormat(PETSC_VIEWER_DRAW_(mat->comm));
    255:     }
    256:     MatDestroy(tmat);
    257:   }

    259:   return(0);
    260: }

    264: PetscErrorCode MatGetOrderingList(PetscFList *list)
    265: {
    267:   *list = MatOrderingList;
    268:   return(0);
    269: }