Question on InterpVec for systems of PDEs

[thartland]

I have been working on including rigid body modes as near nullspace vectors to an AMG preconditioner for improved performance in solving a discretized system of PDEs. I have been seeing some seg faults, which I believe are internal to hypre.

I am including a simpler code for a scalar Poisson problem which tries to attach a single near nullspace vector (a constant vector of all 1's).

If I run mpirun -np 4 ./InterpVecTestPoissonAMG -solver 1 -interp_vec_variant 0, for which HYPRE_BoomerAMGSetInterpVecVariant(precond, 0); then the system is solved with no noticeable change in iteration count, etc. However, if I run with mpirun -np 4 ./InterpVecTestPoissonAMG -solver 1 -interp_vec_variant 2 then there is a segmentation fault.

** Process received signal ***
Signal: Segmentation fault (11)
Signal code: Address not mapped (1)
Failing at address: (nil)

The attached code is a variant of ex11.c with some vector interpolation additions taken from new_ij.c. I am trying to use interp_vec_variant 2 as that is a default option used in MFEM for linear systems from discretized elasticity problems. Any insight into these interpolation options and the proper usage in hypre would be greatly appreciated.

/*
   Variant of Example 11 which attempts to set 
   near null-space vectors but shows a seg-fault
   when interp_vec_variant is set to anything but 0 

   Interface:    Linear-Algebraic (IJ)

   Sample run:   mpirun -np 4 InterpVecTestPoissonAMG -solver 1 -interp_vec_variant 0

   Description:  This example solves the 2-D Laplacian eigenvalue
                 problem with zero boundary conditions on an nxn grid.
                 The number of unknowns is N=n^2. The standard 5-point
                 stencil is used, and we solve for the interior nodes
                 only.
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "HYPRE.h"
#include "HYPRE_parcsr_ls.h"
#include "HYPRE_krylov.h"

#define my_min(a,b)  (((a)<(b)) ? (a) : (b))

int main (int argc, char *argv[])
{
   int i;
   int myid, num_procs;
   int N, n;

   int ilower, iupper;
   int local_size, extra;

   int solver_id = 0;
   /* solver_id:
    * 0 --> CG with no preconditioner
    * 1 --> CG with AMG preconditioner
    */

   HYPRE_IJMatrix A;
   HYPRE_ParCSRMatrix parcsr_A;
   HYPRE_IJVector b;
   HYPRE_ParVector par_b;
   HYPRE_IJVector x;
   HYPRE_ParVector par_x;

   HYPRE_Solver solver, precond;

   HYPRE_Int interp_vec_variant; // strategy for interpolating the near-nullspace vectors

   /* set near-nullspace vectors */
   HYPRE_Int num_interp_vecs = 1;
   HYPRE_IJVector  ij_interp_vecs[num_interp_vecs];
   HYPRE_ParVector interp_vecs[num_interp_vecs];

   /* Initialize MPI */
   MPI_Init(&argc, &argv);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

   /* Initialize HYPRE */
   HYPRE_Init();

   /* Default problem parameters */
   n = 33;

   /* Parse command line */
   {
      int arg_index = 0;
      int print_usage = 0;

      while (arg_index < argc)
      {
         if ( strcmp(argv[arg_index], "-n") == 0 )
         {
            arg_index++;
            n = atoi(argv[arg_index++]);
         }
         else if ( strcmp(argv[arg_index], "-solver") == 0 )
         {
            arg_index++;
            solver_id = atoi(argv[arg_index++]);
         }
         else if(strcmp(argv[arg_index], "-interp_vec_variant") == 0)
         {
           arg_index++;
           interp_vec_variant = atoi(argv[arg_index++]);
         }
         else if ( strcmp(argv[arg_index], "-help") == 0 )
         {
            print_usage = 1;
            break;
         }
         else
         {
            arg_index++;
         }
      }

      if ((print_usage) && (myid == 0))
      {
         printf("\n");
         printf("Usage: %s [<options>]\n", argv[0]);
         printf("\n");
         printf("  -n <n>              : problem size in each direction (default: 33)\n");
         printf("  -solver <n>         : linear system solution strategy (default: 0)\n");
         printf("\n");
      }

      if (print_usage)
      {
         MPI_Finalize();
         return (0);
      }
   }

   /* Preliminaries: want at least one processor per row */
   if (n*n < num_procs) n = sqrt(num_procs) + 1;
   N = n*n; /* global number of rows */

   /* Each processor knows only of its own rows - the range is denoted by ilower
      and iupper.  Here we partition the rows. We account for the fact that
      N may not divide evenly by the number of processors. */
   local_size = N/num_procs;
   extra = N - local_size*num_procs;

   ilower = local_size*myid;
   ilower += my_min(myid, extra);

   iupper = local_size*(myid+1);
   iupper += my_min(myid+1, extra);
   iupper = iupper - 1;

   /* How many rows do I have? */
   local_size = iupper - ilower + 1;

   /* Create the matrix.
      Note that this is a square matrix, so we indicate the row partition
      size twice (since number of rows = number of cols) */
   HYPRE_IJMatrixCreate(MPI_COMM_WORLD, ilower, iupper, ilower, iupper, &A);

   /* Choose a parallel csr format storage (see the User's Manual) */
   HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);

   /* Initialize before setting coefficients */
   HYPRE_IJMatrixInitialize(A);

   /* Now go through my local rows and set the matrix entries.
      Each row has at most 5 entries. For example, if n=3:

      A = [M -I 0; -I M -I; 0 -I M]
      M = [4 -1 0; -1 4 -1; 0 -1 4]

      Note that here we are setting one row at a time, though
      one could set all the rows together (see the User's Manual).
   */
   {
      int nnz;
      double values[5];
      int cols[5];

      for (i = ilower; i <= iupper; i++)
      {
         nnz = 0;

         /* The left identity block:position i-n */
         if ((i-n)>=0)
         {
            cols[nnz] = i-n;
            values[nnz] = -1.0;
            nnz++;
         }

         /* The left -1: position i-1 */
         if (i%n)
         {
            cols[nnz] = i-1;
            values[nnz] = -1.0;
            nnz++;
         }

         /* Set the diagonal: position i */
         cols[nnz] = i;
         values[nnz] = 4.0;
         nnz++;

         /* The right -1: position i+1 */
         if ((i+1)%n)
         {
            cols[nnz] = i+1;
            values[nnz] = -1.0;
            nnz++;
         }

         /* The right identity block:position i+n */
         if ((i+n)< N)
         {
            cols[nnz] = i+n;
            values[nnz] = -1.0;
            nnz++;
         }

         /* Set the values for row i */
         HYPRE_IJMatrixSetValues(A, 1, &nnz, &i, cols, values);
      }
   }

   /* Assemble after setting the coefficients */
   HYPRE_IJMatrixAssemble(A);
   /* Get the parcsr matrix object to use */
   HYPRE_IJMatrixGetObject(A, (void**) &parcsr_A);

   /* Create sample rhs and solution vectors */
   HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&b);
   HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
   HYPRE_IJVectorInitialize(b);
   for (i = 0; i < num_interp_vecs; i++)
   {
     HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper, &ij_interp_vecs[i]);
     HYPRE_IJVectorSetObjectType(ij_interp_vecs[i], HYPRE_PARCSR);
     HYPRE_IJVectorInitialize(ij_interp_vecs[i]);
   }
   double rand_val;
   double *rhs_values;
   int    *rows;
   double *interp_vec_values;
   {
     rhs_values = (double *) calloc(local_size, sizeof(double));
     interp_vec_values = (double *) calloc(local_size, sizeof(double));
     rows = (int *) calloc(local_size, sizeof(int));
     for (i = 0; i < local_size; i++)
     {
       rand_val = (double)rand()/RAND_MAX * 2.0 -1.0;
       rhs_values[i] = rand_val;
       rows[i] = i + ilower;
       interp_vec_values[i] = 1.0;
     }
   }

   /* ------ b rhs of linear system A x =b ------- */
   HYPRE_IJVectorSetValues(b, local_size, rows, rhs_values);
   free(rhs_values);
   HYPRE_IJVectorAssemble(b);
   HYPRE_IJVectorGetObject(b, (void **) &par_b);

   for (i = 0; i < num_interp_vecs; i++)
   {
     HYPRE_IJVectorSetValues(ij_interp_vecs[i], local_size, rows, interp_vec_values);
     HYPRE_IJVectorAssemble(ij_interp_vecs[i]);
     HYPRE_IJVectorGetObject(ij_interp_vecs[i], (void **) &(interp_vecs[i]));
   }
   free(interp_vec_values);

   /* ------- x --------- */
   HYPRE_IJVectorCreate(MPI_COMM_WORLD, ilower, iupper,&x);
   HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR);
   HYPRE_IJVectorInitialize(x);
   double * x_values;
   {
     x_values = (double *) calloc(local_size, sizeof(double));
     for (i = 0; i < local_size; i++)
     {
       x_values[i] = 0.0;
     }
   }

   HYPRE_IJVectorSetValues(x, local_size, rows, x_values);
   free(x_values);
   free(rows);

   HYPRE_IJVectorAssemble(x);
   HYPRE_IJVectorGetObject(x, (void **) &par_x);

   double tolerance = 1e-10;
   int max_iterations = 5000;

   HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);

   /* Set some parameters (See Reference Manual for more parameters) */
   HYPRE_PCGSetMaxIter(solver, max_iterations); /* max iterations */
   HYPRE_PCGSetTol(solver, tolerance); /* conv. tolerance */
   HYPRE_PCGSetTwoNorm(solver, 1); /* use the two norm as the stopping criteria */
   HYPRE_PCGSetPrintLevel(solver, 2); /* prints out the iteration info */
   HYPRE_PCGSetLogging(solver, 1); /* needed to get run info later */

   /* AMG preconditioner */
   if (solver_id == 1)
   {
     HYPRE_BoomerAMGCreate(&precond);

     // AMG coarsening options:
     int coarsen_type = 10;   // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
     int agg_levels   = 1;    // number of aggressive coarsening levels
     double theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8

     // AMG interpolation options:
     int interp_type  = 6;    // 6 = extended+i, 0 = classical
     int Pmax         = 4;    // max number of elements per row in P

     // AMG relaxation options:
     int relax_type   = 8;    // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
     int relax_sweeps = 1;    // relaxation sweeps on each level

     // Additional options:
     int print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
     int max_levels   = 25;   // max number of levels in AMG hierarchy

     HYPRE_BoomerAMGSetCoarsenType(precond, coarsen_type);
     HYPRE_BoomerAMGSetAggNumLevels(precond, agg_levels);

     HYPRE_BoomerAMGSetRelaxType(precond, relax_type);

     // l1-GS
     int relax_down = 13; int relax_up = 14;
     HYPRE_BoomerAMGSetCycleRelaxType(precond, relax_down, 1);
     HYPRE_BoomerAMGSetCycleRelaxType(precond, relax_up, 2);

     HYPRE_BoomerAMGSetNumSweeps(precond, relax_sweeps);
     HYPRE_BoomerAMGSetStrongThreshold(precond, theta);
     HYPRE_BoomerAMGSetInterpType(precond, interp_type);
     HYPRE_BoomerAMGSetPMaxElmts(precond, Pmax);

     HYPRE_BoomerAMGSetPrintLevel(precond, print_level);
     HYPRE_BoomerAMGSetMaxLevels(precond, max_levels);

     // Use as a preconditioner (one V-cycle, zero tolerance)
     HYPRE_BoomerAMGSetMaxIter(precond, 1);
     HYPRE_BoomerAMGSetTol(precond, 0.0);

     /* ----------- */
     HYPRE_Int Q_max = 4;
     HYPRE_BoomerAMGSetInterpVecVariant(precond, interp_vec_variant);
     HYPRE_BoomerAMGSetInterpVecQMax(precond, Q_max);
     HYPRE_BoomerAMGSetInterpVectors(precond, num_interp_vecs, interp_vecs); 

     /* Set the CG preconditioner */
     HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve,
                         (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);
   }
   /* Now setup and solve! */
   HYPRE_ParCSRPCGSetup(solver, parcsr_A, par_b, par_x);
   HYPRE_ParCSRPCGSolve(solver, parcsr_A, par_b, par_x);

   int num_iterations;
   double final_res_norm;
   /* Run info - needed logging turned on */
   HYPRE_PCGGetNumIterations(solver, &num_iterations);
   HYPRE_PCGGetFinalRelativeResidualNorm(solver, &final_res_norm);
   if (myid == 0)
   {
     printf("\n");
     printf("Iterations = %d\n", num_iterations);
     printf("Final Relative Residual Norm = %e\n", final_res_norm);
     printf("\n");
   }

   /* free memory */
   HYPRE_ParCSRPCGDestroy(solver);
   if (solver_id == 1)
   {
     /* destroy the preconditioner */
     HYPRE_BoomerAMGDestroy(precond);
   }

   /* Clean up */
   for (i = 0; i < num_interp_vecs; i++)
   {
     HYPRE_IJVectorDestroy(ij_interp_vecs[i]);
   }
   HYPRE_IJMatrixDestroy(A);
   HYPRE_IJVectorDestroy(b);
   HYPRE_IJVectorDestroy(x);

   /* Finalize HYPRE */
   HYPRE_Finalize();

   /* Finalize MPI*/
   MPI_Finalize();

   return(0);
}

[victorapm]

@thartland, sorry for not responding to this before. Has your issue been resolved?