1 static char help[] = "Tests MATCENTERING matrix type.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **argv) 6 { 7 PetscInt n; 8 Mat C; 9 Vec x, y; 10 PetscReal norm; 11 PetscMPIInt size; 12 13 PetscFunctionBeginUser; 14 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 15 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 16 17 /* Create a parallel vector with 10*size total entries, and fill it with 1s. */ 18 n = 10 * size; 19 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 20 PetscCall(VecSetSizes(x, PETSC_DECIDE, n)); 21 PetscCall(VecSetFromOptions(x)); 22 PetscCall(VecSet(x, 1.0)); 23 24 /* Create a corresponding n x n centering matrix and use it to create a mean-centered y = C * x. */ 25 PetscCall(VecDuplicate(x, &y)); 26 PetscCall(MatCreateCentering(PETSC_COMM_WORLD, PETSC_DECIDE, n, &C)); 27 PetscCall(MatMult(C, x, y)); 28 29 /* Verify that the centered vector y has norm 0. */ 30 PetscCall(VecNorm(y, NORM_2, &norm)); 31 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Vector norm after MatMult() with centering matrix applied to vector of ones is %f.\n", (double)norm)); 32 33 /* Now repeat, but using MatMultTranspose(). */ 34 PetscCall(MatMultTranspose(C, x, y)); 35 PetscCall(VecNorm(y, NORM_2, &norm)); 36 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Vector norm after MatMultTranspose() with centering matrix applied to vector of ones is %f.\n", (double)norm)); 37 38 /* Clean up. */ 39 PetscCall(VecDestroy(&x)); 40 PetscCall(VecDestroy(&y)); 41 PetscCall(MatDestroy(&C)); 42 PetscCall(PetscFinalize()); 43 return 0; 44 } 45 46 /*TEST 47 48 test: 49 suffix: 1 50 nsize: 1 51 output_file: output/ex247.out 52 53 test: 54 suffix: 2 55 nsize: 2 56 output_file: output/ex247.out 57 58 TEST*/ 59