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