1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 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; 11 12 PetscFunctionBegin; 13 ierr = VecSum(xx,&sum);CHKERRQ(ierr); 14 ierr = VecGetSize(xx,&size);CHKERRQ(ierr); 15 mean = sum / (PetscScalar)size; 16 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 17 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 18 for (i=0; i<m; i++) { 19 y[i] = x[i] - mean; 20 } 21 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 22 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 /*@ 27 MatCreateCentering - Creates a new matrix object that implements the (symmetric and idempotent) centering matrix, I - (1/N) * ones*ones' 28 29 Collective on Mat 30 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) 37 38 Output Parameter: 39 . C - the matrix 40 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. 48 49 Level: intermediate 50 51 .seealso: MatCreateLRC(), MatCreateComposite() 52 @*/ 53 PetscErrorCode MatCreateCentering(MPI_Comm comm,PetscInt n,PetscInt N,Mat *C) 54 { 55 PetscErrorCode ierr; 56 PetscMPIInt size; 57 58 PetscFunctionBegin; 59 ierr = MatCreate(comm,C);CHKERRQ(ierr); 60 ierr = MatSetSizes(*C,n,n,N,N);CHKERRQ(ierr); 61 ierr = MPI_Comm_size(comm,&size);CHKERRMPI(ierr); 62 ierr = PetscObjectChangeTypeName((PetscObject)*C,MATCENTERING);CHKERRQ(ierr); 63 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 ierr = MatSetUp(*C);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72