adding the missing makefile and updating some of the files, cleanup, etc

miniapps/CMakeLists.txt

@@ -37,3 +37,4 @@ add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
 add_subdirectory(hdiv-linear-solver)
+add_subdirectory(interiorpointsolver)

miniapps/interiorpointsolver/IPsolver.cpp

@@ -0,0 +1,964 @@
+#include "mfem.hpp"
+#include "IPsolver.hpp"
+#include "Problem.hpp"
+#include <fstream>
+#include <iostream>
+#include <cstdlib>
+
+using namespace std;
+using namespace mfem;
+
+
+ParInteriorPointSolver::ParInteriorPointSolver(ParGeneralOptProblem * problem_) 
+                     : problem(problem_), 
+                       block_offsetsumlz(5), block_offsetsuml(4), block_offsetsx(3),
+                       Huu(nullptr), Hum(nullptr), Hmu(nullptr), 
+                       Hmm(nullptr), Wmm(nullptr), D(nullptr), 
+                       Ju(nullptr), Jm(nullptr), JuT(nullptr), JmT(nullptr), 
+                       saveIterates(false)
+{
+   OptTol  = 1.e-2;
+   max_iter = 20;
+   mu_k     = 1.0;
+
+   sMax     = 1.e2;
+   kSig     = 1.e10;   // control deviation from primal Hessian
+   tauMin   = 0.99;     // control rate at which iterates can approach the boundary
+   eta      = 1.e-4;   // backtracking constant
+   thetaMin = 1.e-4;   // allowed violation of the equality constraints
+
+   // constants in line-step A-5.4
+   delta    = 1.0;
+   sTheta   = 1.1;
+   sPhi     = 2.3;
+
+   // control the rate at which the penalty parameter is decreased
+   kMu     = 0.2;
+   thetaMu = 1.5;
+
+   thetaMax = 1.e6; // maximum constraint violation
+   // data for the second order correction
+   kSoc     = 0.99;
+
+   // equation (18)
+   gTheta = 1.e-5;
+   gPhi   = 1.e-5;
+
+   kEps   = 1.e1;
+
+   dimU = problem->GetDimU();
+   dimM = problem->GetDimM();
+   dimC = problem->GetDimC();
+   ckSoc.SetSize(dimC);
+  
+   block_offsetsumlz[0] = 0;
+   block_offsetsumlz[1] = dimU; // u
+   block_offsetsumlz[2] = dimM; // m
+   block_offsetsumlz[3] = dimC; // lambda
+   block_offsetsumlz[4] = dimM; // zl
+   block_offsetsumlz.PartialSum();
+   
+   MPI_Allreduce(&dimU, &dimUGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+   MPI_Allreduce(&dimM, &dimMGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+   MPI_Allreduce(&dimC, &dimCGlb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+
+   for(int i = 0; i < block_offsetsuml.Size(); i++)  
+   { 
+      block_offsetsuml[i] = block_offsetsumlz[i]; 
+   }
+   for(int i = 0; i < block_offsetsx.Size(); i++)    
+   { 
+      block_offsetsx[i] = block_offsetsuml[i] ; 
+   }
+
+  
+   ml = problem->Getml();
+  
+   lk.SetSize(dimC);  lk  = 0.0;
+   zlk.SetSize(dimM); zlk = 0.0;
+
+   linSolver = 0;
+   linSolveTol = 1.e-8;
+   MyRank = Mpi::WorldRank();
+   iAmRoot = MyRank == 0 ? true : false;
+}
+
+double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xl, Vector &xhat, double tau)
+{
+   double alphaMaxloc = 1.0;
+   double alphaTmp;
+   for(int i = 0; i < x.Size(); i++)
+   {   
+      if( xhat(i) < 0. )
+      {
+         alphaTmp = -1. * tau * (x(i) - xl(i)) / xhat(i);
+         alphaMaxloc = min(alphaMaxloc, alphaTmp);
+      } 
+   }
+
+   // alphaMaxloc is the local maximum step size which is
+   // distinct on each MPI process. Need to compute
+   // the global maximum step size 
+   double alphaMaxglb;
+   MPI_Allreduce(&alphaMaxloc, &alphaMaxglb, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+   return alphaMaxglb;
+}
+
+double ParInteriorPointSolver::MaxStepSize(Vector &x, Vector &xhat, double tau)
+{
+   Vector zero(x.Size()); zero = 0.0;
+   return MaxStepSize(x, zero, xhat, tau);
+}
+
+
+void ParInteriorPointSolver::Mult(const Vector &x0, Vector &xf)
+{
+   BlockVector x0block(block_offsetsx); x0block = 0.0;
+   x0block.GetBlock(0).Set(1.0, x0);
+   
+   
+   ParOptProblem * OptProblem = dynamic_cast<ParOptProblem *>(problem);
+   if (dimM > 0)
+   {
+      if (true)//OptProblem == nullptr)
+      { 
+	 // fixed initialization     
+         x0block.GetBlock(1) = 1.0;
+	 x0block.GetBlock(1).Add(1.0, ml);
+      }
+      else
+      {
+         // use g(d) - s = 0
+         x0block.GetBlock(1) = 0.0;
+         x0block.GetBlock(1).Add(1.0, ml);
+         Vector c0(dimC); c0 = 0.0;	 
+	 problem->c(x0block, c0);
+	 Vector dm(dimM); dm = 0.0;
+	 for (int i = 0; i < dimM; i++)
+	 {
+	    dm(i) = max(1.e0 - x0block(dimU + i), c0(i));
+	 }
+	 x0block.GetBlock(1).Add(1.0, dm);
+      }
+   }  
+   
+   BlockVector xfblock(block_offsetsx); xfblock = 0.0;
+   Mult(x0block, xfblock);
+   xf.Set(1.0, xfblock.GetBlock(0));
+}
+
+
+void ParInteriorPointSolver::Mult(const BlockVector &x0, BlockVector &xf)
+{
+   converged = false;
+   
+   BlockVector xk(block_offsetsx), xhat(block_offsetsx); xk = 0; xhat = 0.0;
+   BlockVector Xk(block_offsetsumlz), Xhat(block_offsetsumlz); Xk = 0.0; Xhat = 0.0;
+   BlockVector Xhatuml(block_offsetsuml); Xhatuml = 0.0;
+   Vector zlhat(dimM); zlhat = 0.0;
+
+   xk.GetBlock(0).Set(1.0, x0.GetBlock(0));
+   xk.GetBlock(1).Set(1.0, x0.GetBlock(1));
+   
+   // running estimate of the final values of the Lagrange multipliers
+   lk = 0.0;
+   zlk = 0.0;
+
+   for(int i = 0; i < dimM; i++)
+   {
+      zlk(i) = 1.e1 * mu_k / (xk(i+dimU) - ml(i));
+   }
+
+   Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
+   Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
+   Xk.GetBlock(2).Set(1.0, lk);
+   Xk.GetBlock(3).Set(1.0, zlk);
+
+  /* set theta0 = theta(x0)
+   *     thetaMin
+   *     thetaMax
+   * when theta(xk) < thetaMin and the switching condition holds
+   * then we ask for the Armijo sufficient decrease of the barrier
+   * objective to be satisfied, in order to accept the trial step length alphakl
+   * 
+   * thetaMax controls how the filter is initialized for each log-barrier subproblem
+   * F0 = {(th, phi) s.t. th > thetaMax}
+   * that is the filter does not allow for iterates where the constraint violation
+   * is larger than that of thetaMax
+   */
+   double theta0 = theta(xk);
+   thetaMin = 1.e-4 * max(1.0, theta0);
+   thetaMax = 1.e8  * thetaMin; // 1.e4 * max(1.0, theta0)
+
+   double Eeval, maxBarrierSolves, Eevalmu0;
+   bool printOptimalityError; // control optimality error print to console for log-barrier subproblems
+   
+   maxBarrierSolves = 10;
+   bool smallStep;
+   int numSmallSteps = 0;
+   for(jOpt = 0; jOpt < max_iter; jOpt++)
+   {
+      if(iAmRoot)
+      {
+         cout << "interior-point solve step " << jOpt << endl;
+      }
+      // A-2. Check convergence of overall optimization problem
+      printOptimalityError = false;
+      Eevalmu0 = E(xk, lk, zlk, printOptimalityError);
+      if(Eevalmu0 < OptTol)
+      {
+         converged = true;
+	 if(iAmRoot)
+         {
+            cout << "solved optimization problem :)\n";
+         }
+         break;
+      }
+      
+      if(jOpt > 0) { maxBarrierSolves = 1; }
+      
+      for(int i = 0; i < maxBarrierSolves; i++)
+      {
+         // A-3. Check convergence of the barrier subproblem
+         printOptimalityError = true;
+         Eeval = E(xk, lk, zlk, mu_k, printOptimalityError);
+         if(iAmRoot)
+         {
+            cout << "E = " << Eeval << endl;
+         }
+         if(Eeval < kEps * mu_k)
+         {
+            if(iAmRoot)
+            {
+               cout << "solved barrier subproblem :), for mu = " << mu_k << endl;
+            }
+            // A-3.1. Recompute the barrier parameter
+            mu_k  = max(OptTol / 10., min(kMu * mu_k, pow(mu_k, thetaMu)));
+            // A-3.2. Re-initialize the filter
+            F1.DeleteAll();
+            F2.DeleteAll();
+         }
+         else
+         {
+            break;
+         }
+      }
+    
+      // A-4. Compute the search direction
+      // solve for (uhat, mhat, lhat)
+      if(iAmRoot)
+      {
+         cout << "\n** A-4. IP-Newton solve **\n";
+      }
+      zlhat = 0.0; Xhatuml = 0.0;
+      IPNewtonSolve(xk, lk, zlk, zlhat, Xhatuml, smallStep, mu_k, false); 
+      if (smallStep && numSmallSteps < 5)
+      {
+         numSmallSteps += 1;
+	 xk.GetBlock(0).Add(1.0, Xhatuml.GetBlock(0));
+         xk.GetBlock(1).Add(1.0, Xhatuml.GetBlock(1));
+         lk.Add(1.0,   Xhatuml.GetBlock(2));
+         zlk.Add(1.0, zlhat);
+         continue;	 
+      }
+      numSmallSteps = 0;
+
+      // assign data stack, X = (u, m, l, zl)
+      Xk = 0.0;
+      Xk.GetBlock(0).Set(1.0, xk.GetBlock(0));
+      Xk.GetBlock(1).Set(1.0, xk.GetBlock(1));
+      Xk.GetBlock(2).Set(1.0, lk);
+      Xk.GetBlock(3).Set(1.0, zlk);
+
+      // assign data stack, Xhat = (uhat, mhat, lhat, zlhat)
+      Xhat = 0.0;
+      for(int i = 0; i < 3; i++)
+      {
+         Xhat.GetBlock(i).Set(1.0, Xhatuml.GetBlock(i));
+      }
+      Xhat.GetBlock(3).Set(1.0, zlhat);
+
+      // A-5. Backtracking line search.
+      if(iAmRoot)
+      {
+         cout << "\n** A-5. Linesearch **\n";
+         cout << "mu = " << mu_k << endl;
+      }
+      lineSearch(Xk, Xhat, mu_k);
+
+      if(lineSearchSuccess)
+      {
+         if(iAmRoot)
+         {
+            cout << "lineSearch successful :)\n";
+         }
+         if(!switchCondition || !sufficientDecrease)
+         {
+            F1.Append( (1. - gTheta) * thx0);
+            F2.Append( phx0 - gPhi * thx0);
+         }
+         // ----- A-6: Accept the trial point
+         xk.GetBlock(0).Add(alpha, Xhat.GetBlock(0));
+         xk.GetBlock(1).Add(alpha, Xhat.GetBlock(1));
+         lk.Add(alpha,   Xhat.GetBlock(2));
+         zlk.Add(alphaz, Xhat.GetBlock(3));
+         projectZ(xk, zlk, mu_k);
+      }
+      else
+      {
+         if(iAmRoot)
+         {
+            cout << "lineSearch not successful :(\n";
+            cout << "attempting feasibility restoration with theta = " << thx0 << endl;
+         }
+	 FeasibilityRestoration(xk, lk, zlk, Xk, mu_k);
+         xk.GetBlock(0).Set(1.0, Xk.GetBlock(0));
+         xk.GetBlock(1).Set(1.0, Xk.GetBlock(1));
+         lk.Set(1.0,   Xk.GetBlock(2));
+         zlk.Set(1.0,  Xk.GetBlock(3));
+      }
+      if(jOpt + 1 == max_iter && iAmRoot) 
+      {  
+         cout << "maximum optimization iterations :(\n";
+      }
+   }
+   // done with optimization routine, just reassign data to xf reference so
+   // that the application code has access to the optimal point
+   xf = 0.0;
+   xf.GetBlock(0).Set(1.0, xk.GetBlock(0));
+   xf.GetBlock(1).Set(1.0, xk.GetBlock(1));
+}
+
+void ParInteriorPointSolver::FormIPNewtonMat(BlockVector & x, Vector & l, Vector &zl, BlockOperator &Ak)
+{
+   // WARNING: Huu, Hum, Hmu, Hmm should all be Hessian terms of the Lagrangian, currently we 
+   //          them by Hessian terms of the objective function and neglect the Hessian of l^T c
+
+   Huu = problem->Duuf(x); 
+   Hum = problem->Dumf(x);
+   Hmu = problem->Dmuf(x);
+   Hmm = problem->Dmmf(x);
+
+   Vector DiagLogBar(dimM); DiagLogBar = 0.0;
+   for(int ii = 0; ii < dimM; ii++)
+   {
+      DiagLogBar(ii) = zl(ii) / (x(ii+dimU) - ml(ii));
+   }
+   
+   if (saveIterates)
+   {
+      std::ofstream diagStream;
+      char diagString[100];
+      snprintf(diagString, 100, "logBarrierHessiandata/D%d.dat", jOpt);
+      diagStream.open(diagString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         diagStream << setprecision(30) << DiagLogBar(ii) << endl;
+      }
+      diagStream.close();
+
+      std::ofstream sStream;
+      char sString[100];
+      snprintf(sString, 100, "logBarrierHessiandata/s%d.dat", jOpt);
+      sStream.open(sString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         sStream << setprecision(30) << x(ii+dimU) << endl;
+      }
+      sStream.close();
+      
+      std::ofstream lStream;
+      char lString[100];
+      snprintf(lString, 100, "logBarrierHessiandata/l%d.dat", jOpt);
+      lStream.open(lString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         lStream << setprecision(30) << l(ii) << endl;
+      }
+      lStream.close();
+      
+      std::ofstream zlStream;
+      char zlString[100];
+      snprintf(zlString, 100, "logBarrierHessiandata/zl%d.dat", jOpt);
+      zlStream.open(zlString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimM; ii++)
+      {
+         zlStream << setprecision(30) << zl(ii) << endl;
+      }
+      zlStream.close();
+      
+      std::ofstream dStream;
+      char dString[100];
+      snprintf(dString, 100, "logBarrierHessiandata/d%d.dat", jOpt);
+      dStream.open(dString, ios::out | ios::trunc);
+      for(int ii = 0; ii < dimU; ii++)
+      {
+         dStream << setprecision(30) << x(ii) << endl;
+      }
+      dStream.close();
+   } 
+
+   D = GenerateHypreParMatrixFromDiagonal(problem->GetDofOffsetsM(), DiagLogBar);
+  
+   if(Hmm != nullptr)
+   {
+      Wmm = Hmm;
+      Wmm->Add(1.0, *D);
+   }
+   else
+   {
+      Wmm = D;
+   }
+
+   Ju = problem->Duc(x); JuT = Ju->Transpose();
+   Jm = problem->Dmc(x); JmT = Jm->Transpose();
+   //         IP-Newton system matrix
+   //    Ak = [[H_(u,u)  H_(u,m)   J_u^T]
+   //          [H_(m,u)  W_(m,m)   J_m^T]
+   //          [ J_u      J_m       0  ]]
+
+   Ak.SetBlock(0, 0, Huu);                         Ak.SetBlock(0, 2, JuT);
+                           Ak.SetBlock(1, 1, Wmm); Ak.SetBlock(1, 2, JmT);
+   Ak.SetBlock(2, 0,  Ju); Ak.SetBlock(2, 1,  Jm);
+
+   if(Hum != nullptr) { Ak.SetBlock(0, 1, Hum); Ak.SetBlock(1, 0, Hmu); }
+}
+
+// perturbed KKT system solve
+// determine the search direction
+void ParInteriorPointSolver::IPNewtonSolve(BlockVector &x, Vector &l, Vector &zl, Vector &zlhat, BlockVector &Xhat, bool & smallStep, double mu, bool socSolve)
+{
+   // solve A x = b, where A is the IP-Newton matrix
+   BlockOperator A(block_offsetsuml, block_offsetsuml); BlockVector b(block_offsetsuml); b = 0.0;
+   FormIPNewtonMat(x, l, zl, A);
+
+   //       [grad_u phi + Ju^T l]
+   // b = - [grad_m phi + Jm^T l]
+   //       [          c        ]
+   BlockVector gradphi(block_offsetsx); gradphi = 0.0;
+   BlockVector JTl(block_offsetsx); JTl = 0.0;
+   Dxphi(x, mu, gradphi);
+   
+   (A.GetBlock(0,2)).Mult(l, JTl.GetBlock(0));
+   (A.GetBlock(1,2)).Mult(l, JTl.GetBlock(1));
+
+   for(int ii = 0; ii < 2; ii++)
+   {
+      b.GetBlock(ii).Set(1.0, gradphi.GetBlock(ii));
+      b.GetBlock(ii).Add(1.0, JTl.GetBlock(ii));
+   }
+   if(!socSolve) 
+   {
+      problem->c(x, b.GetBlock(2));
+   }
+   else
+   {
+      b.GetBlock(2).Set(1.0, ckSoc);
+   }
+   b *= -1.0; 
+   Xhat = 0.0;
+
+
+   // Direct solver (default)
+   if(linSolver == 0)
+   {
+      Array2D<HypreParMatrix *> ABlockMatrix(3,3);
+      for(int ii = 0; ii < 3; ii++)
+      {
+      for(int jj = 0; jj < 3; jj++)
+      {
+         if(!A.IsZeroBlock(ii, jj))
+         {
+            ABlockMatrix(ii, jj) = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(ii, jj))));
+	 }
+	 else
+	 {
+	    ABlockMatrix(ii, jj) = nullptr;
+	 }
+      }
+      }
+      
+      HypreParMatrix * Ah = HypreParMatrixFromBlocks(ABlockMatrix);   
+      
+      /* direct solve of the 3x3 IP-Newton linear system */
+      #ifdef MFEM_USE_MUMPS
+        MUMPSSolver ASolver;
+        ASolver.SetPrintLevel(0);
+        ASolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
+        ASolver.SetOperator(*Ah);
+        ASolver.Mult(b, Xhat);
+      #else 
+        #ifdef MFEM_USE_MKL_CPARDISO
+          CPardisoSolver ASolver(MPI_COMM_WORLD);
+          ASolver.SetOperator(*Ah);
+          ASolver.Mult(b, Xhat);
+	#else
+	  MFEM_VERIFY(false, "linSolver 0 will not work unless compiled with MUMPS or MKL");
+        #endif
+      #endif
+
+      delete Ah;
+   }
+   else if(linSolver == 1 || linSolver == 2)
+   {
+     
+      // assuming Jm = -I and Hum = 0, Hmu = 0, Hmm = 0
+      ParOptProblem * tempProblem = dynamic_cast<ParOptProblem *>(problem);
+      MFEM_VERIFY(tempProblem != nullptr, "linSolver option 1 and 2 are only applicable to ParOptProblem's");      
+      
+
+
+      // form A = Huu + Ju^T D Ju, Wmm = D for contact
+      HypreParMatrix * Huuloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(0, 0))));
+      HypreParMatrix * Wmmloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(1, 1))));
+      HypreParMatrix * Juloc  = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(2, 0))));
+      HypreParMatrix * JuTloc = dynamic_cast<HypreParMatrix *>(const_cast<Operator *>(&(A.GetBlock(0, 2))));
+      
+      
+      HypreParMatrix *JuTDJu   = RAP(Wmmloc, Juloc);      // Ju^T D Ju
+      HypreParMatrix *Areduced = ParAdd(Huuloc, JuTDJu);  // Huu + Ju^T D Ju
+      /* prepare the reduced rhs 
+       * breduced = bu + Ju^T (bm + Wmm bl) */
+      Vector breduced(dimU); breduced = 0.0;
+      Vector tempVec(dimM); tempVec = 0.0;
+      Wmmloc->Mult(b.GetBlock(2), tempVec);
+      tempVec.Add(1.0, b.GetBlock(1));
+      JuTloc->Mult(tempVec, breduced);
+      breduced.Add(1.0, b.GetBlock(0));
+      
+      if(linSolver == 1)
+      {
+         // setup the solver for the reduced linear system
+         #ifdef MFEM_USE_MUMPS
+	   MUMPSSolver AreducedSolver;   
+           AreducedSolver.SetPrintLevel(0);
+           AreducedSolver.SetMatrixSymType(MUMPSSolver::MatType::SYMMETRIC_INDEFINITE);
+	   AreducedSolver.SetOperator(*Areduced);
+	   AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
+         #else 
+           #ifdef MFEM_USE_MKL_CPARDISO
+	     CPardisoSolver AreducedSolver(MPI_COMM_WORLD);
+	     AreducedSolver.SetOperator(*Areduced);
+	     AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
+           #else
+	     MFEM_VERIFY(false, "linSolver 1 will not work unless compiled with MUMPS or MKL");
+           #endif
+         #endif
+      }
+      else
+      {
+         HyprePCG AreducedSolver(MPI_COMM_WORLD);
+	 AreducedSolver.SetOperator(*Areduced);
+	 HypreBoomerAMG AreducedPrec;
+         AreducedSolver.SetTol(linSolveTol);
+         AreducedSolver.SetMaxIter(500);
+         AreducedSolver.SetPreconditioner(AreducedPrec);
+         AreducedSolver.SetResidualConvergenceOptions(); // convergence criteria based on residual norm
+	 AreducedSolver.SetPrintLevel(2);
+         AreducedSolver.Mult(breduced, Xhat.GetBlock(0));
+      }
+
+      // now propagate solved uhat to obtain mhat and lhat
+      // xm = Ju xu - bl
+      Juloc->Mult(Xhat.GetBlock(0), Xhat.GetBlock(1));
+      Xhat.GetBlock(1).Add(-1.0, b.GetBlock(2));
+
+      // xl = Wmm xm - bm
+      Wmmloc->Mult(Xhat.GetBlock(1), Xhat.GetBlock(2));
+      Xhat.GetBlock(2).Add(-1.0, b.GetBlock(1));
+
+      delete JuTDJu;
+      delete Areduced;
+   }
+
+   /* backsolve to determine zlhat */
+   for(int ii = 0; ii < dimM; ii++)
+   {
+      zlhat(ii) = -1.*(zl(ii) + (zl(ii) * Xhat(ii + dimU) - mu) / (x(ii + dimU) - ml(ii)) );
+   }
+   
+   Vector smallStepCheckVec(dimU + dimM); smallStepCheckVec = 0.0;
+   for(int ii = 0; ii < dimU + dimM; ii++)
+   {
+      smallStepCheckVec(ii) = abs(Xhat(ii)) / (1. + abs(x(ii)));
+   }
+   double smallStepCheckVal = GlobalLpNorm(infinity(), smallStepCheckVec.Normlinf(), MPI_COMM_WORLD);
+   smallStep = (smallStepCheckVal < 1.e-15) ? true : false;
+   if (smallStep && iAmRoot)
+   {
+      cout << "SMALL STEP\n";
+   }
+   
+   
+   // free memory
+   delete D;
+   delete JuT;
+   delete JmT;
+   if(Hmm != nullptr)
+   {
+      delete Wmm;
+   }
+}
+
+// here Xhat, X will be BlockVectors w.r.t. the 4 partitioning X = (u, m, l, zl)
+
+void ParInteriorPointSolver::lineSearch(BlockVector& X0, BlockVector& Xhat, double mu)
+{
+   double tau  = tauMin;
+   Vector u0   = X0.GetBlock(0);
+   Vector m0   = X0.GetBlock(1);
+   Vector l0   = X0.GetBlock(2);
+   Vector z0   = X0.GetBlock(3);
+   Vector uhat = Xhat.GetBlock(0);
+   Vector mhat = Xhat.GetBlock(1);
+   Vector lhat = Xhat.GetBlock(2);
+   Vector zhat = Xhat.GetBlock(3);
+   double alphaMax  = MaxStepSize(m0, ml, mhat, tau);
+   double alphaMaxz = MaxStepSize(z0, zhat, tau);
+   alphaz = alphaMaxz;
+
+   BlockVector x0(block_offsetsx); x0 = 0.0;
+   x0.GetBlock(0).Set(1.0, u0);
+   x0.GetBlock(1).Set(1.0, m0);
+   
+   BlockVector xhat(block_offsetsx); xhat = 0.0;
+   xhat.GetBlock(0).Set(1.0, uhat);
+   xhat.GetBlock(1).Set(1.0, mhat);
+   
+   BlockVector xtrial(block_offsetsx); xtrial = 0.0;
+   BlockVector Dxphi0(block_offsetsx); Dxphi0 = 0.0;
+   int maxBacktrack = 20;
+   alpha = alphaMax;
+
+   Vector ck0(dimC); ck0 = 0.0;
+   Vector zhatsoc(dimM); zhatsoc = 0.0;
+   BlockVector Xhatumlsoc(block_offsetsuml); Xhatumlsoc = 0.0;
+   BlockVector xhatsoc(block_offsetsx); xhatsoc = 0.0;
+   Vector uhatsoc(dimU); uhatsoc = 0.0;
+   Vector mhatsoc(dimM); mhatsoc = 0.0;
+
+   Dxphi(x0, mu, Dxphi0);
+
+   Dxphi0_xhat = InnerProduct(MPI_COMM_WORLD, Dxphi0, xhat);
+   descentDirection = Dxphi0_xhat < 0. ? true : false;
+   if(descentDirection)
+   {
+      if (iAmRoot)
+      {
+         cout << "is a descent direction for the log-barrier objective\n";
+      }
+   }
+   else
+   {
+      cout << "is not a descent direction for the log-barrier objective\n";
+   }
+   thx0 = theta(x0);
+   phx0 = phi(x0, mu);
+
+   lineSearchSuccess = false;
+   for(int i = 0; i < maxBacktrack; i++)
+   {
+      if (iAmRoot)
+      {
+         cout << "\n--------- alpha = " << alpha << " ---------\n";
+      }
+      // ----- A-5.2. Compute trial point: xtrial = x0 + alpha_i xhat
+      xtrial.Set(1.0, x0);
+      xtrial.Add(alpha, xhat);
+
+      // ------ A-5.3. if not in filter region go to A.5.4 otherwise go to A-5.5.
+      thxtrial = theta(xtrial);
+      phxtrial = phi(xtrial, mu);
+      filterCheck(thxtrial, phxtrial);    
+      if (!inFilterRegion)
+      {
+         if (iAmRoot)
+         {
+            cout << "not in filter region :)\n";
+         }
+         // ------ A.5.4: Check sufficient decrease
+         if(!descentDirection)
+         {
+            switchCondition = false;
+         }
+         else
+         {
+            switchCondition = (alpha * pow(abs(Dxphi0_xhat), sPhi) > delta * pow(thx0, sTheta)) ? true : false;
+         }
+         if (iAmRoot)
+         {
+            cout << "theta(x0) = "     << thx0     << ", thetaMin = "                  << thetaMin             << endl;
+            cout << "theta(xtrial) = " << thxtrial << ", (1-gTheta) *theta(x0) = "     << (1. - gTheta) * thx0 << endl;
+            cout << "phi(xtrial) = "   << phxtrial << ", phi(x0) - gPhi *theta(x0) = " << phx0 - gPhi * thx0   << endl;
+         }      
+         // Case I      
+         if(thx0 <= thetaMin && switchCondition)
+         {
+            sufficientDecrease = (phxtrial <= phx0 + eta * alpha * Dxphi0_xhat) ? true : false;
+            if(sufficientDecrease)
+            {
+               if(iAmRoot) { cout << "Line search successful: sufficient decrease in log-barrier objective.\n"; } 
+               // accept the trial step
+               lineSearchSuccess = true;
+               break;
+            }
+         }
+         else
+         {
+            if(thxtrial <= (1. - gTheta) * thx0 || phxtrial <= phx0 - gPhi * thx0)
+            {
+               if(iAmRoot) { cout << "Line search successful: infeasibility or log-barrier objective decreased.\n"; } 
+               // accept the trial step
+               lineSearchSuccess = true;
+               break;
+            }
+         }
+      }
+      else
+      {
+         if (iAmRoot)
+         {
+            cout << "in filter region :(\n"; 
+         }
+      }
+      alpha *= 0.5;
+
+   } 
+}
+
+
+void ParInteriorPointSolver::projectZ(const Vector &x, Vector &z, double mu)
+{
+   double zi;
+   double mudivmml;
+   for(int i = 0; i < dimM; i++)
+   {
+      zi = z(i);
+      mudivmml = mu / (x(i + dimU) - ml(i));
+      z(i) = max(min(zi, kSig * mudivmml), mudivmml / kSig);
+   }
+}
+
+void ParInteriorPointSolver::filterCheck(double th, double ph)
+{
+   inFilterRegion = false;
+   if(th > thetaMax)
+   {
+      inFilterRegion = true;
+   }
+   else
+   {
+      for(int i = 0; i < F1.Size(); i++)
+      {
+         if(th >= F1[i] && ph >= F2[i])
+         {
+            inFilterRegion = true;
+            break;
+         }
+      }
+   }
+}
+
+void ParInteriorPointSolver::FeasibilityRestoration(const BlockVector & x, const Vector &l, const Vector &zl, BlockVector &X, double mu)
+{
+   ParOptProblem * OptProblem = dynamic_cast<ParOptProblem *>(problem);
+   X.GetBlock(0).Set(1.0, x.GetBlock(0));
+   X.GetBlock(1).Set(1.0, x.GetBlock(1));
+   X.GetBlock(2).Set(1.0, l);
+   X.GetBlock(3).Set(1.0, zl);
+
+   if (OptProblem != nullptr)
+   {
+      Vector g(dimC); g = 0.0;
+      OptProblem->g(x.GetBlock(0), g);
+      Ju = OptProblem->Ddg(x.GetBlock(0));
+      Ju->DropSmallEntries(1.e-16);
+      SparseMatrix JuMerged;
+      Ju->MergeDiagAndOffd(JuMerged);
+      
+      int num_loc_modified_constraints = 0;
+      int num_glb_modified_constraints = 0;
+      for (int i = 0; i < dimC; i++)
+      {
+         if(JuMerged.RowIsEmpty(i))
+	 {
+	    if (g(i) < 1.e-15)
+	    {
+	       if (iAmRoot)
+	       {
+	           cout << "WARNING: LICQ violation detected (g_i and grad(g_i) both zero)\n";
+	       }
+	       continue;
+	    }
+	    X(dimU + i) = g(i); // s_i = \gamma
+            X(dimU + dimM + i) = -1. * mu / g(i); // l_i = -mu/\gamma	 
+	    X(dimU + dimM + dimC + i) = mu / g(i);
+	    num_loc_modified_constraints += 1;
+	 }
+      }
+      MPI_Allreduce(&num_loc_modified_constraints, &num_glb_modified_constraints, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+      if (num_glb_modified_constraints == 0)
+      {
+         if(iAmRoot)
+         {
+            cout << "trying feasibility restoration with no null rows in Jacobian\n";
+            cout << "exiting\n";
+         }
+         exit(1);
+      }   
+   }
+}
+
+
+
+double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, double mu, bool printEeval)
+{
+   double E1, E2, E3;
+   double sc, sd;
+   BlockVector gradL(block_offsetsx); gradL = 0.0; // stationarity grad L = grad f + J^T l - z
+   Vector cx(dimC); cx = 0.0;     // feasibility c = c(x)
+   Vector comp(dimM); comp = 0.0; // complementarity M Z - mu 1
+
+   DxL(x, l, zl, gradL);
+   E1 = GlobalLpNorm(infinity(), gradL.Normlinf(), MPI_COMM_WORLD);
+
+   problem->c(x, cx);
+   E2 = GlobalLpNorm(infinity(), cx.Normlinf(), MPI_COMM_WORLD); 
+
+   for(int ii = 0; ii < dimM; ii++) 
+   { 
+      comp(ii) = x(dimU + ii) * zl(ii) - mu;
+   }
+   E3 = GlobalLpNorm(infinity(), comp.Normlinf(), MPI_COMM_WORLD); 
+   double ll1, zl1;
+
+   zl1 = GlobalLpNorm(1, zl.Norml1(), MPI_COMM_WORLD); 
+   ll1 = GlobalLpNorm(1, l.Norml1(), MPI_COMM_WORLD);
+   sc = max(sMax, zl1 / (double(dimMGlb)) ) / sMax;
+   sd = max(sMax, (ll1 + zl1) / (double(dimCGlb + dimMGlb))) / sMax;
+   if(iAmRoot)
+   {
+      cout << "evaluating optimality error for mu = " << mu << endl;
+      cout << "stationarity measure = "    << E1 / sd << endl;
+      cout << "feasibility measure  = "    << E2      << endl;
+      cout << "complimentarity measure = " << E3 / sc << endl;
+   }
+   return max(max(E1 / sd, E2), E3 / sc);
+}
+
+double ParInteriorPointSolver::E(const BlockVector &x, const Vector &l, const Vector &zl, bool printEeval)
+{
+  return E(x, l, zl, 0.0, printEeval);
+}
+
+double ParInteriorPointSolver::theta(const BlockVector &x)
+{
+  Vector cx(dimC); cx = 0.0;
+  problem->c(x, cx);
+  return GlobalLpNorm(2, cx.Norml2(), MPI_COMM_WORLD);
+}
+
+// log-barrier objective
+double ParInteriorPointSolver::phi(const BlockVector &x, double mu)
+{
+   double fx = problem->CalcObjective(x); 
+   double logBarrierLoc = 0.0;
+   for(int i = 0; i < dimM; i++) 
+   { 
+     logBarrierLoc += log(x(dimU+i) - ml(i));
+   }
+   double logBarrierGlb;
+   MPI_Allreduce(&logBarrierLoc, &logBarrierGlb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+   return fx - mu * logBarrierGlb;
+}
+
+// gradient of log-barrier objective with respect to x = (u, m)
+void ParInteriorPointSolver::Dxphi(const BlockVector &x, double mu, BlockVector &y)
+{
+   problem->CalcObjectiveGrad(x, y);
+   
+   for(int i = 0; i < dimM; i++) 
+   { 
+      y(dimU + i) -= mu / (x(dimU + i));
+   } 
+}
+
+// Lagrangian function evaluation
+// L(x, l, zl) = f(x) + l^T c(x) - zl^T m
+double ParInteriorPointSolver::L(const BlockVector &x, const Vector &l, const Vector &zl)
+{
+   double fx = problem->CalcObjective(x);
+   Vector cx(dimC); problem->c(x, cx);
+   return (fx + InnerProduct(MPI_COMM_WORLD, cx, l) - InnerProduct(MPI_COMM_WORLD, x.GetBlock(1), zl));
+}
+
+void ParInteriorPointSolver::DxL(const BlockVector &x, const Vector &l, const Vector &zl, BlockVector &y)
+{
+   // evaluate the gradient of the objective with respect to the primal variables x = (u, m)
+   BlockVector gradxf(block_offsetsx); gradxf = 0.0;
+   problem->CalcObjectiveGrad(x, gradxf);
+   
+   HypreParMatrix *Jacu, *Jacm;
+   Jacu = problem->Duc(x); 
+   Jacm = problem->Dmc(x);
+   Jacu->MultTranspose(l, y.GetBlock(0));
+   Jacm->MultTranspose(l, y.GetBlock(1));   
+   
+   y.Add(1.0, gradxf);
+   (y.GetBlock(1)).Add(-1.0, zl);
+}
+
+bool ParInteriorPointSolver::GetConverged() const
+{
+   return converged;
+}
+
+void ParInteriorPointSolver::SetTol(double Tol)
+{
+   OptTol = Tol;
+}
+
+void ParInteriorPointSolver::SetMaxIter(int max_it)
+{
+   max_iter = max_it;
+}
+
+void ParInteriorPointSolver::SetBarrierParameter(double mu_0)
+{
+   mu_k = mu_0;
+}
+
+void ParInteriorPointSolver::SaveIterates(bool save)
+{
+   saveIterates = save;
+}
+
+
+void ParInteriorPointSolver::SetLinearSolver(int LinSolver)
+{
+   linSolver = LinSolver;
+}
+
+void ParInteriorPointSolver::SetLinearSolveTol(double Tol)
+{
+  linSolveTol = Tol;
+}
+
+void ParInteriorPointSolver::GetLagrangeMultiplier(Vector & y)
+{
+  y.SetSize(dimM); y = 0.;
+  y.Set(1.0, zlk);
+}
+
+
+
+ParInteriorPointSolver::~ParInteriorPointSolver() 
+{
+   F1.DeleteAll();
+   F2.DeleteAll();
+   block_offsetsx.DeleteAll();
+   block_offsetsumlz.DeleteAll();
+   block_offsetsuml.DeleteAll();
+   ml.SetSize(0);
+}

miniapps/interiorpointsolver/IPsolver.hpp

@@ -0,0 +1,83 @@
+#include "mfem.hpp"
+#include "Problem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+#ifndef PARIPSOLVER 
+#define PARIPSOLVER
+
+class ParInteriorPointSolver
+{
+protected:
+    ParGeneralOptProblem* problem;
+    double OptTol;
+    int  max_iter;
+    double mu_k; // \mu_k
+    Vector lk, zlk;
+
+    double sMax, kSig, tauMin, eta, thetaMin, delta, sTheta, sPhi, kMu, thetaMu;
+    double thetaMax, kSoc, gTheta, gPhi, kEps;
+	
+    // filter
+    Array<double> F1, F2;
+	
+    // quantities computed in lineSearch
+    double alpha, alphaz;
+    double thx0, thxtrial;
+    double phx0, phxtrial;
+    bool descentDirection, switchCondition, sufficientDecrease, lineSearchSuccess, inFilterRegion;
+    double Dxphi0_xhat;
+
+    int dimU, dimM, dimC;
+    int dimUGlb, dimMGlb, dimCGlb;
+    Array<int> block_offsetsumlz, block_offsetsuml, block_offsetsx;
+    Vector ml;
+
+    Vector ckSoc;
+    HypreParMatrix * Huu, * Hum, * Hmu, * Hmm, * Wmm, *D, * Ju, * Jm, * JuT, * JmT;
+    
+    int jOpt;
+    bool converged;
+    
+    int MyRank;
+    bool iAmRoot;
+
+    bool saveLogBarrierIterates;
+    bool saveIterates;
+    int linSolver;
+    double linSolveTol;
+public:
+    ParInteriorPointSolver(ParGeneralOptProblem*);
+    double MaxStepSize(Vector& , Vector& , Vector& , double);
+    double MaxStepSize(Vector& , Vector& , double);
+    void Mult(const BlockVector& , BlockVector&);
+    void Mult(const Vector&, Vector &); 
+    void GetLagrangeMultiplier(Vector &);
+    void FormIPNewtonMat(BlockVector& , Vector& , Vector& , BlockOperator &);
+    void IPNewtonSolve(BlockVector& , Vector& , Vector& , Vector&, BlockVector& , bool &, double, bool);
+    void lineSearch(BlockVector& , BlockVector& , double);
+    void projectZ(const Vector & , Vector &, double);
+    void filterCheck(double, double);
+    double E(const BlockVector &, const Vector &, const Vector &, double, bool);
+    double E(const BlockVector &, const Vector &, const Vector &, bool);
+    bool GetConverged() const;
+    double theta(const BlockVector &);
+    double phi(const BlockVector &, double);
+    void Dxphi(const BlockVector &, double, BlockVector &);
+    double L(const BlockVector &, const Vector &, const Vector &);
+    void DxL(const BlockVector &, const Vector &, const Vector &, BlockVector &);
+    void SetTol(double);
+    void SetMaxIter(int);
+    void SetBarrierParameter(double);    
+    void SaveIterates(bool);
+    void SetLinearSolver(int);
+    void SetLinearSolveTol(double);
+    void FeasibilityRestoration(const BlockVector &, const Vector &, const Vector &, BlockVector &, double); 
+    virtual ~ParInteriorPointSolver();
+};
+
+#endif

miniapps/interiorpointsolver/Problem.cpp

@@ -0,0 +1,404 @@
+#include "mfem.hpp"
+#include "Problem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+ParGeneralOptProblem::ParGeneralOptProblem() : block_offsetsx(3) {}
+
+void ParGeneralOptProblem::Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_)
+{
+  dofOffsetsU = new HYPRE_BigInt[2];
+  dofOffsetsM = new HYPRE_BigInt[2];
+  for(int i = 0; i < 2; i++)
+  {
+    dofOffsetsU[i] = dofOffsetsU_[i];
+    dofOffsetsM[i] = dofOffsetsM_[i];
+  }
+  dimU = dofOffsetsU[1] - dofOffsetsU[0];
+  dimM = dofOffsetsM[1] - dofOffsetsM[0];
+  dimC = dimM;
+  
+  block_offsetsx[0] = 0;
+  block_offsetsx[1] = dimU;
+  block_offsetsx[2] = dimM;
+  block_offsetsx.PartialSum();
+  
+  MPI_Allreduce(&dimU, &dimUglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&dimM, &dimMglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+}
+
+
+void ParGeneralOptProblem::CalcObjectiveGrad(const BlockVector &x, BlockVector &y) const
+{
+   Duf(x, y.GetBlock(0));
+   Dmf(x, y.GetBlock(1));
+}
+
+ParGeneralOptProblem::~ParGeneralOptProblem()
+{
+   block_offsetsx.DeleteAll();
+}
+
+
+// min E(d) s.t. g(d) >= 0
+// min_(d,s) E(d) s.t. c(d,s) := g(d) - s = 0, s >= 0
+ParOptProblem::ParOptProblem() : ParGeneralOptProblem()
+{
+}
+
+void ParOptProblem::Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_)
+{
+  dofOffsetsU = new HYPRE_BigInt[2];
+  dofOffsetsM = new HYPRE_BigInt[2];
+  for(int i = 0; i < 2; i++)
+  {
+    dofOffsetsU[i] = dofOffsetsU_[i];
+    dofOffsetsM[i] = dofOffsetsM_[i];
+  }
+
+  dimU = dofOffsetsU[1] - dofOffsetsU[0];
+  dimM = dofOffsetsM[1] - dofOffsetsM[0];
+  dimC = dimM;
+  
+  block_offsetsx[0] = 0;
+  block_offsetsx[1] = dimU;
+  block_offsetsx[2] = dimM;
+  block_offsetsx.PartialSum();
+  
+  MPI_Allreduce(&dimU, &dimUglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&dimM, &dimMglb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  
+  ml.SetSize(dimM); ml = 0.0;
+  Vector negIdentDiag(dimM);
+  negIdentDiag = -1.0;
+  Ih = GenerateHypreParMatrixFromDiagonal(dofOffsetsM, negIdentDiag);
+}
+
+
+double ParOptProblem::CalcObjective(const BlockVector &x) const { return E(x.GetBlock(0)); }
+
+void ParOptProblem::Duf(const BlockVector &x, Vector &y) const { DdE(x.GetBlock(0), y); }
+
+void ParOptProblem::Dmf(const BlockVector &x, Vector &y) const { y = 0.0; }
+
+HypreParMatrix * ParOptProblem::Duuf(const BlockVector &x) 
+{ 
+   return DddE(x.GetBlock(0)); 
+}
+
+HypreParMatrix * ParOptProblem::Dumf(const BlockVector &x) { return nullptr; }
+
+HypreParMatrix * ParOptProblem::Dmuf(const BlockVector &x) { return nullptr; }
+
+HypreParMatrix * ParOptProblem::Dmmf(const BlockVector &x) { return nullptr; }
+
+void ParOptProblem::c(const BlockVector &x, Vector &y) const // c(u,m) = g(u) - m 
+{
+   g(x.GetBlock(0), y);
+   y.Add(-1.0, x.GetBlock(1));  
+}
+
+HypreParMatrix * ParOptProblem::Duc(const BlockVector &x) 
+{ 
+   return Ddg(x.GetBlock(0)); 
+}
+
+HypreParMatrix * ParOptProblem::Dmc(const BlockVector &x) 
+{ 
+   return Ih;
+} 
+
+ParOptProblem::~ParOptProblem() 
+{
+  delete[] dofOffsetsU;
+  delete[] dofOffsetsM;
+  delete Ih;
+}
+
+
+
+
+
+// Obstacle Problem, no essential boundary conditions enforced
+// Hessian of energy term is K + M (stiffness + mass)
+ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *Vh_, 
+                                       double (*fSource)(const Vector &),
+				       double (*obstacleSource)(const Vector &)) :
+                                       ParOptProblem(), Vh(Vh_), J(nullptr) 
+{
+   Init(Vh->GetTrueDofOffsets(), Vh->GetTrueDofOffsets());
+   cout << "dimU = " << dimU;
+   f.SetSize(dimU); f = 0.0;
+   psi.SetSize(dimU); psi = 0.0;
+
+
+   Kform = new ParBilinearForm(Vh);
+   Kform->AddDomainIntegrator(new MassIntegrator);
+   Kform->AddDomainIntegrator(new DiffusionIntegrator);
+   Kform->Assemble();
+   Kform->Finalize();
+   Kform->FormSystemMatrix(ess_tdof_list, K);
+
+   FunctionCoefficient fcoeff(fSource);
+   fform = new ParLinearForm(Vh);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+   fform->Assemble();
+   Vector F(dimU);
+   fform->ParallelAssemble(F);
+   f.SetSize(dimU);
+   f.Set(1.0, F);
+
+   Vector iDiag(dimU); iDiag = 1.0;
+   SparseMatrix * Jacg = new SparseMatrix(iDiag);
+   
+   J = new HypreParMatrix(MPI_COMM_WORLD, dimUglb, dofOffsetsU, Jacg);
+   HypreStealOwnership(*J, *Jacg);
+   delete Jacg;
+}
+
+// Obstacle Problem, essential boundary conditions enforced
+// Hessian of energy term is K (stiffness)
+ParObstacleProblem::ParObstacleProblem(ParFiniteElementSpace *Vh_, 
+				       double (*fSource)(const Vector &),
+				       double (*obstacleSource)(const Vector &),
+				       Array<int> tdof_list, Vector &xDC) : ParOptProblem(), 
+	                                                                    Vh(Vh_), J(nullptr)
+{
+   Init(Vh->GetTrueDofOffsets(), Vh->GetTrueDofOffsets());
+   f.SetSize(dimU); f = 0.0;
+   psi.SetSize(dimU); psi = 0.0;
+   // elastic energy functional terms	
+   ess_tdof_list = tdof_list;
+   Kform = new ParBilinearForm(Vh);
+   Kform->AddDomainIntegrator(new DiffusionIntegrator);
+   Kform->Assemble();
+   Kform->Finalize();
+   Kform->FormSystemMatrix(ess_tdof_list, K);
+
+   FunctionCoefficient fcoeff(fSource);
+   fform = new ParLinearForm(Vh);
+   fform->AddDomainIntegrator(new DomainLFIntegrator(fcoeff));
+   fform->Assemble();
+   Vector F(dimU);
+   fform->ParallelAssemble(F);
+   f.SetSize(dimU);
+   f.Set(1.0, F);
+   Kform->EliminateVDofsInRHS(ess_tdof_list, xDC, f);
+   
+   // obstacle constraints --  
+   Vector iDiag(dimU); iDiag = 1.0;
+   for(int i = 0; i < ess_tdof_list.Size(); i++)
+   {
+     iDiag(ess_tdof_list[i]) = 0.0;
+   }
+   SparseMatrix * Jacg = new SparseMatrix(iDiag);
+
+   J = new HypreParMatrix(MPI_COMM_WORLD, dimUglb, dofOffsetsU, Jacg);
+   HypreStealOwnership(*J, *Jacg);
+   delete Jacg;
+
+   FunctionCoefficient psi_fc(obstacleSource);
+   ParGridFunction psi_gf(Vh);
+   psi_gf.ProjectCoefficient(psi_fc);
+   psi.Set(1.0, (*psi_gf.GetTrueDofs()));
+   for(int i = 0; i < ess_tdof_list.Size(); i++)
+   {
+     psi(ess_tdof_list[i]) = xDC(ess_tdof_list[i]) - 1.e-8;
+   }
+}
+
+
+
+double ParObstacleProblem::E(const Vector &d) const
+{
+   Vector Kd(K.Height()); Kd = 0.0;
+   MFEM_VERIFY(d.Size() == K.Width(), "ParObstacleProblem::E - Inconsistent dimensions");
+   K.Mult(d, Kd);
+   return 0.5 * InnerProduct(MPI_COMM_WORLD, d, Kd) - InnerProduct(MPI_COMM_WORLD, f, d);
+}
+
+void ParObstacleProblem::DdE(const Vector &d, Vector &gradE) const
+{
+   gradE.SetSize(K.Height());
+   MFEM_VERIFY(d.Size() == K.Width(), "ParObstacleProblem::DdE - Inconsistent dimensions");
+   K.Mult(d, gradE);
+   MFEM_VERIFY(f.Size() == K.Height(), "ParObstacleProblem::DdE - Inconsistent dimensions");
+   gradE.Add(-1.0, f);
+}
+
+HypreParMatrix * ParObstacleProblem::DddE(const Vector &d)
+{
+   return &K; 
+}
+
+// g(d) = d >= \psi
+void ParObstacleProblem::g(const Vector &d, Vector &gd) const
+{
+   MFEM_VERIFY(d.Size() == J->Width(), "ParObstacleProblem::g - Inconsistent dimensions");
+   J->Mult(d, gd);
+   MFEM_VERIFY(gd.Size() == J->Height(), "ParObstacleProblem::g - Inconsistent dimensions");
+   gd.Add(-1.0, psi);
+}
+
+HypreParMatrix * ParObstacleProblem::Ddg(const Vector &d)
+{
+   return J;
+}
+
+ParObstacleProblem::~ParObstacleProblem()
+{
+   delete Kform;
+   delete fform;
+   delete J;
+}
+
+
+
+ReducedProblem::ReducedProblem(ParOptProblem * problem_, HYPRE_Int * constraintMask)
+{
+  problem = problem_;
+  J = nullptr;
+  P = nullptr;
+  
+  int nprocs = Mpi::WorldSize();
+  int myrank = Mpi::WorldRank();
+
+  HYPRE_BigInt * dofOffsets = problem->GetDofOffsetsU();
+
+  // given a constraint mask, lets update the constraintOffsets
+  // from the original problem
+  int nLocConstraints = 0;
+  int nProblemConstraints = problem->GetDimM();
+  for (int i = 0; i < nProblemConstraints; i++)
+  {
+    if (constraintMask[i] == 1)
+    {
+      nLocConstraints += 1;
+    }
+  }
+
+  HYPRE_BigInt * constraintOffsets_reduced;
+  constraintOffsets_reduced = offsetsFromLocalSizes(nLocConstraints);
+
+
+  for (int i = 0; i < 2; i++)
+  {
+    cout << "constraintOffsetsReduced_" << i << " = " << constraintOffsets_reduced[i] << ", (rank = " << myrank << ")\n";
+  }
+
+  HYPRE_BigInt * constraintOffsets;
+  constraintOffsets = offsetsFromLocalSizes(nProblemConstraints);
+  for (int i = 0; i < 2; i++)
+  {
+    cout << "constraintOffsets_" << i << " = " << constraintOffsets[i] << ", (rank = " << myrank << ")\n";
+  }
+  
+  
+  P = GenerateProjector(constraintOffsets, constraintOffsets_reduced, constraintMask);
+
+  Init(dofOffsets, constraintOffsets_reduced);
+  delete[] constraintOffsets_reduced;
+  delete[] constraintOffsets;
+}
+
+ReducedProblem::ReducedProblem(ParOptProblem * problem_, HypreParVector & constraintMask)
+{
+  problem = problem_;
+  J = nullptr;
+  P = nullptr;
+  
+  int nprocs = Mpi::WorldSize();
+  int myrank = Mpi::WorldRank();
+
+  HYPRE_BigInt * dofOffsets = problem->GetDofOffsetsU();
+
+  // given a constraint mask, lets update the constraintOffsets
+  // from the original problem
+  int nLocConstraints = 0;
+  int nProblemConstraints = problem->GetDimM();
+  for (int i = 0; i < nProblemConstraints; i++)
+  {
+    if (constraintMask[i] == 1)
+    {
+      nLocConstraints += 1;
+    }
+  }
+  cout << "nLocConstraints = " << nLocConstraints << ", (rank = " << myrank << ")\n";
+
+  HYPRE_BigInt * constraintOffsets_reduced;
+  constraintOffsets_reduced = offsetsFromLocalSizes(nLocConstraints);
+
+
+  for (int i = 0; i < 2; i++)
+  {
+    cout << "constraintOffsetsReduced_" << i << " = " << constraintOffsets_reduced[i] << ", (rank = " << myrank << ")\n";
+  }
+
+  HYPRE_BigInt * constraintOffsets;
+  constraintOffsets = offsetsFromLocalSizes(nProblemConstraints);
+  for (int i = 0; i < 2; i++)
+  {
+    cout << "constraintOffsets_" << i << " = " << constraintOffsets[i] << ", (rank = " << myrank << ")\n";
+  }
+  
+  
+  P = GenerateProjector(constraintOffsets, constraintOffsets_reduced, constraintMask);
+
+  Init(dofOffsets, constraintOffsets_reduced);
+  delete[] constraintOffsets_reduced;
+  delete[] constraintOffsets;
+}
+
+// energy objective E(d)
+double ReducedProblem::E(const Vector &d) const
+{
+  return problem->E(d);
+}
+
+
+// gradient of energy objective
+void ReducedProblem::DdE(const Vector &d, Vector & gradE) const
+{
+  problem->DdE(d, gradE);
+}
+
+
+HypreParMatrix * ReducedProblem::DddE(const Vector &d)
+{
+  return problem->DddE(d);
+}
+
+void ReducedProblem::g(const Vector &d, Vector &gd) const
+{
+  Vector gdfull(problem->GetDimM()); gdfull = 0.0;
+  problem->g(d, gdfull);
+  P->Mult(gdfull, gd);
+}
+
+
+HypreParMatrix * ReducedProblem::Ddg(const Vector &d)
+{
+  HypreParMatrix * Jfull = problem->Ddg(d);
+  if (J != nullptr)
+  {
+    delete J; J = nullptr;
+  }
+  J = ParMult(P, Jfull, true);
+  return J;
+}
+
+ReducedProblem::~ReducedProblem()
+{
+  delete P;
+  if (J != nullptr)
+  {
+    delete J;
+  }
+}
+
+

