1 #ifndef lint 2 static char vcid[] = "$Id: mmbaij.c,v 1.4 1996/07/08 22:20:04 bsmith Exp bsmith $"; 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 building scatter-gather */ 63 64 /* ierr = ISCreateSeq(MPI_COMM_SELF,ec*bs,tmp,&from); CHKERRQ(ierr); */ 65 for ( i=0,col=0; i<ec; i++ ) { 66 garray[i] = bs*garray[i]; 67 } 68 ierr = ISCreateBlockSeq(MPI_COMM_SELF,bs,ec,garray,&from);CHKERRQ(ierr); 69 for ( i=0,col=0; i<ec; i++ ) { 70 garray[i] = garray[i]/bs; 71 } 72 73 ierr = ISCreateStrideSeq(MPI_COMM_SELF,ec*bs,0,1,&to);CHKERRQ(ierr); 74 75 /* create temporary global vector to generate scatter context */ 76 /* this is inefficient, but otherwise we must do either 77 1) save garray until the first actual scatter when the vector is known or 78 2) have another way of generating a scatter context without a vector.*/ 79 ierr = VecCreateMPI(mat->comm,baij->n,baij->N,&gvec); CHKERRQ(ierr); 80 81 /* gnerate the scatter context */ 82 ierr = VecScatterCreate(gvec,from,baij->lvec,to,&baij->Mvctx);CHKERRQ(ierr); 83 PLogObjectParent(mat,baij->Mvctx); 84 PLogObjectParent(mat,baij->lvec); 85 PLogObjectParent(mat,from); 86 PLogObjectParent(mat,to); 87 baij->garray = garray; 88 PLogObjectMemory(mat,(ec+1)*sizeof(int)); 89 ierr = ISDestroy(from); CHKERRQ(ierr); 90 ierr = ISDestroy(to); CHKERRQ(ierr); 91 ierr = VecDestroy(gvec); 92 PetscFree(tmp); 93 return 0; 94 } 95 96 97 /* 98 Takes the local part of an already assembled MPIBAIJ matrix 99 and disassembles it. This is to allow new nonzeros into the matrix 100 that require more communication in the matrix vector multiply. 101 Thus certain data-structures must be rebuilt. 102 103 Kind of slow! But that's what application programmers get when 104 they are sloppy. 105 */ 106 int DisAssemble_MPIBAIJ(Mat A) 107 { 108 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *) A->data; 109 Mat B = baij->B,Bnew; 110 Mat_SeqBAIJ *Bbaij = (Mat_SeqBAIJ*)B->data; 111 int ierr,i,j,mbs=Bbaij->mbs,n = baij->N,col,*garray=baij->garray; 112 int k,bs=baij->bs,bs2=baij->bs2,*rvals,*nz,ec,m=Bbaij->m; 113 Scalar *a=Bbaij->a; 114 115 /* free stuff related to matrix-vec multiply */ 116 ierr = VecGetSize(baij->lvec,&ec); /* needed for PLogObjectMemory below */ 117 ierr = VecDestroy(baij->lvec); CHKERRQ(ierr); baij->lvec = 0; 118 ierr = VecScatterDestroy(baij->Mvctx); CHKERRQ(ierr); baij->Mvctx = 0; 119 if (baij->colmap) { 120 PetscFree(baij->colmap); baij->colmap = 0; 121 PLogObjectMemory(A,-Bbaij->nbs*sizeof(int)); 122 } 123 124 /* make sure that B is assembled so we can access its values */ 125 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 126 MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 127 128 /* invent new B and copy stuff over */ 129 nz = (int *) PetscMalloc( mbs*sizeof(int) ); CHKPTRQ(nz); 130 for ( i=0; i<mbs; i++ ) { 131 nz[i] = Bbaij->i[i+1]-Bbaij->i[i]; 132 } 133 ierr = MatCreateSeqBAIJ(MPI_COMM_SELF,baij->bs,m,n,0,nz,&Bnew); CHKERRQ(ierr); 134 PetscFree(nz); 135 136 rvals = (int *) PetscMalloc(bs*sizeof(int)); CHKPTRQ(rvals); 137 for ( i=0; i<mbs; i++ ) { 138 rvals[0] = bs*i; 139 for ( j=1; j<bs; j++ ) { rvals[j] = rvals[j-1] + 1; } 140 for ( j=Bbaij->i[i]; j<Bbaij->i[i+1]; j++ ) { 141 col = garray[Bbaij->j[j]]*bs; 142 for (k=0; k<bs; k++ ) { 143 ierr = MatSetValues(Bnew,bs,rvals,1,&col,a+j*bs2,INSERT_VALUES);CHKERRQ(ierr); 144 col++; 145 } 146 } 147 } 148 PetscFree(baij->garray); baij->garray = 0; 149 PetscFree(rvals); 150 PLogObjectMemory(A,-ec*sizeof(int)); 151 ierr = MatDestroy(B); CHKERRQ(ierr); 152 PLogObjectParent(A,Bnew); 153 baij->B = Bnew; 154 A->was_assembled = PETSC_FALSE; 155 return 0; 156 } 157 158 159