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