Actual source code: centering.c


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

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

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

 23: /*@
 24:    MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones'

 26:    Collective

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

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

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

 46:    Level: intermediate

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

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

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