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