1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2 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 @*/ 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