miniapps/interiorpointsolver/Problem.hpp

@@ -0,0 +1,158 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "utilities.hpp"
+
+using namespace std;
+using namespace mfem;
+
+#ifndef PARPROBLEM_DEFS
+#define PARPROBLEM_DEFS
+
+// abstract ParGeneralOptProblem class
+// of the form
+// min_(u,m) f(u,m) s.t. c(u,m)=0 and m>=ml
+// the primal variable (u, m) is represented as a BlockVector
+class ParGeneralOptProblem
+{
+protected:
+    int dimU, dimM, dimC;
+    int dimUglb, dimMglb;
+    HYPRE_BigInt * dofOffsetsU;
+    HYPRE_BigInt * dofOffsetsM;
+    Array<int> block_offsetsx;
+    Vector ml;
+public:
+    ParGeneralOptProblem();
+    virtual void Init(HYPRE_BigInt * dofOffsetsU_, HYPRE_BigInt * dofOffsetsM_);
+    virtual double CalcObjective(const BlockVector &) const = 0;
+    virtual void Duf(const BlockVector &, Vector &) const = 0;
+    virtual void Dmf(const BlockVector &, Vector &) const = 0;
+    void CalcObjectiveGrad(const BlockVector &, BlockVector &) const;
+    virtual HypreParMatrix * Duuf(const BlockVector &) = 0;
+    virtual HypreParMatrix * Dumf(const BlockVector &) = 0;
+    virtual HypreParMatrix * Dmuf(const BlockVector &) = 0;
+    virtual HypreParMatrix * Dmmf(const BlockVector &) = 0;
+    virtual HypreParMatrix * Duc(const BlockVector &) = 0;
+    virtual HypreParMatrix * Dmc(const BlockVector &) = 0;
+    virtual void c(const BlockVector &, Vector &) const = 0;
+    int GetDimU() const { return dimU; };
+    int GetDimM() const { return dimM; }; 
+    int GetDimC() const { return dimC; };
+    int GetDimUGlb() const { return dimUglb; };
+    int GetDimMGlb() const { return dimMglb; };
+    HYPRE_BigInt * GetDofOffsetsU() const { return dofOffsetsU; };
+    HYPRE_BigInt * GetDofOffsetsM() const { return dofOffsetsM; }; 
+    Vector Getml() const { return ml; };
+    ~ParGeneralOptProblem();
+};
+
+
+// abstract ContactProblem class
+// of the form
+// min_d e(d) s.t. g(d) >= 0
+class ParOptProblem : public ParGeneralOptProblem
+{
+protected:
+    HypreParMatrix * Ih;
+public:
+    ParOptProblem();
+    void Init(HYPRE_BigInt *, HYPRE_BigInt *);
+    
+    // ParGeneralOptProblem methods are defined in terms of
+    // ParOptProblem specific methods: E, DdE, DddE, g, Ddg
+    double CalcObjective(const BlockVector &) const; 
+    void Duf(const BlockVector &, Vector &) const;
+    void Dmf(const BlockVector &, Vector &) const;
+    HypreParMatrix * Duuf(const BlockVector &);
+    HypreParMatrix * Dumf(const BlockVector &);
+    HypreParMatrix * Dmuf(const BlockVector &);
+    HypreParMatrix * Dmmf(const BlockVector &);
+    void c(const BlockVector &, Vector &) const;
+    HypreParMatrix * Duc(const BlockVector &);
+    HypreParMatrix * Dmc(const BlockVector &);
+    
+    // ParOptProblem specific methods:
+    
+    // energy objective function e(d)
+    // input: d an mfem::Vector
+    // output: e(d) a double
+    virtual double E(const Vector &d) const = 0;
+
+    // gradient of energy objective De / Dd
+    // input: d an mfem::Vector,
+    //        gradE an mfem::Vector, which will be the gradient of E at d
+    // output: none    
+    virtual void DdE(const Vector &d, Vector &gradE) const = 0;
+
+    // Hessian of energy objective D^2 e / Dd^2
+    // input:  d, an mfem::Vector
+    // output: The Hessian of the energy objective at d, a pointer to a HypreParMatrix
+    virtual HypreParMatrix * DddE(const Vector &d) = 0;
+
+    // Constraint function g(d) >= 0, e.g., gap function
+    // input: d, an mfem::Vector,
+    //       gd, an mfem::Vector, which upon successfully calling the g method will be
+    //                            the evaluation of the function g at d
+    // output: none
+    virtual void g(const Vector &d, Vector &gd) const = 0;
+
+    // Jacobian of constraint function Dg / Dd, e.g., gap function Jacobian
+    // input:  d, an mfem::Vector,
+    // output: The Jacobain of the constraint function g at d, a pointer to a HypreParMatrix
+    virtual HypreParMatrix * Ddg(const Vector &) = 0;
+    virtual ~ParOptProblem();
+};
+
+
+
+
+
+
+class ParObstacleProblem : public ParOptProblem
+{
+protected:
+   // data to define energy objective function e(d) = 0.5 d^T K d - f^T d, g(d) = d >= \psi
+   // stiffness matrix used to define objective
+   ParBilinearForm *Kform;
+   ParLinearForm   *fform;
+   Array<int> ess_tdof_list; // needed for calls to FormSystemMatrix
+   HypreParMatrix  K;
+   HypreParMatrix  *J;
+   ParFiniteElementSpace *Vh;
+   Vector f;
+   Vector psi;
+public :
+   ParObstacleProblem(ParFiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &));
+   ParObstacleProblem(ParFiniteElementSpace*, double (*fSource)(const Vector &), double (*obstacleSource)(const Vector &), Array<int> tdof_list, Vector &);
+   double E(const Vector &) const;
+   void DdE(const Vector &, Vector &) const;
+   HypreParMatrix* DddE(const Vector &);
+   void g(const Vector &, Vector &) const;
+   HypreParMatrix* Ddg(const Vector &);
+   virtual ~ParObstacleProblem();
+};
+
+
+
+class ReducedProblem : public ParOptProblem
+{
+protected:
+  HypreParMatrix *J;
+  HypreParMatrix *P; // projector
+  ParOptProblem  *problem;
+public:
+  ReducedProblem(ParOptProblem *problem, HYPRE_Int * constraintMask);
+  ReducedProblem(ParOptProblem *problem, HypreParVector & constraintMask);
+  double E(const Vector &) const;
+  void DdE(const Vector &, Vector &) const;
+  HypreParMatrix * DddE(const Vector &);
+  void g(const Vector &, Vector &) const;
+  HypreParMatrix * Ddg(const Vector &);
+  virtual ~ReducedProblem();
+};
+
+
+
+
+#endif

