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