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