1 2 3 #ifndef lint 4 static char vcid[] = "$Id: fdaij.c,v 1.1 1996/10/09 18:09:59 bsmith Exp bsmith $"; 5 #endif 6 7 #include "src/mat/impls/aij/seq/aij.h" 8 #include "src/vec/vecimpl.h" 9 #include "petsc.h" 10 11 12 int MatFDColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 13 { 14 int i,*is,n,nrows,N = mat->N,j,k,m,*rows,ierr,*ci,*cj,ncols,col,flg; 15 int nis = iscoloring->n; 16 IS *isa = iscoloring->is; 17 PetscTruth done; 18 19 c->ncolors = nis; 20 c->ncolumns = (int *) PetscMalloc( nis*sizeof(int) ); CHKPTRQ(c->ncolumns); 21 c->columns = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->columns); 22 c->nrows = (int *) PetscMalloc( nis*sizeof(int) ); CHKPTRQ(c->nrows); 23 c->rows = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->rows); 24 c->columnsforrow = (int **) PetscMalloc( nis*sizeof(int *)); CHKPTRQ(c->columnsforrow); 25 ierr = MatGetColumnIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done); CHKERRQ(ierr); 26 /* 27 Temporary option to allow for debugging/testing 28 */ 29 ierr = OptionsHasName(0,"-matfdcoloring_slow",&flg); 30 for ( i=0; i<nis; i++ ) { 31 ierr = ISGetSize(isa[i],&n); CHKERRQ(ierr); 32 ierr = ISGetIndices(isa[i],&is); CHKERRQ(ierr); 33 c->ncolumns[i] = n; 34 if (n) { 35 c->columns[i] = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(c->columns[i]); 36 PetscMemcpy(c->columns[i],is,n*sizeof(int)); 37 } else { 38 c->columns[i] = 0; 39 } 40 41 if (flg) { /* --------------------------------------------------------------------------------*/ 42 /* crude version requires O(N*N) work */ 43 int *rowhit = (int *) PetscMalloc( N*sizeof(int) ); CHKPTRQ(rowhit); 44 PetscMemzero(rowhit,N*sizeof(int)); 45 /* loop over columns*/ 46 for ( j=0; j<n; j++ ) { 47 col = is[j]; 48 rows = cj + ci[col]; 49 m = ci[col+1] - ci[col]; 50 /* loop over columns marking them in rowhit */ 51 for ( k=0; k<m; k++ ) { 52 rowhit[*rows++] = col + 1; 53 } 54 } 55 /* count the number of hits */ 56 nrows = 0; 57 for ( j=0; j<N; j++ ) { 58 if (rowhit[j]) nrows++; 59 } 60 c->nrows[i] = nrows; 61 c->rows[i] = (int *) PetscMalloc(nrows*sizeof(int)); CHKPTRQ(c->rows[i]); 62 c->columnsforrow[i] = (int *) PetscMalloc(nrows*sizeof(int)); CHKPTRQ(c->columnsforrow[i]); 63 nrows = 0; 64 for ( j=0; j<N; j++ ) { 65 if (rowhit[j]) { 66 c->rows[i][nrows] = j; 67 c->columnsforrow[i][nrows] = rowhit[j] - 1; 68 nrows++; 69 } 70 } 71 PetscFree(rowhit); 72 } else { /*-------------------------------------------------------------------------------*/ 73 /* efficient version, using rowhit as a linked list */ 74 int currentcol,fm,mfm,*rowhit,*columnsforrow; 75 rowhit = (int *) PetscMalloc( (N+1)*sizeof(int) ); CHKPTRQ(rowhit); 76 columnsforrow = (int *) PetscMalloc( (N+1)*sizeof(int) ); CHKPTRQ(columnsforrow); 77 rowhit[N] = N; 78 nrows = 0; 79 /* loop over columns */ 80 for ( j=0; j<n; j++ ) { 81 col = is[j]; 82 rows = cj + ci[col]; 83 m = ci[col+1] - ci[col]; 84 /* loop over columns marking them in rowhit */ 85 fm = N; /* fm points to first entry in linked list */ 86 for ( k=0; k<m; k++ ) { 87 currentcol = *rows++; 88 /* is it already in the list? */ 89 do { 90 mfm = fm; 91 fm = rowhit[fm]; 92 } while (fm < currentcol); 93 /* not in list so add it */ 94 if (fm != currentcol) { 95 nrows++; 96 columnsforrow[currentcol] = col; 97 /* next three lines insert new entry into linked list */ 98 rowhit[mfm] = currentcol; 99 rowhit[currentcol] = fm; 100 fm = currentcol; 101 /* fm points to present position in list since we know the columns are sorted */ 102 } else { 103 SETERRQ(1,"MatFDColoringCreate_SeqAIJ:Invalid coloring"); 104 } 105 106 } 107 } 108 c->nrows[i] = nrows; 109 c->rows[i] = (int *) PetscMalloc(nrows*sizeof(int));CHKPTRQ(c->rows[i]); 110 c->columnsforrow[i] = (int *) PetscMalloc(nrows*sizeof(int));CHKPTRQ(c->columnsforrow[i]); 111 /* now store the linked list of rows into c->rows[i] */ 112 nrows = 0; 113 fm = rowhit[N]; 114 do { 115 c->rows[i][nrows] = fm; 116 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 117 fm = rowhit[fm]; 118 } while (fm < N); 119 PetscFree(rowhit); 120 PetscFree(columnsforrow); 121 } /* ---------------------------------------------------------------------------------------*/ 122 ierr = ISRestoreIndices(isa[i],&is); CHKERRQ(ierr); 123 } 124 ierr = MatRestoreColumnIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done); CHKERRQ(ierr); 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; 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 (*iscoloring)->comm = MPI_COMM_SELF; 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