miniapps/interiorpointsolver/TestProblem1.cpp

@@ -0,0 +1,181 @@
+//                         Example Problem 1
+//
+//
+// Compile with: make ParTestProblem1
+//
+// Sample runs: mpirun -np 4 ./ParTestProblem1
+//
+//
+// Description: This example code demonstrates the use of the MFEM based
+//              interior-point solver to solve the
+//              bound-constrained minimization problem
+//
+//              minimize_(x \in R^n) 1/2 x^T x subject to x - xl ≥ 0 (component-wise).
+//              
+#include "mfem.hpp"
+#include "Problem.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+
+
+// solve min 1/2 x^T K x s.t. J x - xl >= 0
+// where K and J are identity matrices and xl
+// has uniform random values in [-1, 1]
+// for the Lagrangian L(x, s, l, z) = 1/2 x^T x + l^T (x - xl - s) - z^T s
+// the optimal solution is x*_i = max{0, (xl)_i}, z*_i = x*_i
+
+
+class ParEx1Problem : public ParOptProblem
+{
+protected:
+   HypreParMatrix *K;
+   HypreParMatrix *J;
+   Vector xl;
+   //HYPRE_BigInt * dofOffsets;
+public:
+   // create offsets internally only pass problem size
+   //ParEx1Problem(HYPRE_BigInt * offsets);
+   ParEx1Problem(int n);
+   double E(const Vector &) const;
+   void DdE(const Vector &, Vector &) const;
+   HypreParMatrix* DddE(const Vector &);
+   void g(const Vector &, Vector &) const;
+   HypreParMatrix* Ddg(const Vector &);
+   virtual ~ParEx1Problem();
+};
+
+
+
+void mfemIPSolve(ParGeneralOptProblem & problem, Vector &x, Vector &lambda)
+{
+   ParInteriorPointSolver IPoptimizer(&problem);
+
+   int dimX = problem.GetDimU();
+   Vector x0(dimX); x0 = 100.0;
+   x.SetSize(dimX); x = 0.0;
+   
+   double OptTol = 1.e-6;
+   double LinSolveTol = 1.e-10;
+   int linSolveStrategy = 2;
+   int MaxOptIter = 30; 
+   IPoptimizer.SetTol(OptTol);
+   IPoptimizer.SetLinearSolveTol(LinSolveTol);
+   IPoptimizer.SetLinearSolver(linSolveStrategy);
+   IPoptimizer.SetMaxIter(MaxOptIter);
+   IPoptimizer.Mult(x0, x);
+   
+   int dimM = problem.GetDimM();
+   lambda.SetSize(dimM);
+   IPoptimizer.GetLagrangeMultiplier(lambda);
+}
+
+
+
+
+
+
+int main(int argc, char *argv[])
+{
+  // Initialize MPI
+   Mpi::Init();
+   Hypre::Init();   
+
+   int n = 10;
+   OptionsParser args(argc, argv);
+   args.AddOption(&n, "-n", "--n", \
+		   "Size of the optimization problem (dimension of primal variable)");
+   args.ParseCheck();
+   
+   
+   ParEx1Problem problem(n);
+   
+   Vector xOptimal, lambdaOptimal;
+   mfemIPSolve(problem, xOptimal, lambdaOptimal);
+   for(int i = 0; i < xOptimal.Size(); i++)
+   {
+     cout << "optimal (x, z)_" << i << " = (" << xOptimal(i) << ", " << lambdaOptimal(i) << ")\n";
+   }
+
+   Mpi::Finalize();
+   return 0;
+}
+
+
+// Ex1Problem
+// min 1/2 x^T K x such that J x - xl >= 0
+// where K and J are identity matrices
+ParEx1Problem::ParEx1Problem(int n) : ParOptProblem(), K(nullptr), J(nullptr)
+{
+  // generate the parallel partition of the 
+  // variable x and the 
+  int nprocs = Mpi::WorldSize();
+  int myrank = Mpi::WorldRank();
+  
+  
+  HYPRE_BigInt * dofOffsets = new HYPRE_BigInt[2];
+  dofOffsets[0] = HYPRE_BigInt(myrank * n / nprocs);
+  dofOffsets[1] = HYPRE_BigInt((myrank + 1) * n / nprocs);
+  
+  Init(dofOffsets, dofOffsets);
+
+  Vector iDiag(dofOffsets[1] - dofOffsets[0]); iDiag = 1.0;
+  
+  K = GenerateHypreParMatrixFromDiagonal(dofOffsets, iDiag);
+
+  J = GenerateHypreParMatrixFromDiagonal(dofOffsets, iDiag);
+
+  xl.SetSize(dofOffsets[1] - dofOffsets[0]);
+  xl.Randomize(myrank);
+  xl *= 2.0;
+  xl -= 1.0;
+  delete[] dofOffsets;
+}
+
+
+
+double ParEx1Problem::E(const Vector & x) const
+{
+   Vector Kx(K->Height()); Kx = 0.0;
+   MFEM_VERIFY(x.Size() == K->Width(), "ParEx1Problem::E - Inconsistent dimensions");
+   K->Mult(x, Kx);
+   return 0.5 * InnerProduct(MPI_COMM_WORLD, x, Kx);
+}
+
+void ParEx1Problem::DdE(const Vector &x, Vector &gradE) const
+{
+   gradE.SetSize(K->Height());
+   MFEM_VERIFY(x.Size() == K->Width(), "ParEx1Problem::DdE - Inconsistent dimensions");
+   K->Mult(x, gradE);
+}
+
+HypreParMatrix * ParEx1Problem::DddE(const Vector &x)
+{
+   return K; 
+}
+
+// g(x) = x - xl >= 0
+void ParEx1Problem::g(const Vector &x, Vector &gx) const
+{
+   MFEM_VERIFY(x.Size() == J->Width(), "ParEx1Problem::g - Inconsistent dimensions");
+   J->Mult(x, gx);
+   MFEM_VERIFY(gx.Size() == J->Height(), "ParEx1Problem::g - Inconsistent dimensions");
+   gx.Add(-1.0, xl);
+}
+
+HypreParMatrix * ParEx1Problem::Ddg(const Vector &)
+{
+  return J;
+}
+
+ParEx1Problem::~ParEx1Problem()
+{
+   delete K;
+   delete J;
+}
+

