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