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;
19: PetscInitialize(&argc, &args, (char *)0, help);
20: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
21: MPI_Comm_size(PETSC_COMM_WORLD, &size);
22: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
23: n = m;
24: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
25: PetscOptionsHasName(NULL, NULL, "-mats_view", &mats_view);
27: MatCreate(PETSC_COMM_WORLD, &C);
28: MatSetSizes(C, m, n, PETSC_DECIDE, PETSC_DECIDE);
29: MatSetType(C, MATELEMENTAL);
30: MatSetFromOptions(C);
31: MatSetUp(C);
32: MatGetOwnershipIS(C, &isrows, &iscols);
33: ISGetLocalSize(isrows, &nrows);
34: ISGetIndices(isrows, &rows);
35: ISGetLocalSize(iscols, &ncols);
36: ISGetIndices(iscols, &cols);
37: PetscMalloc1(nrows * ncols, &v);
38: #if defined(PETSC_USE_COMPLEX)
39: PetscRandom rand;
40: PetscScalar rval;
41: PetscRandomCreate(PETSC_COMM_WORLD, &rand);
42: PetscRandomSetFromOptions(rand);
43: for (i = 0; i < nrows; i++) {
44: for (j = 0; j < ncols; j++) {
45: PetscRandomGetValue(rand, &rval);
46: v[i * ncols + j] = rval;
47: }
48: }
49: 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: MatSetValues(C, nrows, rows, ncols, cols, v, INSERT_VALUES);
56: ISRestoreIndices(isrows, &rows);
57: ISRestoreIndices(iscols, &cols);
58: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
60: ISDestroy(&isrows);
61: ISDestroy(&iscols);
63: /* Test MatView(), MatDuplicate() and out-of-place MatConvert() */
64: MatDuplicate(C, MAT_COPY_VALUES, &B);
65: if (mats_view) {
66: PetscPrintf(PETSC_COMM_WORLD, "Duplicated C:\n");
67: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
68: }
69: MatDestroy(&B);
70: MatConvert(C, MATMPIDENSE, MAT_INITIAL_MATRIX, &Cdense);
71: MatMultEqual(C, Cdense, 5, &flg);
74: /* Test MatNorm() */
75: MatNorm(C, NORM_1, &Cnorm);
77: /* Test MatTranspose(), MatZeroEntries() and MatGetDiagonal() */
78: MatTranspose(C, MAT_INITIAL_MATRIX, &Ct);
79: MatConjugate(Ct);
80: if (mats_view) {
81: PetscPrintf(PETSC_COMM_WORLD, "C's Transpose Conjugate:\n");
82: MatView(Ct, PETSC_VIEWER_STDOUT_WORLD);
83: }
84: MatZeroEntries(Ct);
85: VecCreate(PETSC_COMM_WORLD, &d);
86: VecSetSizes(d, m > n ? n : m, PETSC_DECIDE);
87: VecSetFromOptions(d);
88: MatGetDiagonal(C, d);
89: if (mats_view) {
90: PetscPrintf(PETSC_COMM_WORLD, "Diagonal of C:\n");
91: VecView(d, PETSC_VIEWER_STDOUT_WORLD);
92: }
93: if (m > n) {
94: MatDiagonalScale(C, NULL, d);
95: } else {
96: MatDiagonalScale(C, d, NULL);
97: }
98: if (mats_view) {
99: PetscPrintf(PETSC_COMM_WORLD, "Diagonal Scaled C:\n");
100: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
101: }
103: /* Test MatAXPY(), MatAYPX() and in-place MatConvert() */
104: MatCreate(PETSC_COMM_WORLD, &B);
105: MatSetSizes(B, m, n, PETSC_DECIDE, PETSC_DECIDE);
106: MatSetType(B, MATELEMENTAL);
107: MatSetFromOptions(B);
108: MatSetUp(B);
109: MatGetOwnershipIS(B, &isrows, &iscols);
110: ISGetLocalSize(isrows, &nrows);
111: ISGetIndices(isrows, &rows);
112: ISGetLocalSize(iscols, &ncols);
113: 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: MatSetValues(B, nrows, rows, ncols, cols, v, INSERT_VALUES);
118: PetscFree(v);
119: ISRestoreIndices(isrows, &rows);
120: ISRestoreIndices(iscols, &cols);
121: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
122: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
123: MatAXPY(B, 2.5, C, SAME_NONZERO_PATTERN);
124: MatAYPX(B, 3.75, C, SAME_NONZERO_PATTERN);
125: MatConvert(B, MATDENSE, MAT_INPLACE_MATRIX, &B);
126: if (mats_view) {
127: PetscPrintf(PETSC_COMM_WORLD, "B after MatAXPY and MatAYPX:\n");
128: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
129: }
130: ISDestroy(&isrows);
131: ISDestroy(&iscols);
132: MatDestroy(&B);
134: /* Test MatMatTransposeMult(): B = C*C^T */
135: MatMatTransposeMult(C, C, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &B);
136: MatScale(C, 2.0);
137: MatMatTransposeMult(C, C, MAT_REUSE_MATRIX, PETSC_DEFAULT, &B);
138: MatMatTransposeMultEqual(C, C, B, 10, &flg);
141: if (mats_view) {
142: PetscPrintf(PETSC_COMM_WORLD, "C MatMatTransposeMult C:\n");
143: MatView(B, PETSC_VIEWER_STDOUT_WORLD);
144: }
146: MatDestroy(&Cdense);
147: PetscFree(v);
148: MatDestroy(&B);
149: MatDestroy(&C);
150: MatDestroy(&Ct);
151: VecDestroy(&d);
152: 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*/