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