miniapps/interiorpointsolver/TestProblem2.cpp

@@ -0,0 +1,161 @@
+//                       Spherical Obstacle Problem
+//
+//
+// Compile with: make ParSphericalObstacleProblem
+//
+// Sample runs: mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 0
+//              mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 1
+//              mpirun -np 4 ./ParSphericalObstacleProblem -linSolver 2
+//
+//
+// Description: This example code demonstrates the use of MFEM to solve the
+//              bound-constrained energy minimization problem
+//
+//                      minimize ||∇u||² subject to u ≥ ϕ in H¹₀.
+
+#include "mfem.hpp"
+#include "Problem.hpp"
+#include "IPsolver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double fRhs(const Vector &);
+double spherical_obstacle(const Vector &);
+double exact_solution_obstacle(const Vector &);
+
+
+int main(int argc, char *argv[])
+{
+  // Initialize MPI
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   int FEorder = 1; // order of the finite elements
+   int linSolver = 2;
+   int maxIPMiters = 30;
+   int ref_levels = 3;
+   OptionsParser args(argc, argv);
+   args.AddOption(&FEorder, "-o", "--order",\
+	 	  "Order of the finite elements.");
+   args.AddOption(&linSolver, "-linSolver", "--linearSolver", \
+        "IP-Newton linear system solution strategy.");
+   args.AddOption(&maxIPMiters, "-IPMiters", "--IPMiters",\
+	 	  "Maximum number of IPM iterations");
+   args.AddOption(&ref_levels, "-r", "--mesh_refinement", \
+		  "Mesh Refinement");
+  
+   args.ParseCheck();
+
+   const char *meshFile = "disk.mesh";
+   Mesh mesh(meshFile, 1, 1);
+   int dim = mesh.Dimension(); // geometric dimension of the meshed domain
+   {
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   
+   FiniteElementCollection *fec = new H1_FECollection(FEorder, dim);
+   ParFiniteElementSpace   *Vh  = new ParFiniteElementSpace(&pmesh, fec);
+   Array<int> boundary_dofs;
+   Vh->GetBoundaryTrueDofs(boundary_dofs);
+   int dimD = Vh->GetTrueVSize();
+   Vector xDC(dimD); xDC = 0.0;
+
+   ParObstacleProblem problem(Vh, &fRhs, &spherical_obstacle, boundary_dofs, xDC);
+   Vector x0(dimD); x0.Set(1.0, xDC);
+   Vector xf(dimD); xf = 0.0;
+
+   ParInteriorPointSolver optimizer(&problem); 
+   optimizer.SetTol(1.e-7);
+   optimizer.SetLinearSolveTol(1.e-9);
+   optimizer.SetLinearSolver(linSolver);
+   optimizer.SetMaxIter(maxIPMiters);
+   optimizer.Mult(x0, xf);
+
+   ParGridFunction d_gf(Vh);
+
+   d_gf.SetFromTrueDofs(xf);
+
+
+   FunctionCoefficient dtrue_fc(exact_solution_obstacle); // analytic solution
+   ParGridFunction dtrue_gf(Vh);
+   dtrue_gf.ProjectCoefficient(dtrue_fc);
+
+   double L2error = d_gf.ComputeL2Error(dtrue_fc);
+   if (myid == 0)
+   {
+     cout << "\n|| u_h - u ||_{L^2} = " << L2error << '\n' << endl;
+   } 
+   
+   ParaViewDataCollection paraview_dc("SphericalObstacleProblem", &pmesh);
+   paraview_dc.SetPrefixPath("ParaView");
+   paraview_dc.SetLevelsOfDetail(FEorder);
+   paraview_dc.SetDataFormat(VTKFormat::BINARY);
+   paraview_dc.SetHighOrderOutput(true);
+   paraview_dc.SetCycle(0);
+   paraview_dc.SetTime(0.0);
+   paraview_dc.RegisterField("u(x,y) (analytic)", &dtrue_gf);
+   paraview_dc.RegisterField("u(x,y) (numerical)", &d_gf);
+   paraview_dc.Save();
+
+   delete Vh;
+   delete fec;
+   return 0;
+}
+
+
+double fRhs(const Vector &x)
+{
+  return 0.;
+}
+
+
+double spherical_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double beta = 0.9;
+
+   double b = r0*beta;
+   double tmp = sqrt(r0*r0 - b*b);
+   double B = tmp + b*b/tmp;
+   double C = -b/tmp;
+
+   if (r > b)
+   {
+      return B + r * C;
+   }
+   else
+   {
+      return sqrt(r0*r0 - r*r);
+   }
+}
+
+double exact_solution_obstacle(const Vector &pt)
+{
+   double x = pt(0), y = pt(1);
+   double r = sqrt(x*x + y*y);
+   double r0 = 0.5;
+   double a =  0.348982574111686;
+   double A = -0.340129705945858;
+
+   if (r > a)
+   {
+      return A * log(r);
+   }
+   else
+   {
+      return sqrt(r0*r0-r*r);
+   }
+} 
+

