Actual source code: ex90.c
1: static char help[] = "Solve multiple shifted linear systems.\n\nInput arguments are:\n\
2: -m size - problem size\n\
3: -nshift nshift - number of shifts\n\
4: -mass (true|false) - test with a mass matrix M\n\
5: -explicitmat (true|false) - build the nested matrix explicitly\n\
6: -cmplx (true|false) - test with complex shifts\n\n";
8: #include <petscksp.h>
10: /* element stiffness for Laplacian */
11: PetscErrorCode FormElementStiffness(PetscReal H, PetscScalar *Ke)
12: {
13: PetscFunctionBeginUser;
14: Ke[0] = H / 6.0;
15: Ke[1] = -.125 * H;
16: Ke[2] = H / 12.0;
17: Ke[3] = -.125 * H;
18: Ke[4] = -.125 * H;
19: Ke[5] = H / 6.0;
20: Ke[6] = -.125 * H;
21: Ke[7] = H / 12.0;
22: Ke[8] = H / 12.0;
23: Ke[9] = -.125 * H;
24: Ke[10] = H / 6.0;
25: Ke[11] = -.125 * H;
26: Ke[12] = -.125 * H;
27: Ke[13] = H / 12.0;
28: Ke[14] = -.125 * H;
29: Ke[15] = H / 6.0;
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: PetscErrorCode FormElementRhs(PetscScalar x, PetscScalar y, PetscReal H, PetscScalar *r)
34: {
35: PetscFunctionBeginUser;
36: r[0] = 0.;
37: r[1] = 0.;
38: r[2] = 0.;
39: r[3] = 0.0;
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /* compute the norm ||(K+alpha*M)*x+beta*M*y-b|| using two workspace vectors r and w */
44: PetscErrorCode ComputeResidualNorm(Mat K, PetscScalar alpha, Mat M, Vec x, Vec b, Vec r, Vec w, PetscScalar beta, Vec y, PetscReal *norm)
45: {
46: PetscFunctionBeginUser;
47: PetscCall(MatMult(K, x, r));
48: if (M) {
49: PetscCall(MatMult(M, x, w));
50: PetscCall(VecAXPY(r, alpha, w));
51: } else PetscCall(VecAXPY(r, alpha, x));
52: if (y) {
53: if (M) {
54: PetscCall(MatMult(M, y, w));
55: PetscCall(VecAXPY(r, beta, w));
56: } else PetscCall(VecAXPY(r, beta, y));
57: }
58: if (b) PetscCall(VecAXPY(r, -1.0, b));
59: PetscCall(VecNorm(r, NORM_2, norm));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: /* print a message with the shift and the residual error */
64: PetscErrorCode PrintError(PetscScalar sigma, PetscScalar sigmai, PetscReal norm, PetscReal normb, PetscBool conjpair, PetscReal tol)
65: {
66: char msg[PETSC_MAX_PATH_LEN];
68: PetscFunctionBeginUser;
69: if (norm / normb < 10 * tol) PetscCall(PetscSNPrintf(msg, sizeof(msg), "Relative residual error < 10 * %g", (double)tol));
70: else PetscCall(PetscSNPrintf(msg, sizeof(msg), "Relative residual error %g", (double)(norm / normb)));
71: #if PetscDefined(USE_COMPLEX)
72: if (PetscImaginaryPart(sigma) == 0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma=%g: %s\n", (double)PetscRealPart(sigma), msg));
73: else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma=%g%+gi: %s\n", (double)PetscRealPart(sigma), (double)PetscImaginaryPart(sigma), msg));
74: #else
75: if (!conjpair) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma=%g: %s\n", (double)sigma, msg));
76: else {
77: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma=%g+%gi: %s\n", (double)sigma, (double)sigmai, msg));
78: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sigma=%g-%gi: %s\n", (double)sigma, (double)sigmai, msg));
79: }
80: #endif
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: int main(int argc, char **args)
85: {
86: Vec b, x, r, b_nest, x_nest, w = NULL;
87: Mat A, K, M = NULL;
88: KSP ksp;
89: PetscInt N, m2, idx[4], count, m = 80, nshift = 4, start, end, ostart, oend;
90: PetscInt *rows;
91: PetscReal h, norm, normb, rtol;
92: PetscScalar Ke[16], v[4], xx, yy, lend = -8e-4, rend = -2e-4;
93: PetscScalar *sigma = NULL, *sigma_imaginary = NULL;
94: PetscBool mass = PETSC_FALSE, explicitmat = PETSC_FALSE, cmplx = PETSC_FALSE;
95: PetscMPIInt rank, size;
96: #if !PetscDefined(USE_COMPLEX)
97: Vec xi;
98: PetscReal normi;
99: #endif
101: PetscFunctionBeginUser;
102: PetscCall(PetscInitialize(&argc, &args, NULL, help));
103: PetscCall(PetscOptionsGetBool(NULL, NULL, "-mass", &mass, NULL));
104: PetscCall(PetscOptionsGetBool(NULL, NULL, "-explicitmat", &explicitmat, NULL));
105: PetscCall(PetscOptionsGetBool(NULL, NULL, "-cmplx", &cmplx, NULL));
106: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
107: N = (m + 1) * (m + 1);
108: m2 = m * m;
109: h = 1.0 / m;
110: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nshift", &nshift, NULL));
111: PetscCheck(nshift > 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "nshift should be at least 1");
112: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
113: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Build the matrices and right-hand side
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Create stiffness matrix
121: */
122: PetscCall(MatCreate(PETSC_COMM_WORLD, &K));
123: PetscCall(MatSetSizes(K, PETSC_DECIDE, PETSC_DECIDE, N, N));
124: PetscCall(MatSetFromOptions(K));
125: PetscCall(MatSeqAIJSetPreallocation(K, 9, NULL));
126: PetscCall(MatMPIAIJSetPreallocation(K, 9, NULL, 8, NULL));
127: start = rank * (m2 / size) + ((m2 % size) < rank ? (m2 % size) : rank);
128: end = start + m2 / size + ((m2 % size) > rank);
130: PetscCall(FormElementStiffness(h * h, Ke));
131: for (PetscInt i = start; i < end; i++) {
132: /* node numbers for the four corners of element */
133: idx[0] = (m + 1) * (i / m) + (i % m);
134: idx[1] = idx[0] + 1;
135: idx[2] = idx[1] + m + 1;
136: idx[3] = idx[2] - 1;
137: PetscCall(MatSetValues(K, 4, idx, 4, idx, Ke, ADD_VALUES));
138: }
139: PetscCall(MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY));
140: PetscCall(MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY));
142: if (mass) {
143: /*
144: Create mass matrix
145: */
146: PetscCall(MatCreate(PETSC_COMM_WORLD, &M));
147: PetscCall(MatSetSizes(M, PETSC_DECIDE, PETSC_DECIDE, N, N));
148: PetscCall(MatSetFromOptions(M));
149: PetscCall(MatGetOwnershipRange(M, &ostart, &oend));
150: for (PetscInt i = ostart; i < oend; i++) PetscCall(MatSetValue(M, i, i, 3.0, INSERT_VALUES));
151: PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
152: PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
153: }
155: /*
156: Assemble right-hand side
157: */
158: PetscCall(MatCreateVecs(K, &r, &b));
159: if (M) PetscCall(VecDuplicate(b, &w));
160: for (PetscInt i = start; i < end; i++) {
161: /* location of lower-left corner of element */
162: xx = h * (i % m);
163: yy = h * (i / m);
164: /* node numbers for the four corners of element */
165: idx[0] = (m + 1) * (i / m) + (i % m);
166: idx[1] = idx[0] + 1;
167: idx[2] = idx[1] + m + 1;
168: idx[3] = idx[2] - 1;
169: PetscCall(FormElementRhs(xx, yy, h * h, v));
170: PetscCall(VecSetValues(b, 4, idx, v, ADD_VALUES));
171: }
172: PetscCall(VecAssemblyBegin(b));
173: PetscCall(VecAssemblyEnd(b));
175: /*
176: Modify matrices and right-hand side for Dirichlet boundary conditions
177: */
178: PetscCall(PetscMalloc1(4 * m, &rows));
179: for (PetscInt i = 0; i < m + 1; i++) {
180: rows[i] = i; /* bottom */
181: rows[3 * m - 1 + i] = m * (m + 1) + i; /* top */
182: }
183: count = m + 1; /* left side */
184: for (PetscInt i = m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
185: count = 2 * m; /* right side */
186: for (PetscInt i = 2 * m + 1; i < m * (m + 1); i += m + 1) rows[count++] = i;
187: for (PetscInt i = 0; i < 4 * m; i++) {
188: yy = h * (rows[i] / (m + 1));
189: PetscCall(VecSetValues(b, 1, &rows[i], &yy, INSERT_VALUES));
190: }
191: PetscCall(MatZeroRows(K, 4 * m, rows, 1.0, NULL, NULL));
192: if (mass) PetscCall(MatZeroRows(M, 4 * m, rows, 2.0, NULL, NULL));
193: PetscCall(PetscFree(rows));
195: PetscCall(VecAssemblyBegin(b));
196: PetscCall(VecAssemblyEnd(b));
197: PetscCall(VecNorm(b, NORM_2, &normb));
199: /*
200: Generate shifts
201: */
202: PetscCall(PetscCalloc2(nshift, &sigma, nshift, &sigma_imaginary));
203: for (PetscInt i = 0; i < nshift; i++) {
204: sigma[i] = lend + (rend - lend) / (nshift + 1) * (i + 1);
205: if (cmplx) {
206: #if PetscDefined(USE_COMPLEX)
207: sigma[i] += 1e-4 * (i + 1) * PETSC_i;
208: #else
209: if (i == 0 && nshift > 1) {
210: sigma[i + 1] = sigma[i];
211: sigma_imaginary[i] = 1e-4;
212: sigma_imaginary[i + 1] = -1e-4;
213: i++;
214: }
215: #endif
216: }
217: }
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Solve the shifted linear systems
221: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
223: PetscCall(MatCreateNestFromMultipleShifts(K, nshift, sigma, sigma_imaginary, M, explicitmat, SUBSET_NONZERO_PATTERN, &A));
224: PetscCall(MatCreateVecNestFromMultipleShifts(A, b, &b_nest));
225: PetscCall(MatCreateVecNestFromMultipleShifts(A, NULL, &x_nest));
227: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
228: PetscCall(KSPSetOperators(ksp, A, A));
229: PetscCall(KSPSetFromOptions(ksp));
230: PetscCall(KSPGetTolerances(ksp, &rtol, NULL, NULL, NULL));
231: PetscCall(KSPSolve(ksp, b_nest, x_nest));
233: /*
234: Check residual norm for each shifted linear system
235: */
236: for (PetscInt i = 0; i < nshift; i++) {
237: PetscCall(VecNestGetSubVec(x_nest, i, &x));
238: #if PetscDefined(USE_COMPLEX)
239: PetscCall(ComputeResidualNorm(K, sigma[i], M, x, b, r, w, 0.0, NULL, &norm));
240: PetscCall(PrintError(sigma[i], 0.0, norm, normb, PETSC_FALSE, rtol));
241: #else
242: if (sigma_imaginary[i] == 0.0) {
243: PetscCall(ComputeResidualNorm(K, sigma[i], M, x, b, r, w, 0.0, NULL, &norm));
244: PetscCall(PrintError(sigma[i], 0.0, norm, normb, PETSC_FALSE, rtol));
245: } else {
246: PetscCall(VecNestGetSubVec(x_nest, i + 1, &xi));
247: PetscCall(ComputeResidualNorm(K, sigma[i], M, x, b, r, w, -sigma_imaginary[i], xi, &norm));
248: PetscCall(ComputeResidualNorm(K, sigma[i], M, xi, NULL, r, w, sigma_imaginary[i], x, &normi));
249: norm = PetscHypotReal(norm, normi);
250: PetscCall(PrintError(sigma[i], sigma_imaginary[i], norm, normb, PETSC_TRUE, rtol));
251: i++;
252: }
253: #endif
254: }
256: /*
257: Clean up
258: */
259: PetscCall(KSPDestroy(&ksp));
260: PetscCall(MatDestroy(&K));
261: PetscCall(MatDestroy(&M));
262: PetscCall(MatDestroy(&A));
263: PetscCall(VecDestroy(&b_nest));
264: PetscCall(VecDestroy(&x_nest));
265: PetscCall(VecDestroy(&w));
266: PetscCall(VecDestroy(&r));
267: PetscCall(VecDestroy(&b));
268: PetscCall(PetscFree2(sigma, sigma_imaginary));
269: PetscCall(PetscFinalize());
270: return 0;
271: }
273: /*TEST
275: testset:
276: nsize: 2
277: args: -ksp_type {{gmres eksm}} -explicitmat {{0 1}}
278: output_file: output/ex90_1.out
279: test:
280: suffix: 1
281: test:
282: suffix: 2
283: args: -mass
285: test:
286: args: -ksp_type {{gmres eksm}} -explicitmat {{0 1}} -cmplx
287: suffix: 3
288: requires: !complex
290: TEST*/