1c1fff1c6SRichard Tran Mills static char help[] = "Tests MATCENTERING matrix type.\n\n";
2c1fff1c6SRichard Tran Mills
3c1fff1c6SRichard Tran Mills #include <petscmat.h>
4c1fff1c6SRichard Tran Mills
main(int argc,char ** argv)5d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
6d71ae5a4SJacob Faibussowitsch {
7c1fff1c6SRichard Tran Mills PetscInt n;
8c1fff1c6SRichard Tran Mills Mat C;
9c1fff1c6SRichard Tran Mills Vec x, y;
10c1fff1c6SRichard Tran Mills PetscReal norm;
11c1fff1c6SRichard Tran Mills PetscMPIInt size;
12c1fff1c6SRichard Tran Mills
13327415f7SBarry Smith PetscFunctionBeginUser;
14*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &argv, NULL, help));
159566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
16c1fff1c6SRichard Tran Mills
17c1fff1c6SRichard Tran Mills /* Create a parallel vector with 10*size total entries, and fill it with 1s. */
18c1fff1c6SRichard Tran Mills n = 10 * size;
199566063dSJacob Faibussowitsch PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
209566063dSJacob Faibussowitsch PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
219566063dSJacob Faibussowitsch PetscCall(VecSetFromOptions(x));
229566063dSJacob Faibussowitsch PetscCall(VecSet(x, 1.0));
23c1fff1c6SRichard Tran Mills
24c1fff1c6SRichard Tran Mills /* Create a corresponding n x n centering matrix and use it to create a mean-centered y = C * x. */
259566063dSJacob Faibussowitsch PetscCall(VecDuplicate(x, &y));
269566063dSJacob Faibussowitsch PetscCall(MatCreateCentering(PETSC_COMM_WORLD, PETSC_DECIDE, n, &C));
279566063dSJacob Faibussowitsch PetscCall(MatMult(C, x, y));
28c1fff1c6SRichard Tran Mills
29c1fff1c6SRichard Tran Mills /* Verify that the centered vector y has norm 0. */
309566063dSJacob Faibussowitsch PetscCall(VecNorm(y, NORM_2, &norm));
319566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Vector norm after MatMult() with centering matrix applied to vector of ones is %f.\n", (double)norm));
32c1fff1c6SRichard Tran Mills
33c1fff1c6SRichard Tran Mills /* Now repeat, but using MatMultTranspose(). */
349566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(C, x, y));
359566063dSJacob Faibussowitsch PetscCall(VecNorm(y, NORM_2, &norm));
369566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Vector norm after MatMultTranspose() with centering matrix applied to vector of ones is %f.\n", (double)norm));
37c1fff1c6SRichard Tran Mills
38c1fff1c6SRichard Tran Mills /* Clean up. */
399566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x));
409566063dSJacob Faibussowitsch PetscCall(VecDestroy(&y));
419566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C));
429566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
43b122ec5aSJacob Faibussowitsch return 0;
44c1fff1c6SRichard Tran Mills }
45c1fff1c6SRichard Tran Mills
46c1fff1c6SRichard Tran Mills /*TEST
47c1fff1c6SRichard Tran Mills
48c1fff1c6SRichard Tran Mills test:
49c1fff1c6SRichard Tran Mills suffix: 1
50c1fff1c6SRichard Tran Mills nsize: 1
51c1fff1c6SRichard Tran Mills output_file: output/ex247.out
52c1fff1c6SRichard Tran Mills
53c1fff1c6SRichard Tran Mills test:
54c1fff1c6SRichard Tran Mills suffix: 2
55c1fff1c6SRichard Tran Mills nsize: 2
56c1fff1c6SRichard Tran Mills output_file: output/ex247.out
57c1fff1c6SRichard Tran Mills
58c1fff1c6SRichard Tran Mills TEST*/
59