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