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