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 */
 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(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>
 93: #include <petsc/private/vecimpl.h>

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

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

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

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

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

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