1 #define PETSCMAT_DLL 2 3 #include "../src/mat/impls/aij/seq/aij.h" 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "MatFDColoringCreate_SeqAIJ" 7 /* 8 This routine is shared by AIJ and BAIJ matrices, since it operators only on the nonzero structure of the elements or blocks. 9 This is why it has the ugly code with the MatGetBlockSize() 10 */ 11 PetscErrorCode MatFDColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 12 { 13 PetscErrorCode ierr; 14 PetscInt i,n,nrows,N,j,k,m,*rows,*ci,*cj,ncols,col; 15 const PetscInt *is; 16 PetscInt nis = iscoloring->n,*rowhit,*columnsforrow,l,bs = 1; 17 IS *isa; 18 PetscTruth done,flg = PETSC_FALSE; 19 PetscTruth flg1,flg2; 20 21 PetscFunctionBegin; 22 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();"); 23 24 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 25 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 26 ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 27 ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 28 if (flg1 || flg2) { 29 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 30 } 31 32 N = mat->cmap->N/bs; 33 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 34 c->N = mat->cmap->N/bs; 35 c->m = mat->rmap->N/bs; 36 c->rstart = 0; 37 38 c->ncolors = nis; 39 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 40 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 41 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 42 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 43 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 44 45 ierr = MatGetColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 46 if (!done) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name); 47 48 /* 49 Temporary option to allow for debugging/testing 50 */ 51 ierr = PetscOptionsGetTruth(PETSC_NULL,"-matfdcoloring_slow",&flg,PETSC_NULL);CHKERRQ(ierr); 52 53 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 54 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 55 56 for (i=0; i<nis; i++) { 57 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 58 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 59 c->ncolumns[i] = n; 60 if (n) { 61 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 62 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 63 } else { 64 c->columns[i] = 0; 65 } 66 67 if (!flg) { /* ------------------------------------------------------------------------------*/ 68 /* fast, crude version requires O(N*N) work */ 69 ierr = PetscMemzero(rowhit,N*sizeof(PetscInt));CHKERRQ(ierr); 70 /* loop over columns*/ 71 for (j=0; j<n; j++) { 72 col = is[j]; 73 rows = cj + ci[col]; 74 m = ci[col+1] - ci[col]; 75 /* loop over columns marking them in rowhit */ 76 for (k=0; k<m; k++) { 77 rowhit[*rows++] = col + 1; 78 } 79 } 80 /* count the number of hits */ 81 nrows = 0; 82 for (j=0; j<N; j++) { 83 if (rowhit[j]) nrows++; 84 } 85 c->nrows[i] = nrows; 86 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 87 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 88 nrows = 0; 89 for (j=0; j<N; j++) { 90 if (rowhit[j]) { 91 c->rows[i][nrows] = j; 92 c->columnsforrow[i][nrows] = rowhit[j] - 1; 93 nrows++; 94 } 95 } 96 } else { /*-------------------------------------------------------------------------------*/ 97 /* slow version, using rowhit as a linked list */ 98 PetscInt currentcol,fm,mfm; 99 rowhit[N] = N; 100 nrows = 0; 101 /* loop over columns */ 102 for (j=0; j<n; j++) { 103 col = is[j]; 104 rows = cj + ci[col]; 105 m = ci[col+1] - ci[col]; 106 /* loop over columns marking them in rowhit */ 107 fm = N; /* fm points to first entry in linked list */ 108 for (k=0; k<m; k++) { 109 currentcol = *rows++; 110 /* is it already in the list? */ 111 do { 112 mfm = fm; 113 fm = rowhit[fm]; 114 } while (fm < currentcol); 115 /* not in list so add it */ 116 if (fm != currentcol) { 117 nrows++; 118 columnsforrow[currentcol] = col; 119 /* next three lines insert new entry into linked list */ 120 rowhit[mfm] = currentcol; 121 rowhit[currentcol] = fm; 122 fm = currentcol; 123 /* fm points to present position in list since we know the columns are sorted */ 124 } else { 125 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Detected invalid coloring"); 126 } 127 128 } 129 } 130 c->nrows[i] = nrows; 131 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 132 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 133 /* now store the linked list of rows into c->rows[i] */ 134 nrows = 0; 135 fm = rowhit[N]; 136 do { 137 c->rows[i][nrows] = fm; 138 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 139 fm = rowhit[fm]; 140 } while (fm < N); 141 } /* ---------------------------------------------------------------------------------------*/ 142 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 143 } 144 ierr = MatRestoreColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 145 146 ierr = PetscFree(rowhit);CHKERRQ(ierr); 147 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 148 149 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 150 /* 151 see the version for MPIAIJ 152 */ 153 ierr = VecCreateGhost(((PetscObject)mat)->comm,mat->rmap->n,PETSC_DETERMINE,0,PETSC_NULL,&c->vscale);CHKERRQ(ierr); 154 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 155 for (k=0; k<c->ncolors; k++) { 156 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 157 for (l=0; l<c->nrows[k]; l++) { 158 col = c->columnsforrow[k][l]; 159 c->vscaleforrow[k][l] = col; 160 } 161 } 162 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 163 PetscFunctionReturn(0); 164 } 165