1 /*$Id: fdaij.c,v 1.40 2001/06/21 21:16:21 bsmith Exp $*/ 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 __FUNCT__ 10 #define __FUNCT__ "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; 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 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 24 c->M = mat->M; /* set total rows, columns and local rows */ 25 c->N = mat->N; 26 c->m = mat->M; 27 c->rstart = 0; 28 29 c->ncolors = nis; 30 ierr = PetscMalloc(nis*sizeof(int),&c->ncolumns);CHKERRQ(ierr); 31 ierr = PetscMalloc(nis*sizeof(int*),&c->columns);CHKERRQ(ierr); 32 ierr = PetscMalloc(nis*sizeof(int),&c->nrows);CHKERRQ(ierr); 33 ierr = PetscMalloc(nis*sizeof(int*),&c->rows);CHKERRQ(ierr); 34 ierr = PetscMalloc(nis*sizeof(int*),&c->columnsforrow);CHKERRQ(ierr); 35 36 /* 37 Calls the _SeqAIJ() version of these routines to make sure it does not 38 get the reduced (by inodes) version of I and J 39 */ 40 ierr = MatGetColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 41 42 /* 43 Temporary option to allow for debugging/testing 44 */ 45 ierr = PetscOptionsHasName(PETSC_NULL,"-matfdcoloring_slow",&flg);CHKERRQ(ierr); 46 47 ierr = PetscMalloc((N+1)*sizeof(int),&rowhit);CHKERRQ(ierr); 48 ierr = PetscMalloc((N+1)*sizeof(int),&columnsforrow);CHKERRQ(ierr); 49 50 for (i=0; i<nis; i++) { 51 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 52 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 53 c->ncolumns[i] = n; 54 if (n) { 55 ierr = PetscMalloc(n*sizeof(int),&c->columns[i]);CHKERRQ(ierr); 56 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(int));CHKERRQ(ierr); 57 } else { 58 c->columns[i] = 0; 59 } 60 61 if (!flg) { /* ------------------------------------------------------------------------------*/ 62 /* fast, crude version requires O(N*N) work */ 63 ierr = PetscMemzero(rowhit,N*sizeof(int));CHKERRQ(ierr); 64 /* loop over columns*/ 65 for (j=0; j<n; j++) { 66 col = is[j]; 67 rows = cj + ci[col]; 68 m = ci[col+1] - ci[col]; 69 /* loop over columns marking them in rowhit */ 70 for (k=0; k<m; k++) { 71 rowhit[*rows++] = col + 1; 72 } 73 } 74 /* count the number of hits */ 75 nrows = 0; 76 for (j=0; j<N; j++) { 77 if (rowhit[j]) nrows++; 78 } 79 c->nrows[i] = nrows; 80 ierr = PetscMalloc((nrows+1)*sizeof(int),&c->rows[i]);CHKERRQ(ierr); 81 ierr = PetscMalloc((nrows+1)*sizeof(int),&c->columnsforrow[i]);CHKERRQ(ierr); 82 nrows = 0; 83 for (j=0; j<N; j++) { 84 if (rowhit[j]) { 85 c->rows[i][nrows] = j; 86 c->columnsforrow[i][nrows] = rowhit[j] - 1; 87 nrows++; 88 } 89 } 90 } else { /*-------------------------------------------------------------------------------*/ 91 /* slow version, using rowhit as a linked list */ 92 int currentcol,fm,mfm; 93 rowhit[N] = N; 94 nrows = 0; 95 /* loop over columns */ 96 for (j=0; j<n; j++) { 97 col = is[j]; 98 rows = cj + ci[col]; 99 m = ci[col+1] - ci[col]; 100 /* loop over columns marking them in rowhit */ 101 fm = N; /* fm points to first entry in linked list */ 102 for (k=0; k<m; k++) { 103 currentcol = *rows++; 104 /* is it already in the list? */ 105 do { 106 mfm = fm; 107 fm = rowhit[fm]; 108 } while (fm < currentcol); 109 /* not in list so add it */ 110 if (fm != currentcol) { 111 nrows++; 112 columnsforrow[currentcol] = col; 113 /* next three lines insert new entry into linked list */ 114 rowhit[mfm] = currentcol; 115 rowhit[currentcol] = fm; 116 fm = currentcol; 117 /* fm points to present position in list since we know the columns are sorted */ 118 } else { 119 SETERRQ(PETSC_ERR_PLIB,"Detected invalid coloring"); 120 } 121 122 } 123 } 124 c->nrows[i] = nrows; 125 ierr = PetscMalloc((nrows+1)*sizeof(int),&c->rows[i]);CHKERRQ(ierr); 126 ierr = PetscMalloc((nrows+1)*sizeof(int),&c->columnsforrow[i]);CHKERRQ(ierr); 127 /* now store the linked list of rows into c->rows[i] */ 128 nrows = 0; 129 fm = rowhit[N]; 130 do { 131 c->rows[i][nrows] = fm; 132 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 133 fm = rowhit[fm]; 134 } while (fm < N); 135 } /* ---------------------------------------------------------------------------------------*/ 136 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 137 } 138 ierr = MatRestoreColumnIJ_SeqAIJ(mat,0,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 139 140 ierr = PetscFree(rowhit);CHKERRQ(ierr); 141 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 142 143 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 144 /* 145 see the version for MPIAIJ 146 */ 147 ierr = VecCreateGhost(mat->comm,mat->m,PETSC_DETERMINE,0,PETSC_NULL,&c->vscale);CHKERRQ(ierr) 148 ierr = PetscMalloc(c->ncolors*sizeof(int*),&c->vscaleforrow);CHKERRQ(ierr); 149 for (k=0; k<c->ncolors; k++) { 150 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(int),&c->vscaleforrow[k]);CHKERRQ(ierr); 151 for (l=0; l<c->nrows[k]; l++) { 152 col = c->columnsforrow[k][l]; 153 c->vscaleforrow[k][l] = col; 154 } 155 } 156 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 /* 161 Makes a longer coloring[] array and calls the usual code with that 162 */ 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatColoringPatch_SeqAIJ_Inode" 165 int MatColoringPatch_SeqAIJ_Inode(Mat mat,int nin,int ncolors,int *coloring,ISColoring *iscoloring) 166 { 167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 168 int n = mat->n,ierr,m = a->inode.node_count,j,*ns = a->inode.size,row; 169 int *colorused,i,*newcolor; 170 171 PetscFunctionBegin; 172 ierr = PetscMalloc((n+1)*sizeof(int),&newcolor);CHKERRQ(ierr); 173 /* loop over inodes, marking a color for each column*/ 174 row = 0; 175 for (i=0; i<m; i++){ 176 for (j=0; j<ns[i]; j++) { 177 newcolor[row++] = coloring[i] + j*ncolors; 178 } 179 } 180 181 /* eliminate unneeded colors */ 182 ierr = PetscMalloc(5*ncolors*sizeof(int),&colorused);CHKERRQ(ierr); 183 ierr = PetscMemzero(colorused,5*ncolors*sizeof(int));CHKERRQ(ierr); 184 for (i=0; i<n; i++) { 185 colorused[newcolor[i]] = 1; 186 } 187 188 for (i=1; i<5*ncolors; i++) { 189 colorused[i] += colorused[i-1]; 190 } 191 ncolors = colorused[5*ncolors-1]; 192 for (i=0; i<n; i++) { 193 newcolor[i] = colorused[newcolor[i]]; 194 } 195 ierr = PetscFree(colorused);CHKERRQ(ierr); 196 ierr = PetscFree(coloring);CHKERRQ(ierr); 197 ierr = ISColoringCreate(mat->comm,n,newcolor,iscoloring);CHKERRQ(ierr); 198 199 PetscFunctionReturn(0); 200 } 201 202 203 204 205 206 207