1c1fff1c6SRichard Tran Mills 2c1fff1c6SRichard Tran Mills #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3c1fff1c6SRichard Tran Mills 4c1fff1c6SRichard Tran Mills PetscErrorCode MatMult_Centering(Mat A,Vec xx,Vec yy) 5c1fff1c6SRichard Tran Mills { 6c1fff1c6SRichard Tran Mills PetscScalar *y; 7c1fff1c6SRichard Tran Mills const PetscScalar *x; 8c1fff1c6SRichard Tran Mills PetscScalar sum,mean; 9c1fff1c6SRichard Tran Mills PetscInt i,m=A->rmap->n,size; 10c1fff1c6SRichard Tran Mills 11c1fff1c6SRichard Tran Mills PetscFunctionBegin; 129566063dSJacob Faibussowitsch PetscCall(VecSum(xx,&sum)); 139566063dSJacob Faibussowitsch PetscCall(VecGetSize(xx,&size)); 14c1fff1c6SRichard Tran Mills mean = sum / (PetscScalar)size; 159566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(xx,&x)); 169566063dSJacob Faibussowitsch PetscCall(VecGetArray(yy,&y)); 17c1fff1c6SRichard Tran Mills for (i=0; i<m; i++) { 18c1fff1c6SRichard Tran Mills y[i] = x[i] - mean; 19c1fff1c6SRichard Tran Mills } 209566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(xx,&x)); 219566063dSJacob Faibussowitsch PetscCall(VecRestoreArray(yy,&y)); 22c1fff1c6SRichard Tran Mills PetscFunctionReturn(0); 23c1fff1c6SRichard Tran Mills } 24c1fff1c6SRichard Tran Mills 25c1fff1c6SRichard Tran Mills /*@ 26c1fff1c6SRichard Tran Mills MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones' 27c1fff1c6SRichard Tran Mills 28c1fff1c6SRichard Tran Mills Collective on Mat 29c1fff1c6SRichard Tran Mills 30c1fff1c6SRichard Tran Mills Input Parameters: 31c1fff1c6SRichard Tran Mills + comm - MPI communicator 32c1fff1c6SRichard Tran Mills . n - number of local rows (or PETSC_DECIDE to have calculated if N is given) 33c1fff1c6SRichard Tran Mills This value should be the same as the local size used in creating the 34c1fff1c6SRichard Tran Mills y vector for the matrix-vector product y = Ax. 35c1fff1c6SRichard Tran Mills - N - number of global rows (or PETSC_DETERMINE to have calculated if n is given) 36c1fff1c6SRichard Tran Mills 37c1fff1c6SRichard Tran Mills Output Parameter: 38c1fff1c6SRichard Tran Mills . C - the matrix 39c1fff1c6SRichard Tran Mills 40c1fff1c6SRichard Tran Mills Notes: 41c1fff1c6SRichard Tran Mills The entries of the matrix C are not explicitly stored. Instead, the new matrix 42c1fff1c6SRichard Tran Mills object is a shell that simply performs the action of the centering matrix, i.e., 43c1fff1c6SRichard Tran Mills multiplying C*x subtracts the mean of the vector x from each of its elements. 44c1fff1c6SRichard Tran Mills This is useful for preserving sparsity when mean-centering the columns of a 45c1fff1c6SRichard Tran Mills matrix is required. For instance, to perform principal components analysis with 46c1fff1c6SRichard Tran Mills a matrix A, the composite matrix C*A can be passed to a partial SVD solver. 47c1fff1c6SRichard Tran Mills 48c1fff1c6SRichard Tran Mills Level: intermediate 49c1fff1c6SRichard Tran Mills 50db781477SPatrick Sanan .seealso: `MatCreateLRC()`, `MatCreateComposite()` 51c1fff1c6SRichard Tran Mills @*/ 52c1fff1c6SRichard Tran Mills PetscErrorCode MatCreateCentering(MPI_Comm comm,PetscInt n,PetscInt N,Mat *C) 53c1fff1c6SRichard Tran Mills { 54c1fff1c6SRichard Tran Mills PetscMPIInt size; 55c1fff1c6SRichard Tran Mills 56c1fff1c6SRichard Tran Mills PetscFunctionBegin; 579566063dSJacob Faibussowitsch PetscCall(MatCreate(comm,C)); 589566063dSJacob Faibussowitsch PetscCall(MatSetSizes(*C,n,n,N,N)); 599566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(comm,&size)); 609566063dSJacob Faibussowitsch PetscCall(PetscObjectChangeTypeName((PetscObject)*C,MATCENTERING)); 61c1fff1c6SRichard Tran Mills 62c1fff1c6SRichard Tran Mills (*C)->ops->mult = MatMult_Centering; 63c1fff1c6SRichard Tran Mills (*C)->assembled = PETSC_TRUE; 64c1fff1c6SRichard Tran Mills (*C)->preallocated = PETSC_TRUE; 65*b94d7dedSBarry Smith (*C)->symmetric = PETSC_BOOL3_TRUE; 66*b94d7dedSBarry Smith (*C)->symmetry_eternal = PETSC_TRUE; 67*b94d7dedSBarry Smith (*C)->structural_symmetry_eternal = PETSC_TRUE; 689566063dSJacob Faibussowitsch PetscCall(MatSetUp(*C)); 69c1fff1c6SRichard Tran Mills PetscFunctionReturn(0); 70c1fff1c6SRichard Tran Mills } 71