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