2: #include <petsc/private/matimpl.h>          /*I "petscmat.h" I*/

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


 11: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqDense(Mat,Vec,Vec);
 12: PETSC_INTERN PetscErrorCode MatMultAdd_SeqDense(Mat,Vec,Vec,Vec);

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

 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:   MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
 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:   MatDestroy(&Na->A);
 53:   MatDestroy(&Na->U);
 54:   MatDestroy(&Na->V);
 55:   VecDestroy(&Na->work1);
 56:   VecDestroy(&Na->work2);
 57:   PetscFree(N->data);
 58:   return(0);
 59: }


 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' is not 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(PetscObjectComm((PetscObject)A),N);
 91:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
 92:   PetscObjectChangeTypeName((PetscObject)*N,MATLRC);

 94:   PetscNewLog(*N,&Na);
 95:   (*N)->data = (void*) Na;
 96:   Na->A      = A;

 98:   MatDenseGetLocalMatrix(U,&Na->U);
 99:   MatDenseGetLocalMatrix(V,&Na->V);
100:   PetscObjectReference((PetscObject)A);
101:   PetscObjectReference((PetscObject)Na->U);
102:   PetscObjectReference((PetscObject)Na->V);

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: }