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