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 i,n,nrows,N,j,k,m,ncols,col; 15 const PetscInt *is,*row,*ci,*cj; 16 PetscInt nis=iscoloring->n,*rowhit,bs,bs2,*spidx,nz; 17 IS *isa; 18 PetscBool isBAIJ; 19 Mat_SeqAIJ *spA = (Mat_SeqAIJ*)mat->data; 20 PetscScalar *A_val=spA->a; 21 PetscScalar **valaddrhit; 22 MatEntry *Jentry; 23 24 PetscFunctionBegin; 25 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 26 27 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 28 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 29 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr); 30 if (!isBAIJ) { 31 bs = 1; /* only bs=1 is supported for non SEQBAIJ matrix */ 32 } 33 N = mat->cmap->N/bs; 34 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 35 c->N = mat->cmap->N/bs; 36 c->m = mat->rmap->N/bs; 37 c->rstart = 0; 38 39 c->ncolors = nis; 40 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 41 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 42 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 43 ierr = PetscLogObjectMemory((PetscObject)c,3*nis*sizeof(PetscInt));CHKERRQ(ierr); 44 45 ierr = PetscMalloc(spA->nz*sizeof(MatEntry),&Jentry);CHKERRQ(ierr); 46 ierr = PetscLogObjectMemory((PetscObject)c,spA->nz*sizeof(MatEntry));CHKERRQ(ierr); 47 c->matentry = Jentry; 48 49 if (isBAIJ) { 50 ierr = MatGetColumnIJ_SeqBAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 51 } else { 52 ierr = MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 53 } 54 55 ierr = PetscMalloc2(c->m,PetscInt,&rowhit,c->m,PetscScalar*,&valaddrhit);CHKERRQ(ierr); 56 ierr = PetscMemzero(rowhit,c->m*sizeof(PetscInt));CHKERRQ(ierr); 57 58 nz = 0; 59 for (i=0; i<nis; i++) { /* loop over colors */ 60 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 61 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 62 63 c->ncolumns[i] = n; 64 if (n) { 65 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 66 ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr); 67 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 68 } else { 69 c->columns[i] = 0; 70 } 71 72 /* fast, crude version requires O(N*N) work */ 73 bs2 = bs*bs; 74 nrows = 0; 75 for (j=0; j<n; j++) { /* loop over columns */ 76 col = is[j]; 77 row = cj + ci[col]; 78 m = ci[col+1] - ci[col]; 79 nrows += m; 80 for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */ 81 rowhit[*row] = col + 1; 82 valaddrhit[*row++] = &A_val[bs2*spidx[ci[col] + k]]; 83 } 84 } 85 c->nrows[i] = nrows; /* total num of rows for this color */ 86 87 for (j=0; j<N; j++) { /* loop over rows */ 88 if (rowhit[j]) { 89 Jentry[nz].row = j; /* local row index */ 90 Jentry[nz].col = rowhit[j] - 1; /* local column index */ 91 Jentry[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */ 92 nz++; 93 rowhit[j] = 0.0; /* zero rowhit for reuse */ 94 } 95 } 96 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 97 } 98 if (nz != spA->nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"nz %d != mat->nz %d\n",nz,spA->nz); 99 100 if (isBAIJ) { 101 ierr = MatRestoreColumnIJ_SeqBAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 102 } else { 103 ierr = MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 104 } 105 ierr = PetscFree2(rowhit,valaddrhit);CHKERRQ(ierr); 106 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 107 108 c->ctype = IS_COLORING_GHOSTED; 109 if (isBAIJ) { 110 ierr = PetscMalloc(bs*mat->cmap->N*sizeof(PetscScalar),&c->dy);CHKERRQ(ierr); 111 } 112 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->rmap->n,PETSC_DETERMINE,0,NULL,&c->vscale);CHKERRQ(ierr); 113 PetscFunctionReturn(0); 114 } 115