Actual source code: centering.c
2: #include <petsc/private/matimpl.h>
4: PetscErrorCode MatMult_Centering(Mat A,Vec xx,Vec yy)
5: {
6: PetscErrorCode ierr;
7: PetscScalar *y;
8: const PetscScalar *x;
9: PetscScalar sum,mean;
10: PetscInt i,m=A->rmap->n,size;
13: VecSum(xx,&sum);
14: VecGetSize(xx,&size);
15: mean = sum / (PetscScalar)size;
16: VecGetArrayRead(xx,&x);
17: VecGetArray(yy,&y);
18: for (i=0; i<m; i++) {
19: y[i] = x[i] - mean;
20: }
21: VecRestoreArrayRead(xx,&x);
22: VecRestoreArray(yy,&y);
23: return(0);
24: }
26: /*@
27: MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones'
29: Collective on Mat
31: Input Parameters:
32: + comm - MPI communicator
33: . n - number of local rows (or PETSC_DECIDE to have calculated if N is given)
34: This value should be the same as the local size used in creating the
35: y vector for the matrix-vector product y = Ax.
36: - N - number of global rows (or PETSC_DETERMINE to have calculated if n is given)
38: Output Parameter:
39: . C - the matrix
41: Notes:
42: The entries of the matrix C are not explicitly stored. Instead, the new matrix
43: object is a shell that simply performs the action of the centering matrix, i.e.,
44: multiplying C*x subtracts the mean of the vector x from each of its elements.
45: This is useful for preserving sparsity when mean-centering the columns of a
46: matrix is required. For instance, to perform principal components analysis with
47: a matrix A, the composite matrix C*A can be passed to a partial SVD solver.
49: Level: intermediate
51: .seealso: MatCreateLRC(), MatCreateComposite()
52: @*/
53: PetscErrorCode MatCreateCentering(MPI_Comm comm,PetscInt n,PetscInt N,Mat *C)
54: {
56: PetscMPIInt size;
59: MatCreate(comm,C);
60: MatSetSizes(*C,n,n,N,N);
61: MPI_Comm_size(comm,&size);
62: PetscObjectChangeTypeName((PetscObject)*C,MATCENTERING);
64: (*C)->ops->mult = MatMult_Centering;
65: (*C)->assembled = PETSC_TRUE;
66: (*C)->preallocated = PETSC_TRUE;
67: (*C)->symmetric = PETSC_TRUE;
68: (*C)->symmetric_eternal = PETSC_TRUE;
69: MatSetUp(*C);
70: return(0);
71: }