1 #ifndef lint 2 static char vcid[] = "$Id: mmbaij.c,v 1.6 1996/08/08 14:43:40 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 PLogObjectParent(mat,baij->Mvctx); 86 PLogObjectParent(mat,baij->lvec); 87 PLogObjectParent(mat,from); 88 PLogObjectParent(mat,to); 89 baij->garray = garray; 90 PLogObjectMemory(mat,(ec+1)*sizeof(int)); 91 ierr = ISDestroy(from); CHKERRQ(ierr); 92 ierr = ISDestroy(to); CHKERRQ(ierr); 93 ierr = VecDestroy(gvec); 94 PetscFree(tmp); 95 return 0; 96 } 97 98 99 /* 100 Takes the local part of an already assembled MPIBAIJ matrix 101 and disassembles it. This is to allow new nonzeros into the matrix 102 that require more communication in the matrix vector multiply. 103 Thus certain data-structures must be rebuilt. 104 105 Kind of slow! But that's what application programmers get when 106 they are sloppy. 107 */ 108 int DisAssemble_MPIBAIJ(Mat A) 109 { 110 Mat_MPIBAIJ *baij = (Mat_MPIBAIJ *) A->data; 111 Mat B = baij->B,Bnew; 112 Mat_SeqBAIJ *Bbaij = (Mat_SeqBAIJ*)B->data; 113 int ierr,i,j,mbs=Bbaij->mbs,n = baij->N,col,*garray=baij->garray; 114 int k,bs=baij->bs,bs2=baij->bs2,*rvals,*nz,ec,m=Bbaij->m; 115 Scalar *a=Bbaij->a; 116 117 /* free stuff related to matrix-vec multiply */ 118 ierr = VecGetSize(baij->lvec,&ec); /* needed for PLogObjectMemory below */ 119 ierr = VecDestroy(baij->lvec); CHKERRQ(ierr); baij->lvec = 0; 120 ierr = VecScatterDestroy(baij->Mvctx); CHKERRQ(ierr); baij->Mvctx = 0; 121 if (baij->colmap) { 122 PetscFree(baij->colmap); baij->colmap = 0; 123 PLogObjectMemory(A,-Bbaij->nbs*sizeof(int)); 124 } 125 126 /* make sure that B is assembled so we can access its values */ 127 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 128 MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 129 130 /* invent new B and copy stuff over */ 131 nz = (int *) PetscMalloc( mbs*sizeof(int) ); CHKPTRQ(nz); 132 for ( i=0; i<mbs; i++ ) { 133 nz[i] = Bbaij->i[i+1]-Bbaij->i[i]; 134 } 135 ierr = MatCreateSeqBAIJ(MPI_COMM_SELF,baij->bs,m,n,0,nz,&Bnew); CHKERRQ(ierr); 136 PetscFree(nz); 137 138 rvals = (int *) PetscMalloc(bs*sizeof(int)); CHKPTRQ(rvals); 139 for ( i=0; i<mbs; i++ ) { 140 rvals[0] = bs*i; 141 for ( j=1; j<bs; j++ ) { rvals[j] = rvals[j-1] + 1; } 142 for ( j=Bbaij->i[i]; j<Bbaij->i[i+1]; j++ ) { 143 col = garray[Bbaij->j[j]]*bs; 144 for (k=0; k<bs; k++ ) { 145 ierr = MatSetValues(Bnew,bs,rvals,1,&col,a+j*bs2,INSERT_VALUES);CHKERRQ(ierr); 146 col++; 147 } 148 } 149 } 150 PetscFree(baij->garray); baij->garray = 0; 151 PetscFree(rvals); 152 PLogObjectMemory(A,-ec*sizeof(int)); 153 ierr = MatDestroy(B); CHKERRQ(ierr); 154 PLogObjectParent(A,Bnew); 155 baij->B = Bnew; 156 A->was_assembled = PETSC_FALSE; 157 return 0; 158 } 159 160 161