1 #ifndef lint 2 static char vcid[] = "$Id: mmbaij.c,v 1.2 1996/06/18 16:10:55 balay Exp balay $"; 3 #endif 4 5 6 /* 7 Support for the parallel BAIJ matrix vector multiply 8 */ 9 #include "mpibaij.h" 10 #include "src/vec/vecimpl.h" 11 #include "../seq/baij.h" 12 13 int MatSetUpMultiply_MPIBAIJ(Mat mat) 14 { 15 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *) mat->data; 16 Mat_SeqBAIJ *B = (Mat_SeqBAIJ *) (baij->B->data); 17 int Nbs = baij->Nbs,i,j,*indices,*aj = B->j,ierr,ec = 0,*garray; 18 int col,bs = baij->bs,*tmp; 19 IS from,to; 20 Vec gvec; 21 22 /* For the first stab we make an array as long as the number of columns */ 23 /* mark those columns that are in baij->B */ 24 indices = (int *) PetscMalloc( Nbs*sizeof(int) ); CHKPTRQ(indices); 25 PetscMemzero(indices,Nbs*sizeof(int)); 26 for ( i=0; i<B->mbs; i++ ) { 27 for ( j=0; j<B->ilen[i]; j++ ) { 28 if (!indices[aj[B->i[i] + j]]) ec++; 29 indices[aj[B->i[i] + j] ] = 1; 30 } 31 } 32 33 /* form array of columns we need */ 34 garray = (int *) PetscMalloc( (ec+1)*sizeof(int) ); CHKPTRQ(garray); 35 tmp = (int *) PetscMalloc( (ec*bs+1)*sizeof(int) ); CHKPTRQ(tmp) 36 ec = 0; 37 for ( i=0; i<Nbs; i++ ) { 38 if (indices[i]) garray[ec++] = i; 39 } 40 41 /* make indices now point into garray */ 42 for ( i=0; i<ec; i++ ) { 43 indices[garray[i]] = i; 44 } 45 46 /* compact out the extra columns in B */ 47 for ( i=0; i<B->mbs; i++ ) { 48 for ( j=0; j<B->ilen[i]; j++ ) { 49 aj[B->i[i] + j] = indices[aj[B->i[i] + j]]; 50 } 51 } 52 B->nbs = ec; 53 B->n = ec*B->bs; 54 PetscFree(indices); 55 56 for ( i=0,col=0; i<ec; i++ ) { 57 for ( j=0; j<bs; j++,col++) tmp[col] = garray[i]*bs+j; 58 } 59 /* create local vector that is used to scatter into */ 60 ierr = VecCreateSeq(MPI_COMM_SELF,ec*bs,&baij->lvec); CHKERRQ(ierr); 61 62 /* create two temporary Index sets for build scatter gather */ 63 ierr = ISCreateSeq(MPI_COMM_SELF,ec*bs,tmp,&from); CHKERRQ(ierr); 64 ierr = ISCreateStrideSeq(MPI_COMM_SELF,ec*bs,0,1,&to); CHKERRQ(ierr); 65 66 /* create temporary global vector to generate scatter context */ 67 /* this is inefficient, but otherwise we must do either 68 1) save garray until the first actual scatter when the vector is known or 69 2) have another way of generating a scatter context without a vector.*/ 70 ierr = VecCreateMPI(mat->comm,baij->n,baij->N,&gvec); CHKERRQ(ierr); 71 72 /* gnerate the scatter context */ 73 ierr = VecScatterCreate(gvec,from,baij->lvec,to,&baij->Mvctx); CHKERRQ(ierr); 74 PLogObjectParent(mat,baij->Mvctx); 75 PLogObjectParent(mat,baij->lvec); 76 PLogObjectParent(mat,from); 77 PLogObjectParent(mat,to); 78 baij->garray = garray; 79 PLogObjectMemory(mat,(ec+1)*sizeof(int)); 80 ierr = ISDestroy(from); CHKERRQ(ierr); 81 ierr = ISDestroy(to); CHKERRQ(ierr); 82 ierr = VecDestroy(gvec); 83 PetscFree(tmp); 84 return 0; 85 } 86 87 88 /* 89 Takes the local part of an already assembled MPIBAIJ matrix 90 and disassembles it. This is to allow new nonzeros into the matrix 91 that require more communication in the matrix vector multiply. 92 Thus certain data-structures must be rebuilt. 93 94 Kind of slow! But that's what application programmers get when 95 they are sloppy. 96 */ 97 int DisAssemble_MPIBAIJ(Mat A) 98 { 99 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *) A->data; 100 Mat B = baij->B,Bnew; 101 Mat_SeqBAIJ *Bbaij = (Mat_SeqBAIJ*)B->data; 102 int ierr,i,j,mbs=Bbaij->mbs,n = baij->N,col,*garray=baij->garray; 103 int k,bs=baij->bs,bs2=baij->bs2,*rvals,*nz,ec,m=Bbaij->m; 104 Scalar *a=Bbaij->a; 105 106 /* free stuff related to matrix-vec multiply */ 107 ierr = VecGetSize(baij->lvec,&ec); /* needed for PLogObjectMemory below */ 108 ierr = VecDestroy(baij->lvec); CHKERRQ(ierr); baij->lvec = 0; 109 ierr = VecScatterDestroy(baij->Mvctx); CHKERRQ(ierr); baij->Mvctx = 0; 110 if (baij->colmap) { 111 PetscFree(baij->colmap); baij->colmap = 0; 112 PLogObjectMemory(A,-Bbaij->nbs*sizeof(int)); 113 } 114 115 /* make sure that B is assembled so we can access its values */ 116 ierr = MatAssemblyBegin(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 117 MatAssemblyEnd(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 118 119 /* invent new B and copy stuff over */ 120 nz = (int *) PetscMalloc( mbs*sizeof(int) ); CHKPTRQ(nz); 121 for ( i=0; i<mbs; i++ ) { 122 nz[i] = Bbaij->i[i+1]-Bbaij->i[i]; 123 } 124 ierr = MatCreateSeqBAIJ(MPI_COMM_SELF,baij->bs,m,n,0,nz,&Bnew); CHKERRQ(ierr); 125 PetscFree(nz); 126 127 rvals = (int *) PetscMalloc(bs*sizeof(int)); CHKPTRQ(rvals); 128 for ( i=0; i<mbs; i++ ) { 129 rvals[0] = bs*i; 130 for ( j=1; j<bs; j++ ) { rvals[j] = rvals[j-1] + 1; } 131 for ( j=Bbaij->i[i]; j<Bbaij->i[i+1]; j++ ) { 132 col = garray[Bbaij->j[j]]*bs; 133 for (k=0; k<bs; k++ ) { 134 ierr = MatSetValues(Bnew,bs,rvals,1,&col,a+j*bs2,INSERT_VALUES);CHKERRQ(ierr); 135 col++; 136 } 137 } 138 } 139 PetscFree(baij->garray); baij->garray = 0; 140 PetscFree(rvals); 141 PLogObjectMemory(A,-ec*sizeof(int)); 142 ierr = MatDestroy(B); CHKERRQ(ierr); 143 PLogObjectParent(A,Bnew); 144 baij->B = Bnew; 145 A->was_assembled = PETSC_FALSE; 146 return 0; 147 } 148 149 150