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


syntax highlighted by Code2HTML, v. 0.9.1