Actual source code: ex19.c
2: static char help[] = "Tests reusing MPI parallel matrices and MatGetValues().\n\
3: To test the parallel matrix assembly, this example intentionally lays out\n\
4: the matrix across processors differently from the way it is assembled.\n\
5: This example uses bilinear elements on the unit square. Input arguments are:\n\
6: -m <size> : problem size\n\n";
8: #include <petscmat.h>
10: int FormElementStiffness(PetscReal H, PetscScalar *Ke)
11: {
12: Ke[0] = H / 6.0;
13: Ke[1] = -.125 * H;
14: Ke[2] = H / 12.0;
15: Ke[3] = -.125 * H;
16: Ke[4] = -.125 * H;
17: Ke[5] = H / 6.0;
18: Ke[6] = -.125 * H;
19: Ke[7] = H / 12.0;
20: Ke[8] = H / 12.0;
21: Ke[9] = -.125 * H;
22: Ke[10] = H / 6.0;
23: Ke[11] = -.125 * H;
24: Ke[12] = -.125 * H;
25: Ke[13] = H / 12.0;
26: Ke[14] = -.125 * H;
27: Ke[15] = H / 6.0;
28: return 0;
29: }
31: int main(int argc, char **args)
32: {
33: Mat C;
34: Vec u, b;
35: PetscMPIInt size, rank;
36: PetscInt i, m = 5, N, start, end, M, idx[4];
37: PetscInt j, nrsub, ncsub, *rsub, *csub, mystart, myend;
38: PetscBool flg;
39: PetscScalar one = 1.0, Ke[16], *vals;
40: PetscReal h, norm;
43: PetscInitialize(&argc, &args, (char *)0, help);
44: PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL);
46: N = (m + 1) * (m + 1); /* dimension of matrix */
47: M = m * m; /* number of elements */
48: h = 1.0 / m; /* mesh width */
49: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
50: MPI_Comm_size(PETSC_COMM_WORLD, &size);
52: /* Create stiffness matrix */
53: MatCreate(PETSC_COMM_WORLD, &C);
54: MatSetSizes(C, PETSC_DECIDE, PETSC_DECIDE, N, N);
55: MatSetFromOptions(C);
56: MatSetUp(C);
58: start = rank * (M / size) + ((M % size) < rank ? (M % size) : rank);
59: end = start + M / size + ((M % size) > rank);
61: /* Form the element stiffness for the Laplacian */
62: FormElementStiffness(h * h, Ke);
63: for (i = start; i < end; i++) {
64: /* location of lower left corner of element */
65: /* node numbers for the four corners of element */
66: idx[0] = (m + 1) * (i / m) + (i % m);
67: idx[1] = idx[0] + 1;
68: idx[2] = idx[1] + m + 1;
69: idx[3] = idx[2] - 1;
70: MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
71: }
72: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
75: /* Assemble the matrix again */
76: MatZeroEntries(C);
78: for (i = start; i < end; i++) {
79: /* location of lower left corner of element */
80: /* node numbers for the four corners of element */
81: idx[0] = (m + 1) * (i / m) + (i % m);
82: idx[1] = idx[0] + 1;
83: idx[2] = idx[1] + m + 1;
84: idx[3] = idx[2] - 1;
85: MatSetValues(C, 4, idx, 4, idx, Ke, ADD_VALUES);
86: }
87: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
90: /* Create test vectors */
91: VecCreate(PETSC_COMM_WORLD, &u);
92: VecSetSizes(u, PETSC_DECIDE, N);
93: VecSetFromOptions(u);
94: VecDuplicate(u, &b);
95: VecSet(u, one);
97: /* Check error */
98: MatMult(C, u, b);
99: VecNorm(b, NORM_2, &norm);
100: if (norm > PETSC_SQRT_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Norm of error b %g should be near 0\n", (double)norm);
102: /* Now test MatGetValues() */
103: PetscOptionsHasName(NULL, NULL, "-get_values", &flg);
104: if (flg) {
105: MatGetOwnershipRange(C, &mystart, &myend);
106: nrsub = myend - mystart;
107: ncsub = 4;
108: PetscMalloc1(nrsub * ncsub, &vals);
109: PetscMalloc1(nrsub, &rsub);
110: PetscMalloc1(ncsub, &csub);
111: for (i = myend - 1; i >= mystart; i--) rsub[myend - i - 1] = i;
112: for (i = 0; i < ncsub; i++) csub[i] = 2 * (ncsub - i) + mystart;
113: MatGetValues(C, nrsub, rsub, ncsub, csub, vals);
114: MatView(C, PETSC_VIEWER_STDOUT_WORLD);
115: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "processor number %d: start=%" PetscInt_FMT ", end=%" PetscInt_FMT ", mystart=%" PetscInt_FMT ", myend=%" PetscInt_FMT "\n", rank, start, end, mystart, myend);
116: for (i = 0; i < nrsub; i++) {
117: for (j = 0; j < ncsub; j++) {
118: if (PetscImaginaryPart(vals[i * ncsub + j]) != 0.0) {
119: PetscSynchronizedPrintf(PETSC_COMM_WORLD, " C[%" PetscInt_FMT ", %" PetscInt_FMT "] = %g + %g i\n", rsub[i], csub[j], (double)PetscRealPart(vals[i * ncsub + j]), (double)PetscImaginaryPart(vals[i * ncsub + j]));
120: } else {
121: PetscSynchronizedPrintf(PETSC_COMM_WORLD, " C[%" PetscInt_FMT ", %" PetscInt_FMT "] = %g\n", rsub[i], csub[j], (double)PetscRealPart(vals[i * ncsub + j]));
122: }
123: }
124: }
125: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
126: PetscFree(rsub);
127: PetscFree(csub);
128: PetscFree(vals);
129: }
131: /* Free data structures */
132: VecDestroy(&u);
133: VecDestroy(&b);
134: MatDestroy(&C);
135: PetscFinalize();
136: return 0;
137: }
139: /*TEST
141: test:
142: nsize: 4
144: TEST*/