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: }