Actual source code: fdaij.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/impls/aij/seq/aij.h

      5: EXTERN PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt*[],PetscInt*[],PetscTruth*);
      6: EXTERN PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt*[],PetscInt*[],PetscTruth*);

     10: PetscErrorCode MatFDColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c)
     11: {
     13:   PetscInt       i,*is,n,nrows,N = mat->N,j,k,m,*rows,*ci,*cj,ncols,col;
     14:   PetscInt       nis = iscoloring->n,*rowhit,*columnsforrow,l;
     15:   IS             *isa;
     16:   PetscTruth     done,flg;

     19:   if (!mat->assembled) {
     20:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();");
     21:   }

     23:   ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);
     24:   c->M       = mat->M;  /* set total rows, columns and local rows */
     25:   c->N       = mat->N;
     26:   c->m       = mat->M;
     27:   c->rstart  = 0;

     29:   c->ncolors = nis;
     30:   PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);
     31:   PetscMalloc(nis*sizeof(PetscInt*),&c->columns);
     32:   PetscMalloc(nis*sizeof(PetscInt),&c->nrows);
     33:   PetscMalloc(nis*sizeof(PetscInt*),&c->rows);
     34:   PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);

     36:   /*
     37:       Calls the _SeqAIJ() version of these routines to make sure it does not
     38:      get the reduced (by inodes) version of I and J
     39:   */
     40:   MatGetColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);

     42:   /*
     43:      Temporary option to allow for debugging/testing
     44:   */
     45:   PetscOptionsHasName(PETSC_NULL,"-matfdcoloring_slow",&flg);

     47:   PetscMalloc((N+1)*sizeof(PetscInt),&rowhit);
     48:   PetscMalloc((N+1)*sizeof(PetscInt),&columnsforrow);

     50:   for (i=0; i<nis; i++) {
     51:     ISGetLocalSize(isa[i],&n);
     52:     ISGetIndices(isa[i],&is);
     53:     c->ncolumns[i] = n;
     54:     if (n) {
     55:       PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);
     56:       PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));
     57:     } else {
     58:       c->columns[i]  = 0;
     59:     }

     61:     if (!flg) { /* ------------------------------------------------------------------------------*/
     62:       /* fast, crude version requires O(N*N) work */
     63:       PetscMemzero(rowhit,N*sizeof(PetscInt));
     64:       /* loop over columns*/
     65:       for (j=0; j<n; j++) {
     66:         col  = is[j];
     67:         rows = cj + ci[col];
     68:         m    = ci[col+1] - ci[col];
     69:         /* loop over columns marking them in rowhit */
     70:         for (k=0; k<m; k++) {
     71:           rowhit[*rows++] = col + 1;
     72:         }
     73:       }
     74:       /* count the number of hits */
     75:       nrows = 0;
     76:       for (j=0; j<N; j++) {
     77:         if (rowhit[j]) nrows++;
     78:       }
     79:       c->nrows[i] = nrows;
     80:       PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);
     81:       PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);
     82:       nrows       = 0;
     83:       for (j=0; j<N; j++) {
     84:         if (rowhit[j]) {
     85:           c->rows[i][nrows]          = j;
     86:           c->columnsforrow[i][nrows] = rowhit[j] - 1;
     87:           nrows++;
     88:         }
     89:       }
     90:     } else {  /*-------------------------------------------------------------------------------*/
     91:       /* slow version, using rowhit as a linked list */
     92:       PetscInt currentcol,fm,mfm;
     93:       rowhit[N] = N;
     94:       nrows     = 0;
     95:       /* loop over columns */
     96:       for (j=0; j<n; j++) {
     97:         col   = is[j];
     98:         rows  = cj + ci[col];
     99:         m     = ci[col+1] - ci[col];
    100:         /* loop over columns marking them in rowhit */
    101:         fm    = N; /* fm points to first entry in linked list */
    102:         for (k=0; k<m; k++) {
    103:           currentcol = *rows++;
    104:           /* is it already in the list? */
    105:           do {
    106:             mfm  = fm;
    107:             fm   = rowhit[fm];
    108:           } while (fm < currentcol);
    109:           /* not in list so add it */
    110:           if (fm != currentcol) {
    111:             nrows++;
    112:             columnsforrow[currentcol] = col;
    113:             /* next three lines insert new entry into linked list */
    114:             rowhit[mfm]               = currentcol;
    115:             rowhit[currentcol]        = fm;
    116:             fm                        = currentcol;
    117:             /* fm points to present position in list since we know the columns are sorted */
    118:           } else {
    119:             SETERRQ(PETSC_ERR_PLIB,"Detected invalid coloring");
    120:           }

    122:         }
    123:       }
    124:       c->nrows[i] = nrows;
    125:       PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);
    126:       PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);
    127:       /* now store the linked list of rows into c->rows[i] */
    128:       nrows       = 0;
    129:       fm          = rowhit[N];
    130:       do {
    131:         c->rows[i][nrows]            = fm;
    132:         c->columnsforrow[i][nrows++] = columnsforrow[fm];
    133:         fm                           = rowhit[fm];
    134:       } while (fm < N);
    135:     } /* ---------------------------------------------------------------------------------------*/
    136:     ISRestoreIndices(isa[i],&is);
    137:   }
    138:   MatRestoreColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);

    140:   PetscFree(rowhit);
    141:   PetscFree(columnsforrow);

    143:   /* Optimize by adding the vscale, and scaleforrow[][] fields */
    144:   /*
    145:        see the version for MPIAIJ
    146:   */
    147:   VecCreateGhost(mat->comm,mat->m,PETSC_DETERMINE,0,PETSC_NULL,&c->vscale);CHKERRQ(ierr)
    148:   PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);
    149:   for (k=0; k<c->ncolors; k++) {
    150:     PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);
    151:     for (l=0; l<c->nrows[k]; l++) {
    152:       col = c->columnsforrow[k][l];
    153:       c->vscaleforrow[k][l] = col;
    154:     }
    155:   }
    156:   ISColoringRestoreIS(iscoloring,&isa);
    157:   return(0);
    158: }