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: /*@
 23:   MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones'

 25:   Collective

 27:   Input Parameters:
 28: + comm - MPI communicator
 29: . n    - number of local rows (or `PETSC_DECIDE` to have calculated if `N` is given)
 30:            This value should be the same as the local size used in creating the
 31:            y vector for the matrix-vector product y = Ax.
 32: - N    - number of global rows (or `PETSC_DETERMINE` to have calculated if `n` is given)

 34:   Output Parameter:
 35: . C - the matrix

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

 45:   Level: intermediate

 47: .seealso: [](ch_matrices), `Mat`, `MatCreateLRC()`, `MatCreateComposite()`
 48: @*/
 49: PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C)
 50: {
 51:   PetscMPIInt size;

 53:   PetscFunctionBegin;
 54:   PetscCall(MatCreate(comm, C));
 55:   PetscCall(MatSetSizes(*C, n, n, N, N));
 56:   PetscCallMPI(MPI_Comm_size(comm, &size));
 57:   PetscCall(PetscObjectChangeTypeName((PetscObject)*C, MATCENTERING));

 59:   (*C)->ops->mult                   = MatMult_Centering;
 60:   (*C)->assembled                   = PETSC_TRUE;
 61:   (*C)->preallocated                = PETSC_TRUE;
 62:   (*C)->symmetric                   = PETSC_BOOL3_TRUE;
 63:   (*C)->symmetry_eternal            = PETSC_TRUE;
 64:   (*C)->structural_symmetry_eternal = PETSC_TRUE;
 65:   PetscCall(MatSetUp(*C));
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }