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