1: #define PETSCMAT_DLL

 3:  #include include/private/matimpl.h
 4:  #include src/mat/impls/dense/seq/dense.h

  6: typedef struct {
  7:   Mat         A,U,V;
  8:   Vec         work1,work2;/* Sequential (big) vectors that hold partial products */
  9:   PetscMPIInt nwork;      /* length of work vectors */
 10: } Mat_LRC;



 16: PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y)
 17: {
 18:   Mat_LRC        *Na = (Mat_LRC*)N->data;
 20:   PetscScalar    *w1,*w2;
 21: 
 23:   MatMult(Na->A,x,y);

 25:   /* multiply the local part of V with the local part of x */
 26:   /* note in this call x is treated as a sequential vector  */
 27:   MatMultTranspose_SeqDense(Na->V,x,Na->work1);

 29:   /* Form the sum of all the local multiplies : this is work2 = V'*x =
 30:      sum_{all processors} work1 */

 32:   VecGetArray(Na->work1,&w1);
 33:   VecGetArray(Na->work2,&w2);
 34:   MPI_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,PetscSum_Op,N->comm);
 35:   VecRestoreArray(Na->work1,&w1);
 36:   VecRestoreArray(Na->work2,&w2);

 38:   /* multiply-sub y = y  + U*work2 */
 39:   /* note in this call y is treated as a sequential vector  */
 40:   MatMultAdd_SeqDense(Na->U,Na->work2,y,y);
 41:   return(0);
 42: }

 46: PetscErrorCode MatDestroy_LRC(Mat N)
 47: {
 48:   Mat_LRC        *Na = (Mat_LRC*)N->data;

 52:   PetscObjectDereference((PetscObject)Na->A);
 53:   PetscObjectDereference((PetscObject)Na->U);
 54:   PetscObjectDereference((PetscObject)Na->V);
 55:   VecDestroy(Na->work1);
 56:   VecDestroy(Na->work2);
 57:   PetscFree(Na);
 58:   return(0);
 59: }
 60: 

 64: /*@
 65:       MatCreateLRC - Creates a new matrix object that behaves like A + U*V'

 67:    Collective on Mat

 69:    Input Parameter:
 70: +   A  - the (sparse) matrix
 71: -   U. V - two dense rectangular (tall and skinny) matrices

 73:    Output Parameter:
 74: .   N - the matrix that represents A + U*V'

 76:    Level: intermediate

 78:    Notes: The matrix A + U*V' formed! Rather the new matrix
 79:           object performs the matrix-vector product by first multiplying by
 80:           A and then adding the other term
 81: @*/
 82: PetscErrorCode  MatCreateLRC(Mat A,Mat U, Mat V,Mat *N)
 83: {
 85:   PetscInt       m,n;
 86:   Mat_LRC        *Na;

 89:   MatGetLocalSize(A,&m,&n);
 90:   MatCreate(A->comm,N);
 91:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
 92:   PetscObjectChangeTypeName((PetscObject)*N,MATLRC);
 93: 
 94:   PetscNew(Mat_LRC,&Na);
 95:   Na->A     = A;

 97:   MatDenseGetLocalMatrix(U,&Na->U);
 98:   MatDenseGetLocalMatrix(V,&Na->V);
 99:   PetscObjectReference((PetscObject)A);
100:   PetscObjectReference((PetscObject)Na->U);
101:   PetscObjectReference((PetscObject)Na->V);
102:   (*N)->data = (void*) Na;

104:   VecCreateSeq(PETSC_COMM_SELF,U->cmap.N,&Na->work1);
105:   VecDuplicate(Na->work1,&Na->work2);
106:   Na->nwork              = U->cmap.N;

108:   (*N)->ops->destroy     = MatDestroy_LRC;
109:   (*N)->ops->mult        = MatMult_LRC;
110:   (*N)->assembled        = PETSC_TRUE;
111:   (*N)->cmap.N                = A->cmap.N;
112:   (*N)->rmap.N                = A->cmap.N;
113:   (*N)->cmap.n                = A->cmap.n;
114:   (*N)->rmap.n                = A->cmap.n;
115:   return(0);
116: }