Actual source code: ex222.c
1: static char help[] = "Tests MatComputeOperator() and MatComputeOperatorTranspose()\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **argv)
6: {
7: Mat A, Ae, Aet;
8: char filename[PETSC_MAX_PATH_LEN];
9: char expltype[128], *etype = NULL;
10: PetscInt bs = 1;
11: PetscBool flg, check = PETSC_TRUE;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
16: PetscCall(PetscOptionsGetString(NULL, NULL, "-expl_type", expltype, sizeof(expltype), &flg));
17: if (flg) PetscCall(PetscStrallocpy(expltype, &etype));
18: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
20: if (!flg) {
21: PetscInt M = 13, N = 6;
22: PetscBool sbaij;
24: PetscCall(PetscOptionsGetInt(NULL, NULL, "-M", &M, NULL));
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
26: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, M, N, NULL, &A));
27: PetscCall(MatSetBlockSize(A, bs));
28: PetscCall(MatSetRandom(A, NULL));
29: PetscCall(PetscStrcmp(etype, MATSBAIJ, &sbaij));
30: if (sbaij) {
31: Mat At;
32: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
33: PetscCall(MatAXPY(A, 1.0, At, DIFFERENT_NONZERO_PATTERN));
34: PetscCall(MatDestroy(&At));
35: }
36: } else {
37: PetscViewer viewer;
39: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
40: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
41: PetscCall(MatSetBlockSize(A, bs));
42: PetscCall(MatSetFromOptions(A));
43: PetscCall(MatLoad(A, viewer));
44: PetscCall(PetscViewerDestroy(&viewer));
45: }
46: PetscCall(PetscObjectSetName((PetscObject)A, "Matrix"));
47: PetscCall(MatViewFromOptions(A, NULL, "-view_expl"));
49: PetscCall(MatComputeOperator(A, etype, &Ae));
50: PetscCall(PetscObjectSetName((PetscObject)Ae, "Explicit matrix"));
51: PetscCall(MatViewFromOptions(Ae, NULL, "-view_expl"));
53: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check", &check, NULL));
54: if (check) {
55: Mat A2;
56: PetscReal err, tol = PETSC_SMALL;
58: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tol", &tol, NULL));
59: PetscCall(MatConvert(A, etype, MAT_INITIAL_MATRIX, &A2));
60: PetscCall(MatAXPY(A2, -1.0, Ae, DIFFERENT_NONZERO_PATTERN));
61: PetscCall(MatNorm(A2, NORM_FROBENIUS, &err));
62: if (err > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error %g > %g (type %s)\n", (double)err, (double)tol, etype));
63: PetscCall(MatDestroy(&A2));
64: }
66: PetscCall(MatComputeOperatorTranspose(A, etype, &Aet));
67: PetscCall(PetscObjectSetName((PetscObject)Aet, "Explicit matrix transpose"));
68: PetscCall(MatViewFromOptions(Aet, NULL, "-view_expl"));
70: PetscCall(PetscFree(etype));
71: PetscCall(MatDestroy(&Ae));
72: PetscCall(MatDestroy(&Aet));
73: PetscCall(MatDestroy(&A));
74: PetscCall(PetscFinalize());
75: return 0;
76: }
78: /*TEST
80: test:
81: output_file: output/ex222_null.out
83: testset:
84: suffix: matexpl_rect
85: output_file: output/ex222_null.out
86: nsize: {{1 3}}
87: args: -expl_type {{dense aij baij}}
89: testset:
90: suffix: matexpl_square
91: output_file: output/ex222_null.out
92: nsize: {{1 3}}
93: args: -bs {{1 2 3}} -M 36 -N 36 -expl_type {{dense aij baij sbaij}}
95: TEST*/