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