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

  5: /* 
  6:   Include "petscmat.h" so that we can use matrices.
  7:   automatically includes:
  8:      petsc.h       - base PETSc routines   petscvec.h    - vectors
  9:      petscsys.h    - system routines       petscmat.h    - matrices
 10:      petscis.h     - index sets            petscviewer.h - viewers               
 11: */
 12:  #include petscmat.h


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

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

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

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


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

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

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

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

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

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

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

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

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

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

 95:      U and V are stored as PETSc MPIDENSE (parallel) dense matrices with their rows partitioned the
 96:      same way as x and y are partitioned

 98:      Note: this routine directly access the Vec and Mat data-structures
 99: */
100: PetscErrorCode LowRankUpdate(Mat U,Mat V,Vec x,Vec y,Vec work1,Vec work2,PetscInt nwork)
101: {
102:   Mat            Ulocal = ((Mat_MPIDense*)U->data)->A;
103:   Mat            Vlocal = ((Mat_MPIDense*)V->data)->A;
104:   PetscInt       Nsave = x->map->N;
106:   PetscScalar    *w1,*w2;

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

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

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

126:   return(0);
127: }