1c1fff1c6SRichard Tran Mills #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 2c1fff1c6SRichard Tran Mills 366976f2fSJacob Faibussowitsch static PetscErrorCode MatMult_Centering(Mat A, Vec xx, Vec yy) 4d71ae5a4SJacob Faibussowitsch { 5c1fff1c6SRichard Tran Mills PetscScalar *y; 6c1fff1c6SRichard Tran Mills const PetscScalar *x; 7c1fff1c6SRichard Tran Mills PetscScalar sum, mean; 8c1fff1c6SRichard Tran Mills PetscInt i, m = A->rmap->n, size; 9c1fff1c6SRichard Tran Mills 10c1fff1c6SRichard Tran Mills PetscFunctionBegin; 119566063dSJacob Faibussowitsch PetscCall(VecSum(xx, &sum)); 129566063dSJacob Faibussowitsch PetscCall(VecGetSize(xx, &size)); 13c1fff1c6SRichard Tran Mills mean = sum / (PetscScalar)size; 149566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(xx, &x)); 159566063dSJacob Faibussowitsch PetscCall(VecGetArray(yy, &y)); 16ad540459SPierre Jolivet for (i = 0; i < m; i++) y[i] = x[i] - mean; 179566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(xx, &x)); 189566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(yy, &y)); 193ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 20c1fff1c6SRichard Tran Mills } 21c1fff1c6SRichard Tran Mills 22*7f296bb3SBarry Smith /*MC 23*7f296bb3SBarry Smith MATCENTERING - matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$ 24*7f296bb3SBarry Smith 25*7f296bb3SBarry Smith Level: beginner 26*7f296bb3SBarry Smith 27*7f296bb3SBarry Smith Note: 28*7f296bb3SBarry Smith See `MatCreateCentering()` 29*7f296bb3SBarry Smith 30*7f296bb3SBarry Smith .seealso: `Mat`, `MatCreateCentering()` 31*7f296bb3SBarry Smith M*/ 32*7f296bb3SBarry Smith 33c1fff1c6SRichard Tran Mills /*@ 34*7f296bb3SBarry Smith MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, $ I - (1/N) * 1^T * 1$ 35c1fff1c6SRichard Tran Mills 3611a5261eSBarry Smith Collective 37c1fff1c6SRichard Tran Mills 38c1fff1c6SRichard Tran Mills Input Parameters: 39c1fff1c6SRichard Tran Mills + comm - MPI communicator 402ef1f0ffSBarry Smith . n - number of local rows (or `PETSC_DECIDE` to have calculated if `N` is given) 41c1fff1c6SRichard Tran Mills This value should be the same as the local size used in creating the 42*7f296bb3SBarry Smith `y` vector for the matrix-vector product $y = Ax$. 432ef1f0ffSBarry Smith - N - number of global rows (or `PETSC_DETERMINE` to have calculated if `n` is given) 44c1fff1c6SRichard Tran Mills 45c1fff1c6SRichard Tran Mills Output Parameter: 46c1fff1c6SRichard Tran Mills . C - the matrix 47c1fff1c6SRichard Tran Mills 48c1fff1c6SRichard Tran Mills Notes: 492ef1f0ffSBarry Smith The entries of the matrix `C` are not explicitly stored. Instead, the new matrix 50c1fff1c6SRichard Tran Mills object is a shell that simply performs the action of the centering matrix, i.e., 51*7f296bb3SBarry Smith multiplying $ C*x$ subtracts the mean of the vector `x` from each of its elements. 52c1fff1c6SRichard Tran Mills This is useful for preserving sparsity when mean-centering the columns of a 53c1fff1c6SRichard Tran Mills matrix is required. For instance, to perform principal components analysis with 54*7f296bb3SBarry Smith a matrix `A`, the composite matrix $C*A$ can be passed to a partial SVD solver. 55c1fff1c6SRichard Tran Mills 56c1fff1c6SRichard Tran Mills Level: intermediate 57c1fff1c6SRichard Tran Mills 581cc06b55SBarry Smith .seealso: [](ch_matrices), `Mat`, `MatCreateLRC()`, `MatCreateComposite()` 59c1fff1c6SRichard Tran Mills @*/ 60d71ae5a4SJacob Faibussowitsch PetscErrorCode MatCreateCentering(MPI_Comm comm, PetscInt n, PetscInt N, Mat *C) 61d71ae5a4SJacob Faibussowitsch { 62c1fff1c6SRichard Tran Mills PetscMPIInt size; 63c1fff1c6SRichard Tran Mills 64c1fff1c6SRichard Tran Mills PetscFunctionBegin; 659566063dSJacob Faibussowitsch PetscCall(MatCreate(comm, C)); 669566063dSJacob Faibussowitsch PetscCall(MatSetSizes(*C, n, n, N, N)); 679566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm, &size)); 689566063dSJacob Faibussowitsch PetscCall(PetscObjectChangeTypeName((PetscObject)*C, MATCENTERING)); 69c1fff1c6SRichard Tran Mills 70c1fff1c6SRichard Tran Mills (*C)->ops->mult = MatMult_Centering; 71c1fff1c6SRichard Tran Mills (*C)->assembled = PETSC_TRUE; 72c1fff1c6SRichard Tran Mills (*C)->preallocated = PETSC_TRUE; 73b94d7dedSBarry Smith (*C)->symmetric = PETSC_BOOL3_TRUE; 74b94d7dedSBarry Smith (*C)->symmetry_eternal = PETSC_TRUE; 75b94d7dedSBarry Smith (*C)->structural_symmetry_eternal = PETSC_TRUE; 769566063dSJacob Faibussowitsch PetscCall(MatSetUp(*C)); 773ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 78c1fff1c6SRichard Tran Mills } 79