1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 PetscErrorCode MatMult_Centering(Mat A, Vec xx, Vec yy) { 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(0); 20 } 21 22 /*@ 23 MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones' 24 25 Collective on Mat 26 27 Input Parameters: 28 + comm - MPI communicator 29 . n - number of local rows (or PETSC_DECIDE to have calculated if N is given) 30 This value should be the same as the local size used in creating the 31 y vector for the matrix-vector product y = Ax. 32 - N - number of global rows (or PETSC_DETERMINE to have calculated if n is given) 33 34 Output Parameter: 35 . C - the matrix 36 37 Notes: 38 The entries of the matrix C are not explicitly stored. Instead, the new matrix 39 object is a shell that simply performs the action of the centering matrix, i.e., 40 multiplying C*x subtracts the mean of the vector x from each of its elements. 41 This is useful for preserving sparsity when mean-centering the columns of a 42 matrix is required. For instance, to perform principal components analysis with 43 a matrix A, the composite matrix C*A can be passed to a partial SVD solver. 44 45 Level: intermediate 46 47 .seealso: `MatCreateLRC()`, `MatCreateComposite()` 48 @*/ 49 PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C) { 50 PetscMPIInt size; 51 52 PetscFunctionBegin; 53 PetscCall(MatCreate(comm, C)); 54 PetscCall(MatSetSizes(*C, n, n, N, N)); 55 PetscCallMPI(MPI_Comm_size(comm, &size)); 56 PetscCall(PetscObjectChangeTypeName((PetscObject)*C, MATCENTERING)); 57 58 (*C)->ops->mult = MatMult_Centering; 59 (*C)->assembled = PETSC_TRUE; 60 (*C)->preallocated = PETSC_TRUE; 61 (*C)->symmetric = PETSC_BOOL3_TRUE; 62 (*C)->symmetry_eternal = PETSC_TRUE; 63 (*C)->structural_symmetry_eternal = PETSC_TRUE; 64 PetscCall(MatSetUp(*C)); 65 PetscFunctionReturn(0); 66 } 67