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, (char *)0, 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) {
93: PetscCall(MatDiagonalScale(C, NULL, d));
94: } else {
95: PetscCall(MatDiagonalScale(C, d, NULL));
96: }
97: if (mats_view) {
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n"));
99: PetscCall(MatView(C, PETSC_VIEWER_STDOUT_WORLD));
100: }
102: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
103: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
104: PetscCall(MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE));
105: PetscCall(MatSetType(B, MATELEMENTAL));
106: PetscCall(MatSetFromOptions(B));
107: PetscCall(MatSetUp(B));
108: PetscCall(MatGetOwnershipIS(B, &isrows, &iscols));
109: PetscCall(ISGetLocalSize(isrows, &nrows));
110: PetscCall(ISGetIndices(isrows, &rows));
111: PetscCall(ISGetLocalSize(iscols, &ncols));
112: PetscCall(ISGetIndices(iscols, &cols));
113: for (i = 0; i < nrows; i++) {
114: for (j = 0; j < ncols; j++) v[i * ncols + j] = (PetscReal)(1000 * rows[i] + cols[j]);
115: }
116: PetscCall(MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES));
117: PetscCall(PetscFree(v));
118: PetscCall(ISRestoreIndices(isrows, &rows));
119: PetscCall(ISRestoreIndices(iscols, &cols));
120: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
121: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
122: PetscCall(MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN));
123: PetscCall(MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN));
124: PetscCall(MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B));
125: if (mats_view) {
126: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n"));
127: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
128: }
129: PetscCall(ISDestroy(&isrows));
130: PetscCall(ISDestroy(&iscols));
131: PetscCall(MatDestroy(&B));
133: /* Test MatMatTransposeMult(): B = C*C^T */
134: PetscCall(MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B));
135: PetscCall(MatScale(C, 2.0));
136: PetscCall(MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B));
137: PetscCall(MatMatTransposeMultEqual(C, C, B, 10, &flg));
138: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "MatMatTransposeMult: B != C*B^T");
140: if (mats_view) {
141: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n"));
142: PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD));
143: }
145: PetscCall(MatDestroy(&Cdense));
146: PetscCall(PetscFree(v));
147: PetscCall(MatDestroy(&B));
148: PetscCall(MatDestroy(&C));
149: PetscCall(MatDestroy(&Ct));
150: PetscCall(VecDestroy(&d));
151: PetscCall(PetscFinalize());
152: return 0;
153: }
155: /*TEST
157: test:
158: nsize: 2
159: args: -m 3 -n 2
160: requires: elemental
162: test:
163: suffix: 2
164: nsize: 6
165: args: -m 2 -n 3
166: requires: elemental
168: test:
169: suffix: 3
170: nsize: 1
171: args: -m 2 -n 3
172: requires: elemental
173: output_file: output/ex39_1.out
175: TEST*/