download the original source code.
1 /*
2 Example 5
3
4 Interface: Linear-Algebraic (IJ)
5
6 Compile with: make ex5
7
8 Sample run: mpirun -np 4 ex5
9
10 Description: This example solves the 2-D Laplacian problem with zero boundary
11 conditions on an n x n grid. The number of unknowns is N=n^2.
12 The standard 5-point stencil is used, and we solve for the
13 interior nodes only.
14
15 This example solves the same problem as Example 3. Available
16 solvers are AMG, PCG, and PCG with AMG or Parasails
17 preconditioners. */
18
19 #include <math.h>
20 #include "_hypre_utilities.h"
21 #include "HYPRE_krylov.h"
22 #include "HYPRE.h"
23 #include "HYPRE_parcsr_ls.h"
24
25 #include "vis.c"
26
27 int hypre_FlexGMRESModifyPCAMGExample(void *precond_data, int iterations,
28 double rel_residual_norm);
29
30
31 int main (int argc, char *argv[])
32 {
33 int i;
34 int myid, num_procs;
35 int N, n;
36
37 int ilower, iupper;
38 int local_size, extra;
39
40 int solver_id;
41 int vis, print_system;
42
43 double h, h2;
44
45 HYPRE_IJMatrix A;
46 HYPRE_ParCSRMatrix parcsr_A;
47 HYPRE_IJVector b;
48 HYPRE_ParVector par_b;
49 HYPRE_IJVector x;
50 HYPRE_ParVector par_x;
51
52 HYPRE_Solver solver, precond;
53
54 /* Initialize MPI */
55 MPI_Init(&argc, &argv);
56 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
57 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
58
59 /* Default problem parameters */
60 n = 33;
61 solver_id = 0;
62 vis = 0;
63 print_system = 0;
64
65
66 /* Parse command line */
67 {
68 int arg_index = 0;
69 int print_usage = 0;
70
71 while (arg_index < argc)
72 {
73 if ( strcmp(argv[arg_index], "-n") == 0 )
74 {
75 arg_index++;
76 n = atoi(argv[arg_index++]);
77 }
78 else if ( strcmp(argv[arg_index], "-solver") == 0 )
79 {
80 arg_index++;
81 solver_id = atoi(argv[arg_index++]);
82 }
83 else if ( strcmp(argv[arg_index], "-vis") == 0 )
84 {
85 arg_index++;
86 vis = 1;
87 }
88 else if ( strcmp(argv[arg_index], "-print_system") == 0 )
89 {
90 arg_index++;
91 print_system = 1;
92 }
93 else if ( strcmp(argv[arg_index], "-help") == 0 )
94 {
95 print_usage = 1;
96 break;
97 }
98 else
99 {
100 arg_index++;
101 }
102 }
103
104 if ((print_usage) && (myid == 0))
105 {
106 printf("\n");
107 printf("Usage: %s [<options>]\n", argv[0]);
108 printf("\n");
109 printf(" -n <n> : problem size in each direction (default: 33)\n");
110 printf(" -solver <ID> : solver ID\n");
111 printf(" 0 - AMG (default) \n");
112 printf(" 1 - AMG-PCG\n");
113 printf(" 8 - ParaSails-PCG\n");
114 printf(" 50 - PCG\n");
115 printf(" 61 - AMG-FlexGMRES\n");
116 printf(" -vis : save the solution for GLVis visualization\n");
117 printf(" -print_system : print the matrix and rhs\n");
118 printf("\n");
119 }
120
121 if (print_usage)
122 {
123 MPI_Finalize();
124 return (0);
125 }
126 }
127
128 /* Preliminaries: want at least one processor per row */
129 if (n*n < num_procs) n = sqrt(num_procs) + 1;
130 N = n*n; /* global number of rows */
131 h = 1.0/(n+1); /* mesh size*/
132 h2 = h*h;
133
134 /* Each processor knows only of its own rows - the range is denoted by ilower
135 and upper. Here we partition the rows. We account for the fact that
136 N may not divide evenly by the number of processors. */
137 local_size = N/num_procs;
138 extra = N - local_size*num_procs;
139
140 ilower = local_size*myid;
141 ilower += hypre_min(myid, extra);
142
143 iupper = local_size*(myid+1);
144 iupper += hypre_min(myid+1, extra);
145 iupper = iupper - 1;
146
147 /* How many rows do I have? */
148 local_size = iupper - ilower + 1;
149
150 /* Create the matrix.
151 Note that this is a square matrix, so we indicate the row partition
152 size twice (since number of rows = number of cols) */
153 HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);
154
155 /* Choose a parallel csr format storage (see the User's Manual) */
156 HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);
157
158 /* Initialize before setting coefficients */
159 HYPRE_IJMatrixInitialize(A);
160
161 /* Now go through my local rows and set the matrix entries.
162 Each row has at most 5 entries. For example, if n=3:
163
164 A = [M -I 0; -I M -I; 0 -I M]
165 M = [4 -1 0; -1 4 -1; 0 -1 4]
166
167 Note that here we are setting one row at a time, though
168 one could set all the rows together (see the User's Manual).
169 */
170 {
171 int nnz;
172 double values[5];
173 int cols[5];
174
175 for (i = ilower; i <= iupper; i++)
176 {
177 nnz = 0;
178
179 /* The left identity block:position i-n */
180 if ((i-n)>=0)
181 {
182 cols[nnz] = i-n;
183 values[nnz] = -1.0;
184 nnz++;
185 }
186
187 /* The left -1: position i-1 */
188 if (i%n)
189 {
190 cols[nnz] = i-1;
191 values[nnz] = -1.0;
192 nnz++;
193 }
194
195 /* Set the diagonal: position i */
196 cols[nnz] = i;
197 values[nnz] = 4.0;
198 nnz++;
199
200 /* The right -1: position i+1 */
201 if ((i+1)%n)
202 {
203 cols[nnz] = i+1;
204 values[nnz] = -1.0;
205 nnz++;
206 }
207
208 /* The right identity block:position i+n */
209 if ((i+n)< N)
210 {
211 cols[nnz] = i+n;
212 values[nnz] = -1.0;
213 nnz++;
214 }
215
216 /* Set the values for row i */
217 HYPRE_IJMatrixSetValues(A, 1, &nnz, &i, cols, values);
218 }
219 }
220
221 /* Assemble after setting the coefficients */
222 HYPRE_IJMatrixAssemble(A);
223
224 /* Note: for the testing of small problems, one may wish to read
225 in a matrix in IJ format (for the format, see the output files
226 from the -print_system option).
227 In this case, one would use the following routine:
228 HYPRE_IJMatrixRead( <filename>, MPI_COMM_WORLD,
229 HYPRE_PARCSR, &A );
230 <filename> = IJ.A.out to read in what has been printed out
231 by -print_system (processor numbers are omitted).
232 A call to HYPRE_IJMatrixRead is an *alternative* to the
233 following sequence of HYPRE_IJMatrix calls:
234 Create, SetObjectType, Initialize, SetValues, and Assemble
235 */
236
237
238 /* Get the parcsr matrix object to use */
239 HYPRE_IJMatrixGetObject(A, (void**) &parcsr_A);
240
241
242 /* Create the rhs and solution */
243 HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&b);
244 HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
245 HYPRE_IJVectorInitialize(b);
246
247 HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&x);
248 HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR);
249 HYPRE_IJVectorInitialize(x);
250
251 /* Set the rhs values to h^2 and the solution to zero */
252 {
253 double *rhs_values, *x_values;
254 int *rows;
255
256 rhs_values = (double*) calloc(local_size, sizeof(double));
257 x_values = (double*) calloc(local_size, sizeof(double));
258 rows = (int*) calloc(local_size, sizeof(int));
259
260 for (i=0; i<local_size; i++)
261 {
262 rhs_values[i] = h2;
263 x_values[i] = 0.0;
264 rows[i] = ilower + i;
265 }
266
267 HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values);
268 HYPRE_IJVectorSetValues(x, local_size, rows, x_values);
269
270 free(x_values);
271 free(rhs_values);
272 free(rows);
273 }
274
275
276 HYPRE_IJVectorAssemble(b);
277 /* As with the matrix, for testing purposes, one may wish to read in a rhs:
278 HYPRE_IJVectorRead( <filename>, MPI_COMM_WORLD,
279 HYPRE_PARCSR, &b );
280 as an alternative to the
281 following sequence of HYPRE_IJVectors calls:
282 Create, SetObjectType, Initialize, SetValues, and Assemble
283 */
284 HYPRE_IJVectorGetObject(b, (void **) &par_b);
285
286 HYPRE_IJVectorAssemble(x);
287 HYPRE_IJVectorGetObject(x, (void **) &par_x);
288
289
290 /* Print out the system - files names will be IJ.out.A.XXXXX
291 and IJ.out.b.XXXXX, where XXXXX = processor id */
292 if (print_system)
293 {
294 HYPRE_IJMatrixPrint(A, "IJ.out.A");
295 HYPRE_IJVectorPrint(b, "IJ.out.b");
296 }
297
298
299 /* Choose a solver and solve the system */
300
301 /* AMG */
302 if (solver_id == 0)
303 {
304 int num_iterations;
305 double final_res_norm;
306
307 /* Create solver */
308 HYPRE_BoomerAMGCreate(&solver);
309
310 /* Set some parameters (See Reference Manual for more parameters) */
311 HYPRE_BoomerAMGSetPrintLevel(solver, 3); /* print solve info + parameters */
312 HYPRE_BoomerAMGSetOldDefault(solver); /* Falgout coarsening with modified classical interpolaiton */
313 HYPRE_BoomerAMGSetRelaxType(solver, 3); /* G-S/Jacobi hybrid relaxation */
314 HYPRE_BoomerAMGSetRelaxOrder(solver, 1); /* uses C/F relaxation */
315 HYPRE_BoomerAMGSetNumSweeps(solver, 1); /* Sweeeps on each level */
316 HYPRE_BoomerAMGSetMaxLevels(solver, 20); /* maximum number of levels */
317 HYPRE_BoomerAMGSetTol(solver, 1e-7); /* conv. tolerance */
318
319 /* Now setup and solve! */
320 HYPRE_BoomerAMGSetup(solver, parcsr_A, par_b, par_x);
321 HYPRE_BoomerAMGSolve(solver, parcsr_A, par_b, par_x);
322
323 /* Run info - needed logging turned on */
324 HYPRE_BoomerAMGGetNumIterations(solver, &num_iterations);
325 HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &final_res_norm);
326 if (myid == 0)
327 {
328 printf("\n");
329 printf("Iterations = %d\n", num_iterations);
330 printf("Final Relative Residual Norm = %e\n", final_res_norm);
331 printf("\n");
332 }
333
334 /* Destroy solver */
335 HYPRE_BoomerAMGDestroy(solver);
336 }
337 /* PCG */
338 else if (solver_id == 50)
339 {
340 int num_iterations;
341 double final_res_norm;
342
343 /* Create solver */
344 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
345
346 /* Set some parameters (See Reference Manual for more parameters) */
347 HYPRE_PCGSetMaxIter(solver, 1000); /* max iterations */
348 HYPRE_PCGSetTol(solver, 1e-7); /* conv. tolerance */
349 HYPRE_PCGSetTwoNorm(solver, 1); /* use the two norm as the stopping criteria */
350 HYPRE_PCGSetPrintLevel(solver, 2); /* prints out the iteration info */
351 HYPRE_PCGSetLogging(solver, 1); /* needed to get run info later */
352
353 /* Now setup and solve! */
354 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
355 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
356
357 /* Run info - needed logging turned on */
358 HYPRE_PCGGetNumIterations(solver, &num_iterations);
359 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
360 if (myid == 0)
361 {
362 printf("\n");
363 printf("Iterations = %d\n", num_iterations);
364 printf("Final Relative Residual Norm = %e\n", final_res_norm);
365 printf("\n");
366 }
367
368 /* Destroy solver */
369 HYPRE_ParCSRPCGDestroy(solver);
370 }
371 /* PCG with AMG preconditioner */
372 else if (solver_id == 1)
373 {
374 int num_iterations;
375 double final_res_norm;
376
377 /* Create solver */
378 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
379
380 /* Set some parameters (See Reference Manual for more parameters) */
381 HYPRE_PCGSetMaxIter(solver, 1000); /* max iterations */
382 HYPRE_PCGSetTol(solver, 1e-7); /* conv. tolerance */
383 HYPRE_PCGSetTwoNorm(solver, 1); /* use the two norm as the stopping criteria */
384 HYPRE_PCGSetPrintLevel(solver, 2); /* print solve info */
385 HYPRE_PCGSetLogging(solver, 1); /* needed to get run info later */
386
387 /* Now set up the AMG preconditioner and specify any parameters */
388 HYPRE_BoomerAMGCreate(&precond);
389 HYPRE_BoomerAMGSetPrintLevel(precond, 1); /* print amg solution info */
390 HYPRE_BoomerAMGSetCoarsenType(precond, 6);
391 HYPRE_BoomerAMGSetOldDefault(precond);
392 HYPRE_BoomerAMGSetRelaxType(precond, 6); /* Sym G.S./Jacobi hybrid */
393 HYPRE_BoomerAMGSetNumSweeps(precond, 1);
394 HYPRE_BoomerAMGSetTol(precond, 0.0); /* conv. tolerance zero */
395 HYPRE_BoomerAMGSetMaxIter(precond, 1); /* do only one iteration! */
396
397 /* Set the PCG preconditioner */
398 HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
399 (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);
400
401 /* Now setup and solve! */
402 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
403 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
404
405 /* Run info - needed logging turned on */
406 HYPRE_PCGGetNumIterations(solver, &num_iterations);
407 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
408 if (myid == 0)
409 {
410 printf("\n");
411 printf("Iterations = %d\n", num_iterations);
412 printf("Final Relative Residual Norm = %e\n", final_res_norm);
413 printf("\n");
414 }
415
416 /* Destroy solver and preconditioner */
417 HYPRE_ParCSRPCGDestroy(solver);
418 HYPRE_BoomerAMGDestroy(precond);
419 }
420 /* PCG with Parasails Preconditioner */
421 else if (solver_id == 8)
422 {
423 int num_iterations;
424 double final_res_norm;
425
426 int sai_max_levels = 1;
427 double sai_threshold = 0.1;
428 double sai_filter = 0.05;
429 int sai_sym = 1;
430
431 /* Create solver */
432 HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
433
434 /* Set some parameters (See Reference Manual for more parameters) */
435 HYPRE_PCGSetMaxIter(solver, 1000); /* max iterations */
436 HYPRE_PCGSetTol(solver, 1e-7); /* conv. tolerance */
437 HYPRE_PCGSetTwoNorm(solver, 1); /* use the two norm as the stopping criteria */
438 HYPRE_PCGSetPrintLevel(solver, 2); /* print solve info */
439 HYPRE_PCGSetLogging(solver, 1); /* needed to get run info later */
440
441 /* Now set up the ParaSails preconditioner and specify any parameters */
442 HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &precond);
443
444 /* Set some parameters (See Reference Manual for more parameters) */
445 HYPRE_ParaSailsSetParams(precond, sai_threshold, sai_max_levels);
446 HYPRE_ParaSailsSetFilter(precond, sai_filter);
447 HYPRE_ParaSailsSetSym(precond, sai_sym);
448 HYPRE_ParaSailsSetLogging(precond, 3);
449
450 /* Set the PCG preconditioner */
451 HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSolve,
452 (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSetup, precond);
453
454 /* Now setup and solve! */
455 HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
456 HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);
457
458
459 /* Run info - needed logging turned on */
460 HYPRE_PCGGetNumIterations(solver, &num_iterations);
461 HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
462 if (myid == 0)
463 {
464 printf("\n");
465 printf("Iterations = %d\n", num_iterations);
466 printf("Final Relative Residual Norm = %e\n", final_res_norm);
467 printf("\n");
468 }
469
470 /* Destory solver and preconditioner */
471 HYPRE_ParCSRPCGDestroy(solver);
472 HYPRE_ParaSailsDestroy(precond);
473 }
474 /* Flexible GMRES with AMG Preconditioner */
475 else if (solver_id == 61)
476 {
477 int num_iterations;
478 double final_res_norm;
479 int restart = 30;
480 int modify = 1;
481
482
483 /* Create solver */
484 HYPRE_ParCSRFlexGMRESCreate(MPI_COMM_WORLD, &solver);
485
486 /* Set some parameters (See Reference Manual for more parameters) */
487 HYPRE_FlexGMRESSetKDim(solver, restart);
488 HYPRE_FlexGMRESSetMaxIter(solver, 1000); /* max iterations */
489 HYPRE_FlexGMRESSetTol(solver, 1e-7); /* conv. tolerance */
490 HYPRE_FlexGMRESSetPrintLevel(solver, 2); /* print solve info */
491 HYPRE_FlexGMRESSetLogging(solver, 1); /* needed to get run info later */
492
493
494 /* Now set up the AMG preconditioner and specify any parameters */
495 HYPRE_BoomerAMGCreate(&precond);
496 HYPRE_BoomerAMGSetPrintLevel(precond, 1); /* print amg solution info */
497 HYPRE_BoomerAMGSetCoarsenType(precond, 6);
498 HYPRE_BoomerAMGSetOldDefault(precond);
499 HYPRE_BoomerAMGSetRelaxType(precond, 6); /* Sym G.S./Jacobi hybrid */
500 HYPRE_BoomerAMGSetNumSweeps(precond, 1);
501 HYPRE_BoomerAMGSetTol(precond, 0.0); /* conv. tolerance zero */
502 HYPRE_BoomerAMGSetMaxIter(precond, 1); /* do only one iteration! */
503
504 /* Set the FlexGMRES preconditioner */
505 HYPRE_FlexGMRESSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
506 (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);
507
508
509 if (modify)
510 /* this is an optional call - if you don't call it, hypre_FlexGMRESModifyPCDefault
511 is used - which does nothing. Otherwise, you can define your own, similar to
512 the one used here */
513 HYPRE_FlexGMRESSetModifyPC( solver,
514 (HYPRE_PtrToModifyPCFcn) hypre_FlexGMRESModifyPCAMGExample);
515
516
517 /* Now setup and solve! */
518 HYPRE_ParCSRFlexGMRESSetup(solver, parcsr_A, par_b, par_x);
519 HYPRE_ParCSRFlexGMRESSolve(solver, parcsr_A, par_b, par_x);
520
521 /* Run info - needed logging turned on */
522 HYPRE_FlexGMRESGetNumIterations(solver, &num_iterations);
523 HYPRE_FlexGMRESGetFinalRelativeResidualNorm(solver, &final_res_norm);
524 if (myid == 0)
525 {
526 printf("\n");
527 printf("Iterations = %d\n", num_iterations);
528 printf("Final Relative Residual Norm = %e\n", final_res_norm);
529 printf("\n");
530 }
531
532 /* Destory solver and preconditioner */
533 HYPRE_ParCSRFlexGMRESDestroy(solver);
534 HYPRE_BoomerAMGDestroy(precond);
535
536 }
537 else
538 {
539 if (myid ==0) printf("Invalid solver id specified.\n");
540 }
541
542 /* Save the solution for GLVis visualization, see vis/glvis-ex5.sh */
543 if (vis)
544 {
545 FILE *file;
546 char filename[255];
547
548 int nvalues = local_size;
549 int *rows = (int*) calloc(nvalues, sizeof(int));
550 double *values = (double*) calloc(nvalues, sizeof(double));
551
552 for (i = 0; i < nvalues; i++)
553 rows[i] = ilower + i;
554
555 /* get the local solution */
556 HYPRE_IJVectorGetValues(x, nvalues, rows, values);
557
558 sprintf(filename, "%s.%06d", "vis/ex5.sol", myid);
559 if ((file = fopen(filename, "w")) == NULL)
560 {
561 printf("Error: can't open output file %s\n", filename);
562 MPI_Finalize();
563 exit(1);
564 }
565
566 /* save solution */
567 for (i = 0; i < nvalues; i++)
568 fprintf(file, "%.14e\n", values[i]);
569
570 fflush(file);
571 fclose(file);
572
573 free(rows);
574 free(values);
575
576 /* save global finite element mesh */
577 if (myid == 0)
578 GLVis_PrintGlobalSquareMesh("vis/ex5.mesh", n-1);
579 }
580
581 /* Clean up */
582 HYPRE_IJMatrixDestroy(A);
583 HYPRE_IJVectorDestroy(b);
584 HYPRE_IJVectorDestroy(x);
585
586 /* Finalize MPI*/
587 MPI_Finalize();
588
589 return(0);
590 }
591
592 /*--------------------------------------------------------------------------
593 hypre_FlexGMRESModifyPCAMGExample -
594
595 This is an example (not recommended)
596 of how we can modify things about AMG that
597 affect the solve phase based on how FlexGMRES is doing...For
598 another preconditioner it may make sense to modify the tolerance..
599
600 *--------------------------------------------------------------------------*/
601
602 int hypre_FlexGMRESModifyPCAMGExample(void *precond_data, int iterations,
603 double rel_residual_norm)
604 {
605
606
607 if (rel_residual_norm > .1)
608 {
609 HYPRE_BoomerAMGSetNumSweeps((HYPRE_Solver)precond_data, 10);
610 }
611 else
612 {
613 HYPRE_BoomerAMGSetNumSweeps((HYPRE_Solver)precond_data, 1);
614 }
615
616
617 return 0;
618 }
syntax highlighted by Code2HTML, v. 0.9.1