Actual source code: ex39.c
1: static char help[] = "Tests Elemental interface.\n\n";
3: #include <petscmat.h>
5: int main(int argc, char **args)
6: {
7: Mat Cdense, B, C, Ct;
8: Vec d;
9: PetscInt i, j, m = 5, n, nrows, ncols;
10: const PetscInt *rows, *cols;
11: IS isrows, iscols;
12: PetscScalar *v;
13: PetscMPIInt rank, size;
14: PetscReal Cnorm;
15: PetscBool flg, mats_view = PETSC_FALSE;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, NULL, help));
19: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
20: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
21: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
22: n = m;
23: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
24: PetscCall(PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view));
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));
31: PetscCall(MatGetOwnershipIS(C, &isrows, &iscols));
32: PetscCall(ISGetLocalSize(isrows, &nrows));
33: PetscCall(ISGetIndices(isrows, &rows));
34: PetscCall(ISGetLocalSize(iscols, &ncols));
35: PetscCall(ISGetIndices(iscols, &cols));
36: PetscCall(PetscMalloc1(nrows * ncols, &v));
37: #if defined(PETSC_USE_COMPLEX)
38: PetscRandom rand;
39: PetscScalar rval;
40: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
41: PetscCall(PetscRandomSetFromOptions(rand));
42: for (i = 0; i < nrows; i++) {
43: for (j = 0; j < ncols; j++) {
44: PetscCall(PetscRandomGetValue(rand, &rval));
45: v[i * ncols + j] = rval;
46: }
47: }
48: PetscCall(PetscRandomDestroy(&rand));
49: #else
50: for (i = 0; i < nrows; i++) {
51: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(10000 * rank + 100 * rows[i] + cols[j]);
52: }
53: #endif
54: PetscCall(MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES));
55: PetscCall(ISRestoreIndices(isrows, &rows));
56: PetscCall(ISRestoreIndices(iscols, &cols));
57: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
58: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
59: PetscCall(ISDestroy(&isrows));
60: PetscCall(ISDestroy(&iscols));
62: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
63: PetscCall(MatDuplicate(C, MAT_COPY_VALUES, &B));
64: if (mats_view) {
65: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n"));
66: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
67: }
68: PetscCall(MatDestroy(&B));
69: PetscCall(MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdense));
70: PetscCall(MatMultEqual(C, Cdense, 5, &flg));
71: PetscCheck(flg, PETSC_COMM_SELF, PETSC_ERR_ARG_NOTSAMETYPE, "Cdense != C. MatConvert() fails");
73: /* Test MatNorm() */
74: PetscCall(MatNorm(C, NORM_1, &Cnorm));
76: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
77: PetscCall(MatTranspose(C, MAT_INITIAL_MATRIX, &Ct));
78: PetscCall(MatConjugate(Ct));
79: if (mats_view) {
80: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n"));
81: PetscCall(MatView(Ct, PETSC_VIEWER_STDOUT_WORLD));
82: }
83: PetscCall(MatZeroEntries(Ct));
84: PetscCall(VecCreate(PETSC_COMM_WORLD, &d));
85: PetscCall(VecSetSizes(d, m > n ? n : m, PETSC_DECIDE));
86: PetscCall(VecSetFromOptions(d));
87: PetscCall(MatGetDiagonal(C, d));
88: if (mats_view) {
89: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n"));
90: PetscCall(VecView(d, PETSC_VIEWER_STDOUT_WORLD));
91: }
92: if (m > n) PetscCall(MatDiagonalScale(C, NULL, d));
93: else PetscCall(MatDiagonalScale(C, d, NULL));
94: if (mats_view) {
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
96: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
97: }
99: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
100: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
101: PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
102: PetscCall(MatSetType(B, MATELEMENTAL));
103: PetscCall(MatSetFromOptions(B));
104: PetscCall(MatSetUp(B));
105: PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
106: PetscCall(ISGetLocalSize(isrows, &nrows));
107: PetscCall(ISGetIndices(isrows, &rows));
108: PetscCall(ISGetLocalSize(iscols, &ncols));
109: PetscCall(ISGetIndices(iscols, &cols));
110: for (i = 0; i < nrows; i++) {
111: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
112: }
113: PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
114: PetscCall(PetscFree(v));
115: PetscCall(ISRestoreIndices(isrows, &rows));
116: PetscCall(ISRestoreIndices(iscols, &cols));
117: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
118: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
119: PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
120: PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
121: PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
122: if (mats_view) {
123: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
124: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
125: }
126: PetscCall(ISDestroy(&isrows));
127: PetscCall(ISDestroy(&iscols));
128: PetscCall(MatDestroy(&B));
130: /* Test MatMatTransposeMult(): B = C*C^T */
131: PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &B));
132: PetscCall(MatScale(C, 2.0));
133: PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DETERMINE, &B));
134: PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
135: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTransposeMult: B != C*B^T");
137: if (mats_view) {
138: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
139: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
140: }
142: PetscCall(MatDestroy(&Cdense));
143: PetscCall(PetscFree(v));
144: PetscCall(MatDestroy(&B));
145: PetscCall(MatDestroy(&C));
146: PetscCall(MatDestroy(&Ct));
147: PetscCall(VecDestroy(&d));
148: PetscCall(PetscFinalize());
149: return 0;
150: }
152: /*TEST
154: test:
155: nsize: 2
156: args: -m 3 -n 2
157: requires: elemental
158: output_file: output/empty.out
160: test:
161: suffix: 2
162: nsize: 6
163: args: -m 2 -n 3
164: requires: elemental
165: output_file: output/empty.out
167: test:
168: suffix: 3
169: nsize: 1
170: args: -m 2 -n 3
171: requires: elemental
172: output_file: output/empty.out
174: TEST*/