miniapps/interiorpointsolver/disk.mesh

@@ -0,0 +1,109 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+2
+
+elements
+5
+1 3 4 5 6 7
+1 3 0 1 5 4
+1 3 1 2 6 5
+1 3 3 7 6 2
+1 3 0 4 7 3
+
+boundary
+4
+1 1 0 1
+1 1 2 3
+1 1 1 2
+1 1 3 0
+
+edges
+12
+0 0 1
+0 4 5
+0 7 6
+0 3 2
+1 1 2
+1 5 6
+1 4 7
+1 0 3
+2 0 4
+2 1 5
+2 2 6
+2 3 7
+
+vertices
+8
+
+knotvectors
+3
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+2 3 0 0 0 1 1 1
+
+weights
+1
+1
+1
+1
+1
+1
+1
+1
+0.70710678118655
+1
+1
+0.70710678118655
+0.70710678118655
+1
+1
+0.70710678118655
+1
+1
+1
+1
+1
+0.85355339059327
+0.85355339059327
+0.85355339059327
+0.85355339059327
+
+FiniteElementSpace
+FiniteElementCollection: NURBS2
+VDim: 2
+Ordering: 1
+
+-0.70710678118 -0.70710678118
+0.70710678118 -0.70710678118
+0.70710678118 0.70710678118
+-0.70710678118 0.70710678118
+-0.35355339059 -0.35355339059
+0.35355339059 -0.35355339059
+0.35355339059 0.35355339059
+-0.35355339059 0.35355339059
+0 -1.41421356236
+0 -0.35355339059
+0 0.35355339059
+0 1.41421356236
+1.41421356236 0
+0.35355339059 0
+-0.35355339059 0
+-1.41421356236 0
+-0.530330085885 -0.530330085885
+0.530330085885 -0.530330085885
+0.530330085885 0.530330085885
+-0.530330085885 0.530330085885
+0 0
+0 -0.883883476475
+0.883883476475 0
+0 0.883883476475
+-0.883883476475 0

