Actual source code: ex4.c

  1: /*$Id: ex1.c,v 1.28 2001/04/10 19:35:58 bsmith Exp $*/

  3: static char help[] = "Reads U and V matrices from a file and performs y = V*U'*x.\n\
  4:   -f <input_file> : file to load \n\n";

  6: /* 
  7:   Include "petscmat.h" so that we can use matrices.
  8:   automatically includes:
  9:      petsc.h       - base PETSc routines   petscvec.h    - vectors
 10:      petscsys.h    - system routines       petscmat.h    - matrices
 11:      petscis.h     - index sets            petscviewer.h - viewers               
 12: */
 13:  #include petscmat.h
 14: extern int LowRankUpdate(Mat,Mat,Vec,Vec,Vec,Vec,int);


 19: int main(int argc,char **args)
 20: {
 21:   Mat                   U,V;                /* matrix */
 22:   PetscViewer           fd;               /* viewer */
 23:   char                  file[128];     /* input file name */
 24:   int                   ierr;
 25:   PetscTruth            flg;
 26:   Vec                   x,y,work1,work2;
 27:   int                   i,N,n,M,m;
 28:   PetscScalar           *xx;

 30:   PetscInitialize(&argc,&args,(char *)0,help);

 32:   /* 
 33:      Determine file from which we read the matrix

 35:   */
 36:   PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);
 37:   if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");


 40:   /* 
 41:      Open binary file.  Note that we use PETSC_FILE_RDONLY to indicate
 42:      reading from this file.
 43:   */
 44:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);

 46:   /*
 47:     Load the matrix; then destroy the viewer.
 48:     Note both U and V are stored as tall skinny matrices 
 49:   */
 50:   MatLoad(fd,MATMPIDENSE,&U);
 51:   MatLoad(fd,MATMPIDENSE,&V);
 52:   PetscViewerDestroy(fd);

 54:   MatGetLocalSize(U,&N,&n);
 55:   MatGetLocalSize(V,&M,&m);
 56:   if (N != M) SETERRQ2(1,"U and V matrices must have same number of local rows %d %d",N,M);
 57:   if (n != m) SETERRQ2(1,"U and V matrices must have same number of local columns %d %d",n,m);

 59:   VecCreateMPI(PETSC_COMM_WORLD,N,PETSC_DETERMINE,&x);
 60:   VecDuplicate(x,&y);

 62:   MatGetSize(U,0,&n);
 63:   VecCreateSeq(PETSC_COMM_SELF,n,&work1);
 64:   VecDuplicate(work1,&work2);

 66:   /* put some initial values into x for testing */
 67:   VecGetArray(x,&xx);
 68:   for (i=0; i<N; i++) xx[i] = i;
 69:   VecRestoreArray(x,&xx);
 70:   LowRankUpdate(U,V,x,y,work1,work2,n);
 71:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
 72:   VecView(y,PETSC_VIEWER_STDOUT_WORLD);

 74:   /* 
 75:      Free work space.  All PETSc objects should be destroyed when they
 76:      are no longer needed.
 77:   */
 78:   MatDestroy(U);
 79:   MatDestroy(V);
 80:   VecDestroy(x);
 81:   VecDestroy(y);
 82:   VecDestroy(work1);
 83:   VecDestroy(work2);

 85:   PetscFinalize();
 86:   return 0;
 87: }

 89:  #include src/mat/impls/dense/mpi/mpidense.h

 93: /*
 94:      Computes y = V*U'*x where U and V are  N by n (N >> n). 

 96:      U and V are stored as PETSc MPIDENSE (parallel) dense matrices with their rows partitioned the
 97:      same way as x and y are partitioned
 98: */
 99: int LowRankUpdate(Mat U,Mat V,Vec x,Vec y,Vec work1,Vec work2,int nwork)
100: {
101:   Mat         Ulocal = ((Mat_MPIDense*)U->data)->A;
102:   Mat         Vlocal = ((Mat_MPIDense*)V->data)->A;
103:   int         ierr,Nsave = x->N;
104:   PetscScalar *w1,*w2;


108:   /* First multiply the local part of U with the local part of x */
109:   x->N = x->n; /* this tricks the silly error checking in MatMultTranspose();*/
110:   MatMultTranspose(Ulocal,x,work1);/* note in this call x is treated as a sequential vector  */
111:   x->N = Nsave;

113:   /* Form the sum of all the local multiplies : this is work2 = U'*x = sum_{all processors} work1 */
114:   VecGetArray(work1,&w1);
115:   VecGetArray(work2,&w2);
116:   MPI_Allreduce(w1,w2,nwork,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
117:   VecRestoreArray(work1,&w1);
118:   VecRestoreArray(work2,&w2);

120:   /* multiply y = V*work2 */
121:   y->N = y->n; /* this tricks the silly error checking in MatMult() */
122:   MatMult(Vlocal,work2,y);/* note in this call y is treated as a sequential vector  */
123:   y->N = Nsave;

125:   return(0);
126: }