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