1cb512458SBarry Smith #ifndef lint 2*78b31e54SBarry Smith static char vcid[] = "$Id: mmaij.c,v 1.10 1995/05/27 20:23: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 1344a69424SLois Curfman McInnes int MatSetUpMultiply_MPIAIJ(Mat mat) 148c79f6d3SBarry Smith { 1544a69424SLois Curfman McInnes Mat_MPIAIJ *aij = (Mat_MPIAIJ *) mat->data; 1644a69424SLois 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 */ 24*78b31e54SBarry Smith indices = (int *) PETSCMALLOC( N*sizeof(int) ); CHKPTRQ(indices); 25*78b31e54SBarry Smith PETSCMEMSET(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 */ 33*78b31e54SBarry Smith garray = (int *) PETSCMALLOC( (ec+1)*sizeof(int) ); CHKPTRQ(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; 51*78b31e54SBarry Smith PETSCFREE(indices); 521eb62cbbSBarry Smith 531eb62cbbSBarry Smith /* create local vector that is used to scatter into */ 54*78b31e54SBarry Smith ierr = VecCreateSequential(MPI_COMM_SELF,ec,&aij->lvec); CHKERRQ(ierr); 551eb62cbbSBarry Smith 56d6dfbf8fSBarry Smith /* create two temporary Index sets for build scatter gather */ 57*78b31e54SBarry Smith ierr = ISCreateSequential(MPI_COMM_SELF,ec,garray,&from); CHKERRQ(ierr); 58*78b31e54SBarry Smith ierr = ISCreateStrideSequential(MPI_COMM_SELF,ec,0,1,&to); CHKERRQ(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.*/ 64*78b31e54SBarry Smith ierr = VecCreateMPI(mat->comm,aij->n,aij->N,&gvec); CHKERRQ(ierr); 651eb62cbbSBarry Smith 661eb62cbbSBarry Smith /* gnerate the scatter context */ 67*78b31e54SBarry Smith ierr = VecScatterCtxCreate(gvec,from,aij->lvec,to,&aij->Mvctx); CHKERRQ(ierr); 68a5a9c739SBarry Smith PLogObjectParent(mat,aij->Mvctx); 69a5a9c739SBarry Smith PLogObjectParent(mat,aij->lvec); 709e25ed09SBarry Smith aij->garray = garray; 71*78b31e54SBarry Smith ierr = ISDestroy(from); CHKERRQ(ierr); 72*78b31e54SBarry Smith ierr = ISDestroy(to); CHKERRQ(ierr); 731eb62cbbSBarry Smith ierr = VecDestroy(gvec); 748c79f6d3SBarry Smith return 0; 758c79f6d3SBarry Smith } 769e25ed09SBarry Smith 779e25ed09SBarry Smith 78