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