1 2 #include <../src/mat/impls/aij/seq/aij.h> 3 #include <../src/mat/impls/baij/seq/baij.h> 4 5 /* 6 This routine is shared by SeqAIJ and SeqBAIJ matrices, 7 since it operators only on the nonzero structure of the elements or blocks. 8 */ 9 #undef __FUNCT__ 10 #define __FUNCT__ "MatFDColoringCreate_SeqXAIJ" 11 PetscErrorCode MatFDColoringCreate_SeqXAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 12 { 13 PetscErrorCode ierr; 14 PetscInt bs,nis=iscoloring->n,m=mat->rmap->n; 15 PetscBool isBAIJ; 16 17 PetscFunctionBegin; 18 /* set default brows and bcols for speedup inserting the dense matrix into sparse Jacobian */ 19 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 20 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr); 21 if (isBAIJ) { 22 c->brows = m; 23 c->bcols = 1; 24 } else { /* seqaij matrix */ 25 /* bcols is chosen s.t. dy-array takes 50% of memory space as mat */ 26 Mat_SeqAIJ *spA = (Mat_SeqAIJ*)mat->data; 27 PetscReal mem; 28 PetscInt nz,brows,bcols; 29 30 bs = 1; /* only bs=1 is supported for SeqAIJ matrix */ 31 32 nz = spA->nz; 33 mem = nz*(sizeof(PetscScalar) + sizeof(PetscInt)) + 3*m*sizeof(PetscInt); 34 bcols = (PetscInt)(0.5*mem /(m*sizeof(PetscScalar))); 35 brows = 1000/bcols; 36 if (bcols > nis) bcols = nis; 37 if (brows == 0 || brows > m) brows = m; 38 c->brows = brows; 39 c->bcols = bcols; 40 } 41 42 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 43 c->N = mat->cmap->N/bs; 44 c->m = mat->rmap->N/bs; 45 c->rstart = 0; 46 c->ncolors = nis; 47 c->ctype = IS_COLORING_GHOSTED; 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "MatFDColoringSetUpBlocked_AIJ_Private" 53 PetscErrorCode MatFDColoringSetUpBlocked_AIJ_Private(Mat mat,MatFDColoring c,PetscInt nz) 54 { 55 PetscErrorCode ierr; 56 PetscInt i,j,nrows,nbcols=0,brows=c->brows,bcols=c->bcols,mbs=c->m,*row_start,*nrows_new,nz_new,row_end,nis=c->ncolors; 57 PetscInt *color_start; 58 MatEntry *Jentry_new,*Jentry=c->matentry; 59 60 PetscFunctionBegin; 61 if (brows < 1 || brows > mbs) brows = mbs; 62 ierr = PetscMalloc((nis+1)*sizeof(PetscInt),&color_start);CHKERRQ(ierr); 63 color_start[0] = 0; 64 for (i=0; i<nis; i++) color_start[i+1] = c->nrows[i] + color_start[i]; 65 66 ierr = PetscMalloc(nz*sizeof(MatEntry),&Jentry_new);CHKERRQ(ierr); 67 ierr = PetscMalloc(bcols*sizeof(PetscInt),&row_start);CHKERRQ(ierr); 68 ierr = PetscMalloc(nis*sizeof(PetscInt),&nrows_new);CHKERRQ(ierr); 69 70 nz_new = 0; 71 for (i=0; i<nis; i+=bcols) { /* loop over colors */ 72 if (i + bcols > nis) bcols = nis - i; 73 74 row_end = brows; 75 if (row_end > mbs) row_end = mbs; 76 for (j=0; j<bcols; j++) row_start[j] = 0; 77 while (row_end <= mbs) { /* loop over block rows */ 78 for (j=0; j<bcols; j++) { /* loop over block columns */ 79 nrows = c->nrows[i+j]; 80 nz = color_start[i+j]; 81 while (row_start[j] < nrows) { 82 if (Jentry[nz].row >= row_end) { 83 color_start[i+j] = nz; 84 break; 85 } else { /* copy Jentry[nz] to Jentry_new[nz_new] */ 86 Jentry_new[nz_new].row = Jentry[nz].row + j*mbs; /* index in dy-array */ 87 Jentry_new[nz_new].col = Jentry[nz].col; 88 Jentry_new[nz_new].valaddr = Jentry[nz].valaddr; 89 nz_new++; nz++; row_start[j]++; 90 } 91 } 92 } 93 if (row_end == mbs) break; 94 row_end += brows; 95 if (row_end > mbs) row_end = mbs; 96 } 97 nrows_new[nbcols++] = nz_new; 98 } 99 for (i=nbcols-1; i>0; i--) nrows_new[i] -= nrows_new[i-1]; 100 ierr = PetscFree(c->nrows);CHKERRQ(ierr); 101 c->nrows = nrows_new; 102 103 ierr = PetscFree(Jentry);CHKERRQ(ierr); 104 c->matentry = Jentry_new; 105 ierr = PetscMalloc(c->bcols*mat->rmap->n*sizeof(PetscScalar),&c->dy);CHKERRQ(ierr); 106 ierr = PetscFree(row_start);CHKERRQ(ierr); 107 ierr = PetscFree(color_start);CHKERRQ(ierr); 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "MatFDColoringSetUp_SeqXAIJ" 113 PetscErrorCode MatFDColoringSetUp_SeqXAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 114 { 115 PetscErrorCode ierr; 116 PetscInt i,n,nrows,mbs=c->m,j,k,m,ncols,col,nis=iscoloring->n,*rowhit,bs,bs2,*spidx,nz; 117 const PetscInt *is,*row,*ci,*cj; 118 IS *isa; 119 PetscBool isBAIJ; 120 PetscScalar *A_val,**valaddrhit; 121 MatEntry *Jentry; 122 123 PetscFunctionBegin; 124 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 125 126 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 127 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr); 128 if (isBAIJ) { 129 Mat_SeqBAIJ *spA = (Mat_SeqBAIJ*)mat->data; 130 A_val = spA->a; 131 nz = spA->nz; 132 } else { 133 Mat_SeqAIJ *spA = (Mat_SeqAIJ*)mat->data; 134 A_val = spA->a; 135 nz = spA->nz; 136 bs = 1; /* only bs=1 is supported for SeqAIJ matrix */ 137 } 138 139 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 140 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 141 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 142 ierr = PetscLogObjectMemory((PetscObject)c,3*nis*sizeof(PetscInt));CHKERRQ(ierr); 143 144 ierr = PetscMalloc(nz*sizeof(MatEntry),&Jentry);CHKERRQ(ierr); 145 ierr = PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry));CHKERRQ(ierr); 146 c->matentry = Jentry; 147 148 if (isBAIJ) { 149 ierr = MatGetColumnIJ_SeqBAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 150 } else { 151 ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 152 } 153 154 ierr = PetscMalloc2(c->m,PetscInt,&rowhit,c->m,PetscScalar*,&valaddrhit);CHKERRQ(ierr); 155 ierr = PetscMemzero(rowhit,c->m*sizeof(PetscInt));CHKERRQ(ierr); 156 157 nz = 0; 158 for (i=0; i<nis; i++) { /* loop over colors */ 159 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 160 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 161 162 c->ncolumns[i] = n; 163 if (n) { 164 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 165 ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr); 166 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 167 } else { 168 c->columns[i] = 0; 169 } 170 171 /* fast, crude version requires O(N*N) work */ 172 bs2 = bs*bs; 173 nrows = 0; 174 for (j=0; j<n; j++) { /* loop over columns */ 175 col = is[j]; 176 row = cj + ci[col]; 177 m = ci[col+1] - ci[col]; 178 nrows += m; 179 for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */ 180 rowhit[*row] = col + 1; 181 valaddrhit[*row++] = &A_val[bs2*spidx[ci[col] + k]]; 182 } 183 } 184 c->nrows[i] = nrows; /* total num of rows for this color */ 185 186 for (j=0; j<mbs; j++) { /* loop over rows */ 187 if (rowhit[j]) { 188 Jentry[nz].row = j; /* local row index */ 189 Jentry[nz].col = rowhit[j] - 1; /* local column index */ 190 Jentry[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */ 191 nz++; 192 rowhit[j] = 0.0; /* zero rowhit for reuse */ 193 } 194 } 195 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 196 } 197 198 if (c->bcols > 1) { /* reorder Jentry for faster MatFDColoringApply() */ 199 ierr = MatFDColoringSetUpBlocked_AIJ_Private(mat,c,nz);CHKERRQ(ierr); 200 } 201 202 if (isBAIJ) { 203 ierr = MatRestoreColumnIJ_SeqBAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 204 ierr = PetscMalloc(bs*mat->rmap->n*sizeof(PetscScalar),&c->dy);CHKERRQ(ierr); 205 } else { 206 ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 207 } 208 ierr = PetscFree2(rowhit,valaddrhit);CHKERRQ(ierr); 209 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 210 211 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->rmap->n,PETSC_DETERMINE,0,NULL,&c->vscale);CHKERRQ(ierr); 212 PetscFunctionReturn(0); 213 } 214