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:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
  9:      petscmat.h    - matrices
 10:      petscis.h     - index sets            petscviewer.h - viewers               
 11: */
 12: #include <petscmat.h>
 13: extern PetscErrorCode LowRankUpdate(Mat,Mat,Vec,Vec,Vec,Vec,PetscInt);


 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:   PetscBool             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,&flg);
 36:   if (!flg) SETERRQ(PETSC_COMM_WORLD,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:   MatCreate(PETSC_COMM_WORLD,&U);
 50:   MatSetType(U,MATMPIDENSE);
 51:   MatLoad(U,fd);
 52:   MatCreate(PETSC_COMM_WORLD,&V);
 53:   MatSetType(V,MATMPIDENSE);
 54:   MatLoad(V,fd);
 55:   PetscViewerDestroy(&fd);

 57:   MatGetLocalSize(U,&N,&n);
 58:   MatGetLocalSize(V,&M,&m);
 59:   if (N != M) SETERRQ2(PETSC_COMM_SELF,1,"U and V matrices must have same number of local rows %D %D",N,M);
 60:   if (n != m) SETERRQ2(PETSC_COMM_SELF,1,"U and V matrices must have same number of local columns %D %D",n,m);

 62:   VecCreateMPI(PETSC_COMM_WORLD,N,PETSC_DETERMINE,&x);
 63:   VecDuplicate(x,&y);

 65:   MatGetSize(U,0,&n);
 66:   VecCreateSeq(PETSC_COMM_SELF,n,&work1);
 67:   VecDuplicate(work1,&work2);

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

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

 88:   PetscFinalize();
 89:   return 0;
 90: }

 92: #include <../src/mat/impls/dense/mpi/mpidense.h>

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

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

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

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

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

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

130:   return(0);
131: }