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