1 2 3 /* 4 Support for the parallel AIJ matrix vector multiply 5 */ 6 #include "mpiaij.h" 7 #include "vec/vecimpl.h" 8 #include "../seq/aij.h" 9 10 int MPIAIJSetUpMultiply(Mat mat) 11 { 12 Matimpiaij *aij = (Matimpiaij *) mat->data; 13 Matiaij *B = (Matiaij *) (aij->B->data); 14 int N = aij->N,i,j,*indices,*aj = B->j; 15 int ierr,ec = 0,*garray; 16 IS from,to; 17 Vec gvec; 18 19 /* For the first stab we make an array as long as the number of columns */ 20 /* mark those columns that are in aij->B */ 21 indices = (int *) MALLOC( N*sizeof(int) ); CHKPTR(indices); 22 MEMSET(indices,0,N*sizeof(int)); 23 for ( i=0; i<B->m; i++ ) { 24 for ( j=0; j<B->ilen[i]; j++ ) { 25 if (!indices[aj[B->i[i] - 1 + j]-1]) ec++; 26 indices[aj[B->i[i] - 1 + j]-1] = 1;} 27 } 28 29 /* form array of columns we need */ 30 garray = (int *) MALLOC( ec*sizeof(int) ); CHKPTR(garray); 31 ec = 0; 32 for ( i=0; i<N; i++ ) { 33 if (indices[i]) garray[ec++] = i; 34 } 35 36 /* make indices now point into garray */ 37 for ( i=0; i<ec; i++ ) { 38 indices[garray[i]] = i+1; 39 } 40 41 /* compact out the extra columns in B */ 42 for ( i=0; i<B->m; i++ ) { 43 for ( j=0; j<B->ilen[i]; j++ ) { 44 aj[B->i[i] - 1 + j] = indices[aj[B->i[i] - 1 + j]-1]; 45 } 46 } 47 B->n = ec; 48 FREE(indices); 49 50 /* create local vector that is used to scatter into */ 51 ierr = VecCreateSequential(ec,&aij->lvec); CHKERR(ierr); 52 53 /* create two temporary Index sets for build scatter gather */ 54 ierr = ISCreateSequential(ec,garray,&from); CHKERR(ierr); 55 ierr = ISCreateStrideSequential(ec,0,1,&to); CHKERR(ierr); 56 57 /* create temporary global vector to generate scatter context */ 58 /* this is inefficient, but otherwise we must do either 59 1) save garray until the first actual scatter when the vector is known or 60 2) have another way of generating a scatter context without a vector.*/ 61 ierr = VecCreateMPI(mat->comm,aij->n,aij->N,&gvec); CHKERR(ierr); 62 63 /* gnerate the scatter context */ 64 ierr = VecScatterCtxCreate(gvec,from,aij->lvec,to,&aij->Mvctx); CHKERR(ierr); 65 66 FREE(garray); 67 ierr = ISDestroy(from); CHKERR(ierr); 68 ierr = ISDestroy(to); CHKERR(ierr); 69 ierr = VecDestroy(gvec); 70 return 0; 71 } 72