1: #include <petsc/private/matimpl.h>

  3: static PetscErrorCode MatMult_Centering(Mat A, Vec xx, Vec yy)
  4: {
  5:   PetscScalar       *y;
  6:   const PetscScalar *x;
  7:   PetscScalar        sum, mean;
  8:   PetscInt           i, m = A->rmap->n, size;

 10:   PetscFunctionBegin;
 11:   PetscCall(VecSum(xx, &sum));
 12:   PetscCall(VecGetSize(xx, &size));
 13:   mean = sum / (PetscScalar)size;
 14:   PetscCall(VecGetArrayRead(xx, &x));
 15:   PetscCall(VecGetArray(yy, &y));
 16:   for (i = 0; i < m; i++) y[i] = x[i] - mean;
 17:   PetscCall(VecRestoreArrayRead(xx, &x));
 18:   PetscCall(VecRestoreArray(yy, &y));
 19:   PetscFunctionReturn(PETSC_SUCCESS);
 20: }

 22: /*MC
 23:   MATCENTERING - matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$

 25:   Level: beginner

 27:   Note:
 28:   See `MatCreateCentering()`

 30: .seealso:  `Mat`, `MatCreateCentering()`
 31: M*/

 33: /*@
 34:   MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$

 36:   Collective

 38:   Input Parameters:
 39: + comm - MPI communicator
 40: . n    - number of local rows (or `PETSC_DECIDE` to have calculated if `N` is given)
 41:          This value should be the same as the local size used in creating the
 42:          `y` vector for the matrix-vector product $y = Ax$.
 43: - N    - number of global rows (or `PETSC_DETERMINE` to have calculated if `n` is given)

 45:   Output Parameter:
 46: . C - the matrix

 48:   Notes:
 49:   The entries of the matrix `C` are not explicitly stored. Instead, the new matrix
 50:   object is a shell that simply performs the action of the centering matrix, i.e.,
 51:   multiplying $ C*x$ subtracts the mean of the vector `x` from each of its elements.
 52:   This is useful for preserving sparsity when mean-centering the columns of a
 53:   matrix is required. For instance, to perform principal components analysis with
 54:   a matrix `A`, the composite matrix $C*A$ can be passed to a partial SVD solver.

 56:   Level: intermediate

 58: .seealso: [](ch_matrices), `Mat`, `MatCreateLRC()`, `MatCreateComposite()`
 59: @*/
 60: PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C)
 61: {
 62:   PetscMPIInt size;

 64:   PetscFunctionBegin;
 65:   PetscCall(MatCreate(comm, C));
 66:   PetscCall(MatSetSizes(*C, n, n, N, N));
 67:   PetscCallMPI(MPI_Comm_size(comm, &size));
 68:   PetscCall(PetscObjectChangeTypeName((PetscObject)*C, MATCENTERING));

 70:   (*C)->ops->mult                   = MatMult_Centering;
 71:   (*C)->assembled                   = PETSC_TRUE;
 72:   (*C)->preallocated                = PETSC_TRUE;
 73:   (*C)->symmetric                   = PETSC_BOOL3_TRUE;
 74:   (*C)->symmetry_eternal            = PETSC_TRUE;
 75:   (*C)->structural_symmetry_eternal = PETSC_TRUE;
 76:   PetscCall(MatSetUp(*C));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }