Actual source code: ex38.c
2: static char help[] = "Test interface of Elemental. \n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Mat C, Caij;
9: PetscInt i, j, m = 5, n, nrows, ncols;
10: const PetscInt *rows, *cols;
11: IS isrows, iscols;
12: PetscBool flg, Test_MatMatMult = PETSC_FALSE, mats_view = PETSC_FALSE;
13: PetscScalar *v;
14: PetscMPIInt rank, size;
16: PetscFunctionBeginUser;
17: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
18: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
19: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
20: PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
22: /* Get local block or element size*/
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
24: n = m;
25: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
27: PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
28: PetscCall(MatSetSizes(C, m, n, PETSC_DECIDE, PETSC_DECIDE));
29: PetscCall(MatSetType(C, MATELEMENTAL));
30: PetscCall(MatSetFromOptions(C));
31: PetscCall(MatSetUp(C));
33: PetscCall(PetscOptionsHasName(NULL, NULL, "-row_oriented", &flg));
34: if (flg) PetscCall(MatSetOption(C, MAT_ROW_ORIENTED, PETSC_TRUE));
35: PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
36: PetscCall(PetscOptionsHasName(NULL, NULL, "-Cexp_view_ownership", &flg));
37: if (flg) { /* View ownership of explicit C */
38: IS tmp;
39: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Ownership of explicit C:\n"));
40: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Row index set:\n"));
41: PetscCall(ISOnComm(isrows, PETSC_COMM_WORLD, PETSC_USE_POINTER, &tmp));
42: PetscCall(ISView(tmp, PETSC_VIEWER_STDOUT_WORLD));
43: PetscCall(ISDestroy(&tmp));
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Column index set:\n"));
45: PetscCall(ISOnComm(iscols, PETSC_COMM_WORLD, PETSC_USE_POINTER, &tmp));
46: PetscCall(ISView(tmp, PETSC_VIEWER_STDOUT_WORLD));
47: PetscCall(ISDestroy(&tmp));
48: }
50: /* Set local matrix entries */
51: PetscCall(ISGetLocalSize(isrows, &nrows));
52: PetscCall(ISGetIndices(isrows, &rows));
53: PetscCall(ISGetLocalSize(iscols, &ncols));
54: PetscCall(ISGetIndices(iscols, &cols));
55: PetscCall(PetscMalloc1(nrows * ncols, &v));
56: for (i = 0; i < nrows; i++) {
57: for (j = 0; j < ncols; j++) {
58: /*v[i*ncols+j] = (PetscReal)(rank);*/
59: v[i * ncols + j] = (PetscReal)(rank * 10000 + 100 * rows[i] + cols[j]);
60: }
61: }
62: PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
63: PetscCall(ISRestoreIndices(isrows, &rows));
64: PetscCall(ISRestoreIndices(iscols, &cols));
65: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
68: /* Test MatView() */
69: if (mats_view) PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
71: /* Test MatMissingDiagonal() */
72: PetscCall(MatMissingDiagonal(C, &flg, NULL));
73: PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "MatMissingDiagonal() did not return false!");
75: /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */
76: if (rank == 0) {
77: PetscInt M, N, cols[2];
78: PetscCall(MatGetSize(C, &M, &N));
79: for (i = 0; i < M; i++) {
80: cols[0] = i;
81: v[0] = i + 0.5;
82: cols[1] = i - 1;
83: v[1] = 0.5;
84: if (i) {
85: PetscCall(MatSetValues(C, 1, &i, 2, cols, v, ADD_VALUES));
86: } else {
87: PetscCall(MatSetValues(C, 1, &i, 1, &i, v, ADD_VALUES));
88: }
89: }
90: }
91: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
92: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
94: /* Test MatMult() */
95: PetscCall(MatComputeOperator(C, MATAIJ, &Caij));
96: PetscCall(MatMultEqual(C, Caij, 5, &flg));
97: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultEqual() fails");
98: PetscCall(MatMultTransposeEqual(C, Caij, 5, &flg));
99: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeEqual() fails");
101: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
102: PetscCall(MatMultAddEqual(C, Caij, 5, &flg));
103: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultAddEqual() fails");
104: PetscCall(MatMultTransposeAddEqual(C, Caij, 5, &flg));
105: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "C != Caij. MatMultTransposeAddEqual() fails");
107: /* Test MatMatMult() */
108: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_matmatmult", &Test_MatMatMult));
109: if (Test_MatMatMult) {
110: Mat CCelem, CCaij;
111: PetscCall(MatMatMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCelem));
112: PetscCall(MatMatMult(Caij, Caij, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CCaij));
113: PetscCall(MatMultEqual(CCelem, CCaij, 5, &flg));
114: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "CCelem != CCaij. MatMatMult() fails");
115: PetscCall(MatDestroy(&CCaij));
116: PetscCall(MatDestroy(&CCelem));
117: }
119: PetscCall(PetscFree(v));
120: PetscCall(ISDestroy(&isrows));
121: PetscCall(ISDestroy(&iscols));
122: PetscCall(MatDestroy(&C));
123: PetscCall(MatDestroy(&Caij));
124: PetscCall(PetscFinalize());
125: return 0;
126: }
128: /*TEST
130: test:
131: nsize: 2
132: requires: elemental
133: args: -mat_type elemental -m 2 -n 3
135: test:
136: suffix: 2
137: nsize: 6
138: requires: elemental
139: args: -mat_type elemental -m 2 -n 2
141: test:
142: suffix: 3
143: nsize: 6
144: requires: elemental
145: args: -mat_type elemental -m 2 -n 2 -test_matmatmult
146: output_file: output/ex38_2.out
148: TEST*/