2: #include <petsc/private/matimpl.h>          /*I "petscmat.h" I*/
  3: #include <../src/mat/impls/dense/seq/dense.h>

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



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

 22:   MatMult(Na->A,x,y);

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

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

 31:   VecGetArray(Na->work1,&w1);
 32:   VecGetArray(Na->work2,&w2);
 33:   MPI_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
 34:   VecRestoreArray(Na->work1,&w1);
 35:   VecRestoreArray(Na->work2,&w2);

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

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

 51:   MatDestroy(&Na->A);
 52:   MatDestroy(&Na->U);
 53:   MatDestroy(&Na->V);
 54:   VecDestroy(&Na->work1);
 55:   VecDestroy(&Na->work2);
 56:   PetscFree(N->data);
 57:   return(0);
 58: }


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

 66:    Collective on Mat

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

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

 75:    Level: intermediate

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

 88:   MatGetLocalSize(A,&m,&n);
 89:   MatCreate(PetscObjectComm((PetscObject)A),N);
 90:   MatSetSizes(*N,n,n,PETSC_DECIDE,PETSC_DECIDE);
 91:   PetscObjectChangeTypeName((PetscObject)*N,MATLRC);

 93:   PetscNewLog(*N,&Na);
 94:   (*N)->data = (void*) 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);

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

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