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