1 2 #ifdef PETSC_RCS_HEADER 3 static char vcid[] = "$Id: fdaij.c,v 1.20 1999/05/04 20:31:42 balay Exp balay $"; 4 #endif 5 6 #include "src/mat/impls/aij/seq/aij.h" 7 #include "src/vec/vecimpl.h" 8 #include "petsc.h" 9 10 extern int MatGetColumnIJ_SeqAIJ(Mat,int,PetscTruth,int*,int**,int**,PetscTruth*); 11 extern int MatRestoreColumnIJ_SeqAIJ(Mat,int,PetscTruth,int*,int**,int**,PetscTruth*); 12 13 #undef __FUNC__ 14 #define __FUNC__ "MatFDColoringCreate_SeqAIJ" 15 int MatFDColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 16 { 17 int i,*is,n,nrows,N = mat->N,j,k,m,*rows,ierr,*ci,*cj,ncols,col,flg; 18 int nis = iscoloring->n,*rowhit,*columnsforrow; 19 IS *isa = iscoloring->is; 20 PetscTruth done; 21 22 PetscFunctionBegin; 23 if (!mat->assembled) { 24 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,1,"Matrix must be assembled by calls to MatAssemblyBegin/End();"); 25 } 26 27 c->M = mat->M; /* set total rows, columns and local rows */ 28 c->N = mat->N; 29 c->m = mat->M; 30 c->rstart = 0; 31 32 c->ncolors = nis; 33 c->ncolumns = (int *) PetscMalloc( nis*sizeof(int) );CHKPTRQ(c->ncolumns); 34 c->columns = (int **) PetscMalloc( nis*sizeof(int *));CHKPTRQ(c->columns); 35 c->nrows = (int *) PetscMalloc( nis*sizeof(int) );CHKPTRQ(c->nrows); 36 c->rows = (int **) PetscMalloc( nis*sizeof(int *));CHKPTRQ(c->rows); 37 c->columnsforrow = (int **) PetscMalloc( nis*sizeof(int *));CHKPTRQ(c->columnsforrow); 38 39 /* 40 Calls the _SeqAIJ() version of these routines to make sure it does not 41 get the reduced (by inodes) version of I and J 42 */ 43 ierr = MatGetColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 44 45 /* 46 Temporary option to allow for debugging/testing 47 */ 48 ierr = OptionsHasName(0,"-matfdcoloring_slow",&flg);CHKERRQ(ierr); 49 50 rowhit = (int *) PetscMalloc( (N+1)*sizeof(int) );CHKPTRQ(rowhit); 51 columnsforrow = (int *) PetscMalloc( (N+1)*sizeof(int) );CHKPTRQ(columnsforrow); 52 53 for ( i=0; i<nis; i++ ) { 54 ierr = ISGetSize(isa[i],&n);CHKERRQ(ierr); 55 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 56 c->ncolumns[i] = n; 57 if (n) { 58 c->columns[i] = (int *) PetscMalloc( n*sizeof(int) );CHKPTRQ(c->columns[i]); 59 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(int));CHKERRQ(ierr); 60 } else { 61 c->columns[i] = 0; 62 } 63 64 if (flg) { /* ------------------------------------------------------------------------------*/ 65 /* crude version requires O(N*N) work */ 66 ierr = PetscMemzero(rowhit,N*sizeof(int));CHKERRQ(ierr); 67 /* loop over columns*/ 68 for ( j=0; j<n; j++ ) { 69 col = is[j]; 70 rows = cj + ci[col]; 71 m = ci[col+1] - ci[col]; 72 /* loop over columns marking them in rowhit */ 73 for ( k=0; k<m; k++ ) { 74 rowhit[*rows++] = col + 1; 75 } 76 } 77 /* count the number of hits */ 78 nrows = 0; 79 for ( j=0; j<N; j++ ) { 80 if (rowhit[j]) nrows++; 81 } 82 c->nrows[i] = nrows; 83 c->rows[i] = (int *) PetscMalloc(nrows*sizeof(int));CHKPTRQ(c->rows[i]); 84 c->columnsforrow[i] = (int *) PetscMalloc(nrows*sizeof(int));CHKPTRQ(c->columnsforrow[i]); 85 nrows = 0; 86 for ( j=0; j<N; j++ ) { 87 if (rowhit[j]) { 88 c->rows[i][nrows] = j; 89 c->columnsforrow[i][nrows] = rowhit[j] - 1; 90 nrows++; 91 } 92 } 93 } else { /*-------------------------------------------------------------------------------*/ 94 /* efficient version, using rowhit as a linked list */ 95 int currentcol,fm,mfm; 96 rowhit[N] = N; 97 nrows = 0; 98 /* loop over columns */ 99 for ( j=0; j<n; j++ ) { 100 col = is[j]; 101 rows = cj + ci[col]; 102 m = ci[col+1] - ci[col]; 103 /* loop over columns marking them in rowhit */ 104 fm = N; /* fm points to first entry in linked list */ 105 for ( k=0; k<m; k++ ) { 106 currentcol = *rows++; 107 /* is it already in the list? */ 108 do { 109 mfm = fm; 110 fm = rowhit[fm]; 111 } while (fm < currentcol); 112 /* not in list so add it */ 113 if (fm != currentcol) { 114 nrows++; 115 columnsforrow[currentcol] = col; 116 /* next three lines insert new entry into linked list */ 117 rowhit[mfm] = currentcol; 118 rowhit[currentcol] = fm; 119 fm = currentcol; 120 /* fm points to present position in list since we know the columns are sorted */ 121 } else { 122 SETERRQ(PETSC_ERR_PLIB,0,"Detected invalid coloring"); 123 } 124 125 } 126 } 127 c->nrows[i] = nrows; 128 c->rows[i] = (int *)PetscMalloc((nrows+1)*sizeof(int));CHKPTRQ(c->rows[i]); 129 c->columnsforrow[i] = (int *)PetscMalloc((nrows+1)*sizeof(int));CHKPTRQ(c->columnsforrow[i]); 130 /* now store the linked list of rows into c->rows[i] */ 131 nrows = 0; 132 fm = rowhit[N]; 133 do { 134 c->rows[i][nrows] = fm; 135 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 136 fm = rowhit[fm]; 137 } while (fm < N); 138 } /* ---------------------------------------------------------------------------------------*/ 139 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 140 } 141 ierr = MatRestoreColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 142 143 PetscFree(rowhit); 144 PetscFree(columnsforrow); 145 146 c->scale = (Scalar *) PetscMalloc( 2*N*sizeof(Scalar) );CHKPTRQ(c->scale); 147 c->wscale = c->scale + N; 148 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNC__ 153 #define __FUNC__ "MatColoringPatch_SeqAIJ" 154 int MatColoringPatch_SeqAIJ(Mat mat,int ncolors,int *coloring,ISColoring *iscoloring) 155 { 156 Mat_SeqAIJ *a = (Mat_SeqAIJ *) mat->data; 157 int n = a->n,*sizes,i,**ii,ierr,tag; 158 IS *is; 159 160 PetscFunctionBegin; 161 /* construct the index sets from the coloring array */ 162 sizes = (int *) PetscMalloc( ncolors*sizeof(int) );CHKPTRQ(sizes); 163 ierr = PetscMemzero(sizes,ncolors*sizeof(int));CHKERRQ(ierr); 164 for ( i=0; i<n; i++ ) { 165 sizes[coloring[i]-1]++; 166 } 167 ii = (int **) PetscMalloc( ncolors*sizeof(int*) );CHKPTRQ(ii); 168 ii[0] = (int *) PetscMalloc( n*sizeof(int) );CHKPTRQ(ii[0]); 169 for ( i=1; i<ncolors; i++ ) { 170 ii[i] = ii[i-1] + sizes[i-1]; 171 } 172 ierr = PetscMemzero(sizes,ncolors*sizeof(int));CHKERRQ(ierr); 173 for ( i=0; i<n; i++ ) { 174 ii[coloring[i]-1][sizes[coloring[i]-1]++] = i; 175 } 176 is = (IS *) PetscMalloc( ncolors*sizeof(IS) );CHKPTRQ(is); 177 for ( i=0; i<ncolors; i++ ) { 178 ierr = ISCreateGeneral(PETSC_COMM_SELF,sizes[i],ii[i],is+i);CHKERRQ(ierr); 179 } 180 181 *iscoloring = (ISColoring) PetscMalloc(sizeof(struct _p_ISColoring));CHKPTRQ(*iscoloring); 182 (*iscoloring)->n = ncolors; 183 (*iscoloring)->is = is; 184 ierr = PetscCommDuplicate_Private(mat->comm,&(*iscoloring)->comm,&tag);CHKERRQ(ierr); 185 PetscFree(sizes); 186 PetscFree(ii[0]); 187 PetscFree(ii); 188 PetscFunctionReturn(0); 189 } 190 191 /* 192 Makes a longer coloring[] array and calls the usual code with that 193 */ 194 #undef __FUNC__ 195 #define __FUNC__ "MatColoringPatch_SeqAIJ_Inode" 196 int MatColoringPatch_SeqAIJ_Inode(Mat mat,int ncolors,int *coloring,ISColoring *iscoloring) 197 { 198 Mat_SeqAIJ *a = (Mat_SeqAIJ *) mat->data; 199 int n = a->n,ierr, m = a->inode.node_count,j,*ns = a->inode.size,row; 200 int *colorused,i,*newcolor; 201 202 PetscFunctionBegin; 203 newcolor = (int *) PetscMalloc((n+1)*sizeof(int));CHKPTRQ(newcolor); 204 205 /* loop over inodes, marking a color for each column*/ 206 row = 0; 207 for ( i=0; i<m; i++){ 208 for ( j=0; j<ns[i]; j++) { 209 newcolor[row++] = coloring[i] + j*ncolors; 210 } 211 } 212 213 /* eliminate unneeded colors */ 214 colorused = (int *) PetscMalloc( 5*ncolors*sizeof(int) );CHKPTRQ(colorused); 215 ierr = PetscMemzero(colorused,5*ncolors*sizeof(int));CHKERRQ(ierr); 216 for ( i=0; i<n; i++ ) { 217 colorused[newcolor[i]-1] = 1; 218 } 219 220 for ( i=1; i<5*ncolors; i++ ) { 221 colorused[i] += colorused[i-1]; 222 } 223 ncolors = colorused[5*ncolors-1]; 224 for ( i=0; i<n; i++ ) { 225 newcolor[i] = colorused[newcolor[i]-1]; 226 } 227 PetscFree(colorused); 228 229 ierr = MatColoringPatch_SeqAIJ(mat,ncolors,newcolor,iscoloring);CHKERRQ(ierr); 230 PetscFree(newcolor); 231 232 PetscFunctionReturn(0); 233 } 234 235 236 237 238 239 240