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