1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
MatMult_Centering(Mat A,Vec xx,Vec yy)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;
9
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 }
21
22 /*MC
23 MATCENTERING - matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$
24
25 Level: beginner
26
27 Note:
28 See `MatCreateCentering()`
29
30 .seealso: `Mat`, `MatCreateCentering()`
31 M*/
32
33 /*@
34 MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$
35
36 Collective
37
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)
44
45 Output Parameter:
46 . C - the matrix
47
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.
55
56 Level: intermediate
57
58 .seealso: [](ch_matrices), `Mat`, `MatCreateLRC()`, `MatCreateComposite()`
59 @*/
MatCreateCentering(MPI_Comm comm,PetscInt n,PetscInt N,Mat * C)60 PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C)
61 {
62 PetscMPIInt size;
63
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));
69
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 (*C)->structurally_symmetric = PETSC_BOOL3_TRUE;
77 if (!PetscDefined(USE_COMPLEX)) (*C)->hermitian = PETSC_BOOL3_TRUE;
78 PetscCall(MatSetUp(*C));
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81