1 2 #ifdef PETSC_RCS_HEADER 3 static char vcid[] = "$Id: fdaij.c,v 1.12 1997/07/09 20:53:48 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 #undef __FUNC__ 14 #define __FUNC__ "MatFDColoringCreate_SeqAIJ" /* ADIC Ignore */ 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->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 c->ncolumns = (int *) PetscMalloc( nis*sizeof(int) ); CHKPTRQ(c->ncolumns); 29 c->columns = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->columns); 30 c->nrows = (int *) PetscMalloc( nis*sizeof(int) ); CHKPTRQ(c->nrows); 31 c->rows = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->rows); 32 c->columnsforrow = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->columnsforrow); 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 = OptionsHasName(0,"-matfdcoloring_slow",&flg); 44 45 rowhit = (int *) PetscMalloc( (N+1)*sizeof(int) ); CHKPTRQ(rowhit); 46 columnsforrow = (int *) PetscMalloc( (N+1)*sizeof(int) ); CHKPTRQ(columnsforrow); 47 48 for ( i=0; i<nis; i++ ) { 49 ierr = ISGetSize(isa[i],&n); CHKERRQ(ierr); 50 ierr = ISGetIndices(isa[i],&is); CHKERRQ(ierr); 51 c->ncolumns[i] = n; 52 if (n) { 53 c->columns[i] = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(c->columns[i]); 54 PetscMemcpy(c->columns[i],is,n*sizeof(int)); 55 } else { 56 c->columns[i] = 0; 57 } 58 59 if (flg) { /* ------------------------------------------------------------------------------*/ 60 /* crude version requires O(N*N) work */ 61 PetscMemzero(rowhit,N*sizeof(int)); 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 c->rows[i] = (int *) PetscMalloc(nrows*sizeof(int)); CHKPTRQ(c->rows[i]); 79 c->columnsforrow[i] = (int *) PetscMalloc(nrows*sizeof(int)); CHKPTRQ(c->columnsforrow[i]); 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 /* efficient 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(1,0,"Invalid coloring"); 118 } 119 120 } 121 } 122 c->nrows[i] = nrows; 123 c->rows[i] = (int *)PetscMalloc((nrows+1)*sizeof(int));CHKPTRQ(c->rows[i]); 124 c->columnsforrow[i] = (int *)PetscMalloc((nrows+1)*sizeof(int));CHKPTRQ(c->columnsforrow[i]); 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 PetscFree(rowhit); 139 PetscFree(columnsforrow); 140 141 c->scale = (Scalar *) PetscMalloc( 2*N*sizeof(Scalar) ); CHKPTRQ(c->scale); 142 c->wscale = c->scale + N; 143 144 return 0; 145 } 146 147 #undef __FUNC__ 148 #define __FUNC__ "MatColoringPatch_SeqAIJ" /* ADIC Ignore */ 149 int MatColoringPatch_SeqAIJ(Mat mat,int ncolors,int *coloring,ISColoring *iscoloring) 150 { 151 Mat_SeqAIJ *a = (Mat_SeqAIJ *) mat->data; 152 int n = a->n,*sizes,i,**ii,ierr,tag; 153 IS *is; 154 155 /* construct the index sets from the coloring array */ 156 sizes = (int *) PetscMalloc( ncolors*sizeof(int) ); CHKPTRQ(sizes); 157 PetscMemzero(sizes,ncolors*sizeof(int)); 158 for ( i=0; i<n; i++ ) { 159 sizes[coloring[i]-1]++; 160 } 161 ii = (int **) PetscMalloc( ncolors*sizeof(int*) ); CHKPTRQ(ii); 162 ii[0] = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(ii[0]); 163 for ( i=1; i<ncolors; i++ ) { 164 ii[i] = ii[i-1] + sizes[i-1]; 165 } 166 PetscMemzero(sizes,ncolors*sizeof(int)); 167 for ( i=0; i<n; i++ ) { 168 ii[coloring[i]-1][sizes[coloring[i]-1]++] = i; 169 } 170 is = (IS *) PetscMalloc( ncolors*sizeof(IS) ); CHKPTRQ(is); 171 for ( i=0; i<ncolors; i++ ) { 172 ierr = ISCreateGeneral(PETSC_COMM_SELF,sizes[i],ii[i],is+i); CHKERRQ(ierr); 173 } 174 175 *iscoloring = (ISColoring) PetscMalloc(sizeof(struct _p_ISColoring));CHKPTRQ(*iscoloring); 176 (*iscoloring)->n = ncolors; 177 (*iscoloring)->is = is; 178 PetscCommDup_Private(mat->comm,&(*iscoloring)->comm,&tag); 179 PetscFree(sizes); 180 PetscFree(ii[0]); 181 PetscFree(ii); 182 return 0; 183 } 184 185 /* 186 Makes a longer coloring[] array and calls the usual code with that 187 */ 188 #undef __FUNC__ 189 #define __FUNC__ "MatColoringPatch_SeqAIJ_Inode" /* ADIC Ignore */ 190 int MatColoringPatch_SeqAIJ_Inode(Mat mat,int ncolors,int *coloring,ISColoring *iscoloring) 191 { 192 Mat_SeqAIJ *a = (Mat_SeqAIJ *) mat->data; 193 int n = a->n,ierr, m = a->inode.node_count,j,*ns = a->inode.size,row; 194 int *colorused,i,*newcolor; 195 196 newcolor = (int *) PetscMalloc((n+1)*sizeof(int)); CHKPTRQ(newcolor); 197 198 /* loop over inodes, marking a color for each column*/ 199 row = 0; 200 for ( i=0; i<m; i++){ 201 for ( j=0; j<ns[i]; j++) { 202 newcolor[row++] = coloring[i] + j*ncolors; 203 } 204 } 205 206 /* eliminate unneeded colors */ 207 colorused = (int *) PetscMalloc( 5*ncolors*sizeof(int) ); CHKPTRQ(colorused); 208 PetscMemzero(colorused,5*ncolors*sizeof(int)); 209 for ( i=0; i<n; i++ ) { 210 colorused[newcolor[i]-1] = 1; 211 } 212 213 for ( i=1; i<5*ncolors; i++ ) { 214 colorused[i] += colorused[i-1]; 215 } 216 ncolors = colorused[5*ncolors-1]; 217 for ( i=0; i<n; i++ ) { 218 newcolor[i] = colorused[newcolor[i]-1]; 219 } 220 PetscFree(colorused); 221 222 ierr = MatColoringPatch_SeqAIJ(mat,ncolors,newcolor,iscoloring); CHKERRQ(ierr); 223 PetscFree(newcolor); 224 225 return 0; 226 } 227 228