miniapps/interiorpointsolver/inline-quad.mesh

@@ -0,0 +1,7 @@
+MFEM INLINE mesh v1.0
+
+type = quad
+nx = 4
+ny = 4
+sx = 1.0
+sy = 1.0

miniapps/interiorpointsolver/makefile

@@ -0,0 +1,35 @@
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/interiorpointsolver/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+# Remove built-in rule
+%: %.cpp
+%: %.c
+%: %.o
+
+%.o: %.cpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
+
+FRAMEWORK_SRC = Problem.cpp IPsolver.cpp utilities.cpp
+FRAMEWORK_OBJ = $(FRAMEWORK_SRC:.cpp=.o)
+
+TestProblem1: TestProblem1.o $(FRAMEWORK_OBJ) $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) TestProblem1.o $(FRAMEWORK_OBJ) -o $@ $(MFEM_LIBS)
+
+TestProblem2: TestProblem2.o $(FRAMEWORK_OBJ) $(MFEM_LIB_FILE)
+	$(MFEM_CXX) $(MFEM_FLAGS) TestProblem2.o $(FRAMEWORK_OBJ) -o $@ $(MFEM_LIBS)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+
+.PHONY: clean
+clean:
+	rm -f *.o TestProblem1 TestProblem2

miniapps/interiorpointsolver/utilities.cpp

