1cb512458SBarry Smith #ifndef lint 2*44a69424SLois Curfman McInnes static char vcid[] = "$Id: mmaij.c,v 1.7 1995/03/21 23:19:19 bsmith Exp curfman $"; 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 13*44a69424SLois Curfman McInnes int MatSetUpMultiply_MPIAIJ(Mat mat) 148c79f6d3SBarry Smith { 15*44a69424SLois Curfman McInnes Mat_MPIAIJ *aij = (Mat_MPIAIJ *) mat->data; 16*44a69424SLois Curfman McInnes Mat_AIJ *B = (Mat_AIJ *) (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); 68a5a9c739SBarry Smith PLogObjectParent(mat,aij->Mvctx); 69a5a9c739SBarry Smith PLogObjectParent(mat,aij->lvec); 709e25ed09SBarry Smith aij->garray = garray; 711eb62cbbSBarry Smith ierr = ISDestroy(from); CHKERR(ierr); 721eb62cbbSBarry Smith ierr = ISDestroy(to); CHKERR(ierr); 731eb62cbbSBarry Smith ierr = VecDestroy(gvec); 748c79f6d3SBarry Smith return 0; 758c79f6d3SBarry Smith } 769e25ed09SBarry Smith 779e25ed09SBarry Smith 78