Actual source code: ex4.c

petsc-master 2016-07-29
Report Typos and Errors
  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);if (ierr) return ierr;
 30:   /*
 31:      Determine file from which we read the matrix

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


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

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

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

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

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

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

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

 87:   PetscFinalize();
 88:   return ierr;
 89: }

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

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

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

101:      Note: this routine directly access the Vec and Mat data-structures
102: */
103: PetscErrorCode LowRankUpdate(Mat U,Mat V,Vec x,Vec y,Vec work1,Vec work2,PetscInt nwork)
104: {
105:   Mat            Ulocal = ((Mat_MPIDense*)U->data)->A;
106:   Mat            Vlocal = ((Mat_MPIDense*)V->data)->A;
107:   PetscInt       n;
109:   PetscScalar    *w1,*w2,*array;
110:   Vec            xseq,yseq;
111: 
113:   /* First multiply the local part of U with the local part of x */
114:   /* this tricks the silly error checking in MatMultTranspose();*/
115:   VecGetLocalSize(x,&n);
116:   VecGetArray(x,&array);
117:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,array,&xseq);
118:   MatMultTranspose(Ulocal,xseq,work1); /* note in this call x is treated as a sequential vector  */
119:   VecRestoreArray(x,&array);
120:   VecDestroy(&xseq);
121: 
122:   /* Form the sum of all the local multiplies : this is work2 = U'*x = sum_{all processors} work1 */
123:   VecGetArray(work1,&w1);
124:   VecGetArray(work2,&w2);
125:   MPI_Allreduce(w1,w2,nwork,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
126:   VecRestoreArray(work1,&w1);
127:   VecRestoreArray(work2,&w2);

129:   /* multiply y = V*work2 */
130:   /* this tricks the silly error checking in MatMult() */
131:   VecGetLocalSize(y,&n);
132:   VecGetArray(y,&array);
133:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,array,&yseq);
134:   MatMult(Vlocal,work2,y);
135:   VecRestoreArray(y,&array);
136:   VecDestroy(&yseq);
137:   return(0);
138: }