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