1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 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; 10 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 } 22 23 /*@ 24 MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones' 25 26 Collective 27 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) 34 35 Output Parameter: 36 . C - the matrix 37 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. 45 46 Level: intermediate 47 48 .seealso: `MatCreateLRC()`, `MatCreateComposite()` 49 @*/ 50 PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C) 51 { 52 PetscMPIInt size; 53 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)); 59 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 } 69