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