1cb512458SBarry Smith #ifndef lint 2*9e25ed09SBarry Smith static char vcid[] = "$Id: mmaij.c,v 1.5 1995/03/06 04:05:15 bsmith Exp bsmith $"; 3cb512458SBarry Smith #endif 48c79f6d3SBarry Smith 5d6dfbf8fSBarry Smith 68c79f6d3SBarry Smith /* 78c79f6d3SBarry Smith Support for the parallel AIJ matrix vector multiply 88c79f6d3SBarry Smith */ 98c79f6d3SBarry Smith #include "mpiaij.h" 108c79f6d3SBarry Smith #include "vec/vecimpl.h" 118c79f6d3SBarry Smith #include "../seq/aij.h" 128c79f6d3SBarry Smith 138c79f6d3SBarry Smith int MPIAIJSetUpMultiply(Mat mat) 148c79f6d3SBarry Smith { 158c79f6d3SBarry Smith Matimpiaij *aij = (Matimpiaij *) mat->data; 168c79f6d3SBarry Smith Matiaij *B = (Matiaij *) (aij->B->data); 176abc6512SBarry Smith int N = aij->N,i,j,*indices,*aj = B->j; 181eb62cbbSBarry Smith int ierr,ec = 0,*garray; 191eb62cbbSBarry Smith IS from,to; 201eb62cbbSBarry Smith Vec gvec; 218c79f6d3SBarry Smith 228c79f6d3SBarry Smith /* For the first stab we make an array as long as the number of columns */ 231eb62cbbSBarry Smith /* mark those columns that are in aij->B */ 248c79f6d3SBarry Smith indices = (int *) MALLOC( N*sizeof(int) ); CHKPTR(indices); 258c79f6d3SBarry Smith MEMSET(indices,0,N*sizeof(int)); 26d6dfbf8fSBarry Smith for ( i=0; i<B->m; i++ ) { 27d6dfbf8fSBarry Smith for ( j=0; j<B->ilen[i]; j++ ) { 28d6dfbf8fSBarry Smith if (!indices[aj[B->i[i] - 1 + j]-1]) ec++; 29d6dfbf8fSBarry Smith indices[aj[B->i[i] - 1 + j]-1] = 1;} 301eb62cbbSBarry Smith } 318c79f6d3SBarry Smith 321eb62cbbSBarry Smith /* form array of columns we need */ 331eb62cbbSBarry Smith garray = (int *) MALLOC( ec*sizeof(int) ); CHKPTR(garray); 341eb62cbbSBarry Smith ec = 0; 351eb62cbbSBarry Smith for ( i=0; i<N; i++ ) { 361eb62cbbSBarry Smith if (indices[i]) garray[ec++] = i; 371eb62cbbSBarry Smith } 381eb62cbbSBarry Smith 391eb62cbbSBarry Smith /* make indices now point into garray */ 401eb62cbbSBarry Smith for ( i=0; i<ec; i++ ) { 411eb62cbbSBarry Smith indices[garray[i]] = i+1; 421eb62cbbSBarry Smith } 431eb62cbbSBarry Smith 441eb62cbbSBarry Smith /* compact out the extra columns in B */ 45d6dfbf8fSBarry Smith for ( i=0; i<B->m; i++ ) { 46d6dfbf8fSBarry Smith for ( j=0; j<B->ilen[i]; j++ ) { 47d6dfbf8fSBarry Smith aj[B->i[i] - 1 + j] = indices[aj[B->i[i] - 1 + j]-1]; 481eb62cbbSBarry Smith } 49d6dfbf8fSBarry Smith } 50d6dfbf8fSBarry Smith B->n = ec; 511eb62cbbSBarry Smith FREE(indices); 521eb62cbbSBarry Smith 531eb62cbbSBarry Smith /* create local vector that is used to scatter into */ 541eb62cbbSBarry Smith ierr = VecCreateSequential(ec,&aij->lvec); CHKERR(ierr); 551eb62cbbSBarry Smith 56d6dfbf8fSBarry Smith /* create two temporary Index sets for build scatter gather */ 571eb62cbbSBarry Smith ierr = ISCreateSequential(ec,garray,&from); CHKERR(ierr); 581eb62cbbSBarry Smith ierr = ISCreateStrideSequential(ec,0,1,&to); CHKERR(ierr); 591eb62cbbSBarry Smith 601eb62cbbSBarry Smith /* create temporary global vector to generate scatter context */ 611eb62cbbSBarry Smith /* this is inefficient, but otherwise we must do either 621eb62cbbSBarry Smith 1) save garray until the first actual scatter when the vector is known or 631eb62cbbSBarry Smith 2) have another way of generating a scatter context without a vector.*/ 64d6dfbf8fSBarry Smith ierr = VecCreateMPI(mat->comm,aij->n,aij->N,&gvec); CHKERR(ierr); 651eb62cbbSBarry Smith 661eb62cbbSBarry Smith /* gnerate the scatter context */ 671eb62cbbSBarry Smith ierr = VecScatterCtxCreate(gvec,from,aij->lvec,to,&aij->Mvctx); CHKERR(ierr); 681eb62cbbSBarry Smith 69*9e25ed09SBarry Smith aij->garray = garray; 701eb62cbbSBarry Smith ierr = ISDestroy(from); CHKERR(ierr); 711eb62cbbSBarry Smith ierr = ISDestroy(to); CHKERR(ierr); 721eb62cbbSBarry Smith ierr = VecDestroy(gvec); 738c79f6d3SBarry Smith return 0; 748c79f6d3SBarry Smith } 75*9e25ed09SBarry Smith 76*9e25ed09SBarry Smith 77