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