@@ -0,0 +1,178 @@
+#include "mfem.hpp"
+#include "utilities.hpp"
+
+using namespace mfem;
+
+HypreParMatrix * GenerateHypreParMatrixFromSparseMatrix(HYPRE_BigInt * colOffsetsloc, HYPRE_BigInt * rowOffsetsloc, SparseMatrix * Asparse)
+{
+  int ncols_loc = colOffsetsloc[1] - colOffsetsloc[0];
+  int nrows_loc = rowOffsetsloc[1] - rowOffsetsloc[0];
+  HYPRE_BigInt ncols_glb, nrows_glb;
+  MPI_Allreduce(&nrows_loc, &nrows_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&ncols_loc, &ncols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+  int * AI          = Asparse->GetI();
+  HYPRE_BigInt * AJ = Asparse->GetJ();
+  double * Adata    = Asparse->GetData();
+
+  HypreParMatrix * Ahypre = nullptr;
+  Ahypre = new HypreParMatrix(MPI_COMM_WORLD, nrows_loc, nrows_glb, ncols_glb, AI, AJ, Adata, rowOffsetsloc, colOffsetsloc);
+  return Ahypre;
+}
+
+
+HypreParMatrix * GenerateHypreParMatrixFromDiagonal(HYPRE_BigInt * offsetsloc, 
+		Vector & diag)
+{
+   int n_loc = offsetsloc[1] - offsetsloc[0];
+   int n_glb = 0;
+   MPI_Allreduce(&n_loc, &n_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+   
+   SparseMatrix * Dsparse = new SparseMatrix(n_loc, n_glb);
+   Array<int> cols;
+   Vector entries;
+   cols.SetSize(1);
+   entries.SetSize(1);
+   for(int j = 0; j < n_loc; j++)
+   {
+     cols[0] = offsetsloc[0] + j;
+     entries(0) = diag(j);
+     Dsparse->SetRow(j, cols, entries);
+   }   
+   Dsparse->Finalize();
+   HypreParMatrix * Dhypre = nullptr;
+   Dhypre = GenerateHypreParMatrixFromSparseMatrix(offsetsloc, offsetsloc, Dsparse);
+   delete Dsparse;
+   return Dhypre;   
+}
+
+HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, HYPRE_Int * mask)
+{
+  int n_cols_loc = offsets[1] - offsets[0];
+  int n_cols_glb = 0;
+  MPI_Allreduce(&n_cols_loc, &n_cols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+  int n_rows_loc = reduced_offsets[1] - reduced_offsets[0];
+
+  SparseMatrix * Psparse = new SparseMatrix(n_rows_loc, n_cols_glb);
+  Array<int> cols;
+  Vector entries;
+  cols.SetSize(1);
+  entries.SetSize(1);
+
+  int row = 0;
+  for(int j = 0; j < n_cols_loc; j++)
+  {
+    if (mask[j] == 1)
+    {
+      cols[0] = offsets[0] + j;
+      entries(0) = 1.0;
+      Psparse->SetRow(row, cols, entries);
+      row += 1;
+    }
+  }
+  Psparse->Finalize();
+  HypreParMatrix * Phypre = nullptr;
+  Phypre = GenerateHypreParMatrixFromSparseMatrix(offsets, reduced_offsets, Psparse);
+  delete Psparse;
+  return Phypre;
+}
+
+HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, const HypreParVector & mask)
+{
+  int n_cols_loc = offsets[1] - offsets[0];
+  int n_cols_glb = 0;
+  MPI_Allreduce(&n_cols_loc, &n_cols_glb, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+  int n_rows_loc = reduced_offsets[1] - reduced_offsets[0];
+
+  SparseMatrix * Psparse = new SparseMatrix(n_rows_loc, n_cols_glb);
+  Array<int> cols;
+  Vector entries;
+  cols.SetSize(1);
+  entries.SetSize(1);
+
+  int row = 0;
+  for(int j = 0; j < n_cols_loc; j++)
+  {
+    if (mask(j) > 0.5)
+    {
+      cols[0] = offsets[0] + j;
+      entries(0) = 1.0;
+      Psparse->SetRow(row, cols, entries);
+      row += 1;
+    }
+  }
+  Psparse->Finalize();
+  HypreParMatrix * Phypre = nullptr;
+  Phypre = GenerateHypreParMatrixFromSparseMatrix(offsets, reduced_offsets, Psparse);
+  delete Psparse;
+  return Phypre;
+}
+
+
+
+
+HYPRE_BigInt * offsetsFromLocalSizes(int n)
+{
+  HYPRE_BigInt * offsets = new HYPRE_BigInt[2];
+  
+  int nprocs = Mpi::WorldSize();
+  int myrank = Mpi::WorldRank();
+  
+  if (myrank == 0)
+  {
+    offsets[0] = 0;
+    offsets[1] = n;
+  }
+  else
+  {
+    offsets[0] = 0;
+    offsets[1] = 0;
+  }
+  
+  // receive then send
+  
+  // Receive local size info from processes with rank less than myrank 
+  // Populate that as entries of helper
+  HYPRE_BigInt * helper;
+  if (myrank > 0)
+  {
+    helper = new HYPRE_BigInt[myrank];
+  }
+  int tag;
+  for (int i = 0; i < myrank; i++)
+  {
+    tag = myrank + i * nprocs;
+    MPI_Recv (&(helper[i]), 1, MPI_INT, i, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    offsets[0] += helper[i];
+  }
+
+  if (myrank > 0)
+  {
+    delete[] helper;
+  }
+  offsets[1] = offsets[0] + n;
+  
+  // Send local size info to all processes with rank greater than myrank
+  for (int i = myrank + 1; i < nprocs; i++)
+  {
+    tag = i + myrank * nprocs;
+    MPI_Send (&n, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
+  }
+  return offsets;
+}
+
+
+void HypreToMfemOffsets(HYPRE_BigInt * offsets)
+{
+  if (offsets[1] < offsets[0])
+  {
+    offsets[1] = offsets[0];
+  }
+  else
+  {
+    offsets[1] = offsets[1] + 1;
+  }
+}
+

miniapps/interiorpointsolver/utilities.hpp

@@ -0,0 +1,23 @@
+#include "mfem.hpp"
+
+using namespace mfem;
+#ifndef UTILITY_FUNCTIONS
+#define UTILITY_FUNCTIONS
+
+void HypreToMfemOffsets(HYPRE_BigInt * offsets);
+
+HypreParMatrix * GenerateHypreParMatrixFromSparseMatrix(HYPRE_BigInt * colOffsetsloc, HYPRE_BigInt * rowOffsetsloc, SparseMatrix * Asparse);
+
+HypreParMatrix * GenerateHypreParMatrixFromDiagonal(HYPRE_BigInt * offsetsloc, 
+		Vector & diag);
+
+
+HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, HYPRE_Int * mask);
+
+HypreParMatrix * GenerateProjector(HYPRE_BigInt * offsets, HYPRE_BigInt * reduced_offsets, const HypreParVector & mask);
+
+
+HYPRE_BigInt * offsetsFromLocalSizes(int n);
+
+
+#endif