18c79f6d3SBarry Smith 2*d6dfbf8fSBarry Smith 38c79f6d3SBarry Smith /* 48c79f6d3SBarry Smith Support for the parallel AIJ matrix vector multiply 58c79f6d3SBarry Smith */ 68c79f6d3SBarry Smith #include "mpiaij.h" 78c79f6d3SBarry Smith #include "vec/vecimpl.h" 88c79f6d3SBarry Smith #include "../seq/aij.h" 98c79f6d3SBarry Smith 108c79f6d3SBarry Smith int MPIAIJSetUpMultiply(Mat mat) 118c79f6d3SBarry Smith { 128c79f6d3SBarry Smith Matimpiaij *aij = (Matimpiaij *) mat->data; 138c79f6d3SBarry Smith Matiaij *B = (Matiaij *) (aij->B->data); 141eb62cbbSBarry Smith int N = aij->N,i,j,*indices,n = B->i[B->m] - 1,*aj = B->j; 151eb62cbbSBarry Smith int ierr,ec = 0,*garray; 161eb62cbbSBarry Smith IS from,to; 171eb62cbbSBarry Smith Vec gvec; 188c79f6d3SBarry Smith 198c79f6d3SBarry Smith /* For the first stab we make an array as long as the number of columns */ 201eb62cbbSBarry Smith /* mark those columns that are in aij->B */ 218c79f6d3SBarry Smith indices = (int *) MALLOC( N*sizeof(int) ); CHKPTR(indices); 228c79f6d3SBarry Smith MEMSET(indices,0,N*sizeof(int)); 23*d6dfbf8fSBarry Smith for ( i=0; i<B->m; i++ ) { 24*d6dfbf8fSBarry Smith for ( j=0; j<B->ilen[i]; j++ ) { 25*d6dfbf8fSBarry Smith if (!indices[aj[B->i[i] - 1 + j]-1]) ec++; 26*d6dfbf8fSBarry Smith indices[aj[B->i[i] - 1 + j]-1] = 1;} 271eb62cbbSBarry Smith } 288c79f6d3SBarry Smith 291eb62cbbSBarry Smith /* form array of columns we need */ 301eb62cbbSBarry Smith garray = (int *) MALLOC( ec*sizeof(int) ); CHKPTR(garray); 311eb62cbbSBarry Smith ec = 0; 321eb62cbbSBarry Smith for ( i=0; i<N; i++ ) { 331eb62cbbSBarry Smith if (indices[i]) garray[ec++] = i; 341eb62cbbSBarry Smith } 351eb62cbbSBarry Smith 361eb62cbbSBarry Smith /* make indices now point into garray */ 371eb62cbbSBarry Smith for ( i=0; i<ec; i++ ) { 381eb62cbbSBarry Smith indices[garray[i]] = i+1; 391eb62cbbSBarry Smith } 401eb62cbbSBarry Smith 411eb62cbbSBarry Smith /* compact out the extra columns in B */ 42*d6dfbf8fSBarry Smith for ( i=0; i<B->m; i++ ) { 43*d6dfbf8fSBarry Smith for ( j=0; j<B->ilen[i]; j++ ) { 44*d6dfbf8fSBarry Smith aj[B->i[i] - 1 + j] = indices[aj[B->i[i] - 1 + j]-1]; 451eb62cbbSBarry Smith } 46*d6dfbf8fSBarry Smith } 47*d6dfbf8fSBarry Smith B->n = ec; 481eb62cbbSBarry Smith FREE(indices); 491eb62cbbSBarry Smith 501eb62cbbSBarry Smith /* create local vector that is used to scatter into */ 511eb62cbbSBarry Smith ierr = VecCreateSequential(ec,&aij->lvec); CHKERR(ierr); 521eb62cbbSBarry Smith 53*d6dfbf8fSBarry Smith /* create two temporary Index sets for build scatter gather */ 541eb62cbbSBarry Smith ierr = ISCreateSequential(ec,garray,&from); CHKERR(ierr); 551eb62cbbSBarry Smith ierr = ISCreateStrideSequential(ec,0,1,&to); CHKERR(ierr); 561eb62cbbSBarry Smith 571eb62cbbSBarry Smith /* create temporary global vector to generate scatter context */ 581eb62cbbSBarry Smith /* this is inefficient, but otherwise we must do either 591eb62cbbSBarry Smith 1) save garray until the first actual scatter when the vector is known or 601eb62cbbSBarry Smith 2) have another way of generating a scatter context without a vector.*/ 61*d6dfbf8fSBarry Smith ierr = VecCreateMPI(mat->comm,aij->n,aij->N,&gvec); CHKERR(ierr); 621eb62cbbSBarry Smith 631eb62cbbSBarry Smith /* gnerate the scatter context */ 641eb62cbbSBarry Smith ierr = VecScatterCtxCreate(gvec,from,aij->lvec,to,&aij->Mvctx); CHKERR(ierr); 651eb62cbbSBarry Smith 661eb62cbbSBarry Smith FREE(garray); 671eb62cbbSBarry Smith ierr = ISDestroy(from); CHKERR(ierr); 681eb62cbbSBarry Smith ierr = ISDestroy(to); CHKERR(ierr); 691eb62cbbSBarry Smith ierr = VecDestroy(gvec); 708c79f6d3SBarry Smith return 0; 718c79f6d3SBarry Smith } 72