text:
/* multilarge.c
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multilarge.h>
#include <gsl/gsl_blas.h>

gsl_multilarge_linear_workspace *
gsl_multilarge_linear_alloc(const gsl_multilarge_linear_type *T,
                            const size_t p)
{
  gsl_multilarge_linear_workspace *w;

  w = calloc(1, sizeof(gsl_multilarge_linear_workspace));
  if (w == NULL)
    {
      GSL_ERROR_NULL("failed to allocate space for workspace",
                     GSL_ENOMEM);
    }

  w->type = T;

  w->state = w->type->alloc(p);
  if (w->state == NULL)
    {
      gsl_multilarge_linear_free(w);
      GSL_ERROR_NULL("failed to allocate space for multilarge state",
                     GSL_ENOMEM);
    }

  w->p = p;

  /* initialize newly allocated state */
  gsl_multilarge_linear_reset(w);

  return w;
}

void
gsl_multilarge_linear_free(gsl_multilarge_linear_workspace *w)
{
  RETURN_IF_NULL(w);

  if (w->state)
    w->type->free(w->state);

  free(w);
}

const char *
gsl_multilarge_linear_name(const gsl_multilarge_linear_workspace *w)
{
  return w->type->name;
}

int
gsl_multilarge_linear_reset(gsl_multilarge_linear_workspace *w)
{
  int status = w->type->reset(w->state);
  return status;
}

int
gsl_multilarge_linear_accumulate(gsl_matrix * X, gsl_vector * y,
                                 gsl_multilarge_linear_workspace * w)
{
  int status = w->type->accumulate(X, y, w->state);
  return status;
}

int
gsl_multilarge_linear_solve(const double lambda, gsl_vector * c,
                            double * rnorm, double * snorm,
                            gsl_multilarge_linear_workspace * w)
{
  int status = w->type->solve(lambda, c, rnorm, snorm, w->state);
  return status;
}

int
gsl_multilarge_linear_rcond(double *rcond, gsl_multilarge_linear_workspace * w)
{
  int status = w->type->rcond(rcond, w->state);
  return status;
}

int
gsl_multilarge_linear_lcurve(gsl_vector * reg_param, gsl_vector * rho,
                             gsl_vector * eta,
                             gsl_multilarge_linear_workspace * w)
{
  const size_t len = reg_param->size;

  if (len != rho->size)
    {
      GSL_ERROR ("reg_param and rho have different sizes", GSL_EBADLEN);
    }
  else if (len != eta->size)
    {
      GSL_ERROR ("reg_param and eta have different sizes", GSL_EBADLEN);
    }
  else
    {
      int status = w->type->lcurve(reg_param, rho, eta, w->state);
      return status;
    }
}

/*
gsl_multilarge_linear_wstdform1()
  Using regularization matrix
L = diag(l_1,l_2,...,l_p), transform to Tikhonov standard form:

X~ = sqrt(W) X L^{-1}
y~ = sqrt(W) y
c~ = L c

Inputs: L    - Tikhonov matrix as a vector of diagonal elements p-by-1;
               or NULL for L = I
        X    - least squares matrix n-by-p
        y    - right hand side vector n-by-1
        w    - weight vector n-by-1; or NULL for W = I
        Xs   - least squares matrix in standard form X~ n-by-p
        ys   - right hand side vector in standard form y~ n-by-1
        work - workspace

Return: success/error

Notes:
1) It is allowed for X = Xs and y = ys
*/

int
gsl_multilarge_linear_wstdform1 (const gsl_vector * L,
                                 const gsl_matrix * X,
                                 const gsl_vector * w,
                                 const gsl_vector * y,
                                 gsl_matrix * Xs,
                                 gsl_vector * ys,
                                 gsl_multilarge_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  (void) work;

  if (L != NULL && p != L->size)
    {
      GSL_ERROR("L vector does not match X", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("weight vector does not match X", GSL_EBADLEN);
    }
  else if (n != Xs->size1 || p != Xs->size2)
    {
      GSL_ERROR("Xs matrix dimensions do not match X", GSL_EBADLEN);
    }
  else if (n != ys->size)
    {
      GSL_ERROR("ys vector must be length n", GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;

      /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
      status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
      if (status)
        return status;

      if (L != NULL)
        {
          size_t j;

          /* construct X~ = sqrt(W) X * L^{-1} matrix */
          for (j = 0; j < p; ++j)
            {
              gsl_vector_view Xj = gsl_matrix_column(Xs, j);
              double lj = gsl_vector_get(L, j);

              if (lj == 0.0)
                {
                  GSL_ERROR("L matrix is singular", GSL_EDOM);
                }

              gsl_vector_scale(&Xj.vector, 1.0 / lj);
            }
        }

      return status;
    }
}

int
gsl_multilarge_linear_stdform1 (const gsl_vector * L,
                                const gsl_matrix * X,
                                const gsl_vector * y,
                                gsl_matrix * Xs,
                                gsl_vector * ys,
                                gsl_multilarge_linear_workspace * work)
{
  int status;

  status = gsl_multilarge_linear_wstdform1(L, X, NULL, y, Xs, ys, work);

  return status;
}

int
gsl_multilarge_linear_L_decomp (gsl_matrix * L, gsl_vector * tau)
{
  const size_t m = L->size1;
  const size_t p = L->size2;

  if (m < p)
    {
      GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
    }
  else
    {
      int status;

      status = gsl_multifit_linear_L_decomp(L, tau);

      return status;
    }
}

int
gsl_multilarge_linear_wstdform2 (const gsl_matrix * LQR,
                                 const gsl_vector * Ltau,
                                 const gsl_matrix * X,
                                 const gsl_vector * w,
                                 const gsl_vector * y,
                                 gsl_matrix * Xs,
                                 gsl_vector * ys,
                                 gsl_multilarge_linear_workspace * work)
{
  const size_t m = LQR->size1;
  const size_t n = X->size1;
  const size_t p = X->size2;

  (void) Ltau;

  if (p != work->p)
    {
      GSL_ERROR("X has wrong number of columns", GSL_EBADLEN);
    }
  else if (p != LQR->size2)
    {
      GSL_ERROR("LQR and X matrices have different numbers of columns", GSL_EBADLEN);
    }
  else if (n != y->size)
    {
      GSL_ERROR("y vector does not match X", GSL_EBADLEN);
    }
  else if (w != NULL && n != w->size)
    {
      GSL_ERROR("weights vector must be length n", GSL_EBADLEN);
    }
  else if (m < p)
    {
      GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
    }
  else if (n != Xs->size1 || p != Xs->size2)
    {
      GSL_ERROR("Xs matrix must be n-by-p", GSL_EBADLEN);
    }
  else if (n != ys->size)
    {
      GSL_ERROR("ys vector must have length n", GSL_EBADLEN);
    }
  else
    {
      int status;
      size_t i;
      gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p);

      /* compute Xs = sqrt(W) X and ys = sqrt(W) y */
      status = gsl_multifit_linear_applyW(X, w, y, Xs, ys);
      if (status)
        return status;

      /* compute X~ = X R^{-1} using QR decomposition of L */
      for (i = 0; i < n; ++i)
        {
          gsl_vector_view v = gsl_matrix_row(Xs, i);

          /* solve: R^T y = X_i */
          gsl_blas_dtrsv(CblasUpper, CblasTrans, CblasNonUnit, &R.matrix, &v.vector);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_multilarge_linear_stdform2 (const gsl_matrix * LQR,
                                const gsl_vector * Ltau,
                                const gsl_matrix * X,
                                const gsl_vector * y,
                                gsl_matrix * Xs,
                                gsl_vector * ys,
                                gsl_multilarge_linear_workspace * work)
{
  int status;

  status = gsl_multilarge_linear_wstdform2(LQR, Ltau, X, NULL, y, Xs, ys, work);

  return status;
}

/*
gsl_multilarge_linear_genform1()
  Backtransform regularized solution vector using matrix
L = diag(L)
*/

int
gsl_multilarge_linear_genform1 (const gsl_vector * L,
                                const gsl_vector * cs,
                                gsl_vector * c,
                                gsl_multilarge_linear_workspace * work)
{
  if (L->size != work->p)
    {
      GSL_ERROR("L vector does not match workspace", GSL_EBADLEN);
    }
  else if (L->size != cs->size)
    {
      GSL_ERROR("cs vector does not match L", GSL_EBADLEN);
    }
  else if (L->size != c->size)
    {
      GSL_ERROR("c vector does not match L", GSL_EBADLEN);
    }
  else
    {
      /* compute true solution vector c = L^{-1} c~ */
      gsl_vector_memcpy(c, cs);
      gsl_vector_div(c, L);

      return GSL_SUCCESS;
    }
}

int
gsl_multilarge_linear_genform2 (const gsl_matrix * LQR,
                                const gsl_vector * Ltau,
                                const gsl_vector * cs,
                                gsl_vector * c,
                                gsl_multilarge_linear_workspace * work)
{
  const size_t m = LQR->size1;
  const size_t p = LQR->size2;

  (void) Ltau;
  (void) work;

  if (p != c->size)
    {
      GSL_ERROR("c vector does not match LQR", GSL_EBADLEN);
    }
  else if (m < p)
    {
      GSL_ERROR("m < p not yet supported", GSL_EBADLEN);
    }
  else if (p != cs->size)
    {
      GSL_ERROR("cs vector size does not match c", GSL_EBADLEN);
    }
  else
    {
      int s;
      gsl_matrix_const_view R = gsl_matrix_const_submatrix(LQR, 0, 0, p, p); /* R factor of L */

      /* solve R c = cs for true solution c, using QR decomposition of L */
      gsl_vector_memcpy(c, cs);
      s = gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, &R.matrix, c);

      return s;
    }
}

const gsl_matrix *
gsl_multilarge_linear_matrix_ptr (const gsl_multilarge_linear_workspace * work)
{
  return work->type->matrix_ptr(work->state);
}

const gsl_vector *
gsl_multilarge_linear_rhs_ptr (const gsl_multilarge_linear_workspace * work)
{
  return work->type->rhs_ptr(work->state);
}


text:
#if !defined(_MUSIC_GSL_)
#define _MUSIC_GSL_

#define WIN32

#include <gsl/gsl_combination.h>
#include <gsl/gsl_sf.h>

#endif

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/* siman/gsl_siman.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_SIMAN_H__
#define __GSL_SIMAN_H__
#include <stdlib.h>
#include <gsl/gsl_rng.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* types for the function pointers passed to gsl_siman_solve */

typedef double (*gsl_siman_Efunc_t) (void *xp);
typedef void (*gsl_siman_step_t) (const gsl_rng *r, void *xp, double step_size);
typedef double (*gsl_siman_metric_t) (void *xp, void *yp);
typedef void (*gsl_siman_print_t) (void *xp);
typedef void (*gsl_siman_copy_t) (void *source, void *dest);
typedef void * (*gsl_siman_copy_construct_t) (void *xp);
typedef void (*gsl_siman_destroy_t) (void *xp);

/* this structure contains all the information needed to structure the
   search, beyond the energy function, the step function and the
   initial guess. */

typedef struct {
  int n_tries;          /* how many points to try for each step */
  int iters_fixed_T;    /* how many iterations at each temperature? */
  double step_size;     /* max step size in the random walk */
  /* the following parameters are for the Boltzmann distribution */
  double k, t_initial, mu_t, t_min;
} gsl_siman_params_t;

/* prototype for the workhorse function */

void gsl_siman_solve(const gsl_rng * r,
                     void *x0_p, gsl_siman_Efunc_t Ef,
                     gsl_siman_step_t take_step,
                     gsl_siman_metric_t distance,
                     gsl_siman_print_t print_position,
                     gsl_siman_copy_t copyfunc,
                     gsl_siman_copy_construct_t copy_constructor,
                     gsl_siman_destroy_t destructor,
                     size_t element_size,
                     gsl_siman_params_t params);

void
gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
                      gsl_siman_step_t take_step,
                      gsl_siman_metric_t distance,
                      gsl_siman_print_t print_position,
                      size_t element_size,
                      gsl_siman_params_t params);

__END_DECLS

#endif /* __GSL_SIMAN_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
static char help[] =
"Solves structured-grid Poisson problem in 1D, 2D, 3D.  Option prefix fsh_.\n"
"Equation is\n"
"    - cx u_xx - cy u_yy - cz u_zz = f,\n"
"subject to Dirichlet boundary conditions.  Solves three different problems\n"
"where exact solution is known.  Uses DMDA and SNES.  Equation is put in form\n"
"F(u) = - grad^2 u - f.  Call-backs fully-rediscretize for the supplied grid.\n"
"Defaults to 2D, a SNESType of KSPONLY, and a KSPType of CG.\n\n";

#include <petsc.h>
#include "poissonfunctions.h"

// exact solutions  u(x,y),  for boundary condition and error calculation

static double u_exact_1Dmanupoly(double x, double y, double z, void *ctx) {
    return x*x * (1.0 - x*x);
}

static double u_exact_2Dmanupoly(double x, double y, double z, void *ctx) {
    return x*x * (1.0 - x*x) * y*y *(y*y - 1.0);
}

static double u_exact_3Dmanupoly(double x, double y, double z, void *ctx) {
    return x*x * (1.0 - x*x) * y*y * (y*y - 1.0) * z*z * (z*z - 1.0);
}

static double u_exact_1Dmanuexp(double x, double y, double z, void *ctx) {
    return - exp(x);
}

static double u_exact_2Dmanuexp(double x, double y, double z, void *ctx) {
    return - x * exp(y);
}

static double u_exact_3Dmanuexp(double x, double y, double z, void *ctx) {
    return - x * exp(y + z);
}

static double zero(double x, double y, double z, void *ctx) {
    return 0.0;
}

// right-hand-side functions  f(x,y) = - laplacian u

static double f_rhs_1Dmanupoly(double x, double y, double z, void *ctx) {
    PoissonCtx* user = (PoissonCtx*)ctx;
    return user->cx * 12.0 * x*x - 2.0;
}

static double f_rhs_2Dmanupoly(double x, double y, double z, void *ctx) {
    PoissonCtx* user = (PoissonCtx*)ctx;
    double aa, bb, ddaa, ddbb;
    aa = x*x * (1.0 - x*x);
    bb = y*y * (y*y - 1.0);
    ddaa = 2.0 * (1.0 - 6.0 * x*x);
    ddbb = 2.0 * (6.0 * y*y - 1.0);
    return - (user->cx * ddaa * bb + user->cy * aa * ddbb);
}

static double f_rhs_3Dmanupoly(double x, double y, double z, void *ctx) {
    PoissonCtx* user = (PoissonCtx*)ctx;
    double aa, bb, cc, ddaa, ddbb, ddcc;
    aa = x*x * (1.0 - x*x);
    bb = y*y * (y*y - 1.0);
    cc = z*z * (z*z - 1.0);
    ddaa = 2.0 * (1.0 - 6.0 * x*x);
    ddbb = 2.0 * (6.0 * y*y - 1.0);
    ddcc = 2.0 * (6.0 * z*z - 1.0);
    return - (user->cx * ddaa * bb * cc + user->cy * aa * ddbb * cc + user->cz * aa * bb * ddcc);
}

static double f_rhs_1Dmanuexp(double x, double y, double z, void *ctx) {
    return exp(x);
}

static double f_rhs_2Dmanuexp(double x, double y, double z, void *ctx) {
    return x * exp(y);  // note  f = - (u_xx + u_yy) = - u
}

static double f_rhs_3Dmanuexp(double x, double y, double z, void *ctx) {
    return 2.0 * x * exp(y + z);  // note  f = - laplacian u = - 2 u
}

// functions simply to put u_exact()=g_bdry() into a grid
// these are irritatingly-dimension-dependent inside ...
extern PetscErrorCode Form1DUExact(DMDALocalInfo*, Vec, PoissonCtx*);
extern PetscErrorCode Form2DUExact(DMDALocalInfo*, Vec, PoissonCtx*);
extern PetscErrorCode Form3DUExact(DMDALocalInfo*, Vec, PoissonCtx*);

//STARTPTRARRAYS
// arrays of pointers to functions:   ..._ptr[DIMS]
static void* residual_ptr[3]
    = {&Poisson1DFunctionLocal, &Poisson2DFunctionLocal, &Poisson3DFunctionLocal};

static void* jacobian_ptr[3]
    = {&Poisson1DJacobianLocal, &Poisson2DJacobianLocal, &Poisson3DJacobianLocal};

static void* getuexact_ptr[3]
    = {&Form1DUExact, &Form2DUExact, &Form3DUExact};
//ENDPTRARRAYS

typedef enum {MANUPOLY, MANUEXP, ZERO} ProblemType;
static const char* ProblemTypes[] = {"manupoly","manuexp","zero",
                                     "ProblemType", "", NULL};

// more arrays of pointers to functions:   ..._ptr[DIMS][PROBLEMS]
static void* g_bdry_ptr[3][3]
    = {{&u_exact_1Dmanupoly, &u_exact_1Dmanuexp, &zero},
       {&u_exact_2Dmanupoly, &u_exact_2Dmanuexp, &zero},
       {&u_exact_3Dmanupoly, &u_exact_3Dmanuexp, &zero}};

static void* f_rhs_ptr[3][3]
    = {{&f_rhs_1Dmanupoly, &f_rhs_1Dmanuexp, &zero},
       {&f_rhs_2Dmanupoly, &f_rhs_2Dmanuexp, &zero},
       {&f_rhs_3Dmanupoly, &f_rhs_3Dmanuexp, &zero}};

static const char* InitialTypes[] = {"zeros","random",
                                     "InitialType", "", NULL};

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    DM             da, da_after;
    SNES           snes;
    KSP            ksp;
    Vec            u_initial, u, u_exact;
    PoissonCtx     user;
    DMDALocalInfo  info;
    double         errinf, normconst2h, err2h;
    char           gridstr[99];
    PetscErrorCode (*getuexact)(DMDALocalInfo*,Vec,PoissonCtx*);

    // fish defaults:
    int            dim = 2;                  // 2D
    ProblemType    problem = MANUEXP;        // manufactured problem using exp()
    InitialType    initial = ZEROS;          // set u=0 for initial iterate
    PetscBool      gonboundary = PETSC_TRUE; // initial iterate has u=g on boundary

    // get options and configure context
    user.Lx = 1.0;
    user.Ly = 1.0;
    user.Lz = 1.0;
    user.cx = 1.0;
    user.cy = 1.0;
    user.cz = 1.0;
    PetscInitialize(&argc,&argv,NULL,help);
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"fsh_", "options for fish.c", ""); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-cx",
         "set coefficient of x term u_xx in equation",
         "fish.c",user.cx,&user.cx,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-cy",
         "set coefficient of y term u_yy in equation",
         "fish.c",user.cy,&user.cy,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-cz",
         "set coefficient of z term u_zz in equation",
         "fish.c",user.cz,&user.cz,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsInt("-dim",
         "dimension of problem (=1,2,3 only)",
         "fish.c",dim,&dim,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsBool("-initial_gonboundary",
         "set initial iterate to have correct boundary values",
         "fish.c",gonboundary,&gonboundary,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-initial_type",
         "type of initial iterate",
         "fish.c",InitialTypes,(PetscEnum)initial,(PetscEnum*)&initial,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsReal("-Lx",
         "set Lx in domain ([0,Lx] x [0,Ly] x [0,Lz], etc.)",
         "fish.c",user.Lx,&user.Lx,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-Ly",
         "set Ly in domain ([0,Lx] x [0,Ly] x [0,Lz], etc.)",
         "fish.c",user.Ly,&user.Ly,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsReal("-Lz",
         "set Ly in domain ([0,Lx] x [0,Ly] x [0,Lz], etc.)",
         "fish.c",user.Lz,&user.Lz,NULL);CHKERRQ(ierr);
    ierr = PetscOptionsEnum("-problem",
         "problem type; determines exact solution and RHS",
         "fish.c",ProblemTypes,(PetscEnum)problem,(PetscEnum*)&problem,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
    user.g_bdry = g_bdry_ptr[dim-1][problem];
    user.f_rhs = f_rhs_ptr[dim-1][problem];
    if ( user.cx <= 0.0 || user.cy <= 0.0 || user.cz <= 0.0 ) {
        SETERRQ(PETSC_COMM_WORLD,2,"positivity required for coefficients cx,cy,cz\n");
    }
    if ((problem == MANUEXP) && ( user.cx != 1.0 || user.cy != 1.0 || user.cz != 1.0)) {
        SETERRQ(PETSC_COMM_WORLD,3,"cx=cy=cz=1 required for problem MANUEXP\n");
    }

//STARTCREATE
    // DMDA in chosen dimension
    switch (dim) {
        case 1:
            ierr = DMDACreate1d(PETSC_COMM_WORLD,
                DM_BOUNDARY_NONE,3,1,1, NULL, &da); CHKERRQ(ierr);
            break;
        case 2:
            ierr = DMDACreate2d(PETSC_COMM_WORLD,
                DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,
                3,3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
            break;
        case 3:
            ierr = DMDACreate3d(PETSC_COMM_WORLD,
                DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,
                DMDA_STENCIL_STAR,
                3,3,3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
                1,1,NULL,NULL,NULL,&da); CHKERRQ(ierr);
            break;
        default:
            SETERRQ(PETSC_COMM_WORLD,1,"invalid dim for DMDA creation\n");
    }
    ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);  // call BEFORE SetUniformCoordinates
    ierr = DMDASetUniformCoordinates(da,0.0,user.Lx,0.0,user.Ly,0.0,user.Lz); CHKERRQ(ierr);

    // create SNES supply call-backs
    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
    ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)(residual_ptr[dim-1]),&user); CHKERRQ(ierr);
    ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)(jacobian_ptr[dim-1]),&user); CHKERRQ(ierr);

    // default to KSPONLY+CG because problem is linear and SPD
    ierr = SNESSetType(snes,SNESKSPONLY); CHKERRQ(ierr);
    ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
    ierr = KSPSetType(ksp,KSPCG); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

    // initial iterate and solve
    ierr = DMGetGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = InitialState(da, initial, gonboundary, u_initial, &user); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,u_initial); CHKERRQ(ierr);
//ENDCREATE

//STARTGETSOLUTION
    // -snes_grid_sequence could change grid resolution
    ierr = DMRestoreGlobalVector(da,&u_initial); CHKERRQ(ierr);
    ierr = DMDestroy(&da); CHKERRQ(ierr);

    // evaluate error and report
    ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr);  // SNES owns u; do not destroy it
    ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
    ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr); // SNES owns da_after; do not destroy it
    ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
    getuexact = getuexact_ptr[dim-1];
    ierr = (*getuexact)(&info,u_exact,&user); CHKERRQ(ierr);
    ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr);   // u <- u + (-1.0) uexact
    ierr = VecDestroy(&u_exact); CHKERRQ(ierr);      // no longer needed
    ierr = VecNorm(u,NORM_INFINITY,&errinf); CHKERRQ(ierr);
    ierr = VecNorm(u,NORM_2,&err2h); CHKERRQ(ierr);
//ENDGETSOLUTION

    switch (dim) {
        case 1:
            normconst2h = PetscSqrtReal((double)(info.mx-1));
            snprintf(gridstr,99,"%d point 1D",info.mx);
            break;
        case 2:
            normconst2h = PetscSqrtReal((double)(info.mx-1)*(info.my-1));
            snprintf(gridstr,99,"%d x %d point 2D",info.mx,info.my);
            break;
        case 3:
            normconst2h = PetscSqrtReal((double)(info.mx-1)*(info.my-1)*(info.mz-1));
            snprintf(gridstr,99,"%d x %d x %d point 3D",info.mx,info.my,info.mz);
            break;
        default:
            SETERRQ(PETSC_COMM_WORLD,4,"invalid dim value in final report\n");
    }
    err2h /= normconst2h; // like continuous L2
    ierr = PetscPrintf(PETSC_COMM_WORLD,
                "problem %s on %s grid:\n"
                "  error |u-uexact|_inf = %.3e, |u-uexact|_h = %.3e\n",
                ProblemTypes[problem],gridstr,errinf,err2h); CHKERRQ(ierr);

    // destroy what we explicitly Created
    ierr = SNESDestroy(&snes); CHKERRQ(ierr);
    return PetscFinalize();
}

PetscErrorCode Form1DUExact(DMDALocalInfo *info, Vec u, PoissonCtx* user) {
  PetscErrorCode ierr;
  int          i;
  double       xmax[1], xmin[1], hx, x, *au;
  ierr = DMDAGetBoundingBox(info->da,xmin,xmax); CHKERRQ(ierr);
  hx = (xmax[0] - xmin[0]) / (info->mx - 1);
  ierr = DMDAVecGetArray(info->da, u, &au);CHKERRQ(ierr);
  for (i=info->xs; i<info->xs+info->xm; i++) {
      x = xmin[0] + i * hx;
      au[i] = user->g_bdry(x,0.0,0.0,user);
  }
  ierr = DMDAVecRestoreArray(info->da, u, &au);CHKERRQ(ierr);
  return 0;
}

PetscErrorCode Form2DUExact(DMDALocalInfo *info, Vec u, PoissonCtx* user) {
    PetscErrorCode ierr;
    int     i, j;
    double  xymin[2], xymax[2], hx, hy, x, y, **au;
    ierr = DMDAGetBoundingBox(info->da,xymin,xymax); CHKERRQ(ierr);
    hx = (xymax[0] - xymin[0]) / (info->mx - 1);
    hy = (xymax[1] - xymin[1]) / (info->my - 1);
    ierr = DMDAVecGetArray(info->da, u, &au);CHKERRQ(ierr);
    for (j=info->ys; j<info->ys+info->ym; j++) {
        y = xymin[1] + j * hy;
        for (i=info->xs; i<info->xs+info->xm; i++) {
            x = xymin[0] + i * hx;
            au[j][i] = user->g_bdry(x,y,0.0,user);
        }
    }
    ierr = DMDAVecRestoreArray(info->da, u, &au);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode Form3DUExact(DMDALocalInfo *info, Vec u, PoissonCtx* user) {
    PetscErrorCode ierr;
    int    i, j, k;
    double xyzmin[3], xyzmax[3], hx, hy, hz, x, y, z, ***au;
    ierr = DMDAGetBoundingBox(info->da,xyzmin,xyzmax); CHKERRQ(ierr);
    hx = (xyzmax[0] - xyzmin[0]) / (info->mx - 1);
    hy = (xyzmax[1] - xyzmin[1]) / (info->my - 1);
    hz = (xyzmax[2] - xyzmin[2]) / (info->mz - 1);
    ierr = DMDAVecGetArray(info->da, u, &au);CHKERRQ(ierr);
    for (k=info->zs; k<info->zs+info->zm; k++) {
        z = xyzmin[2] + k * hz;
        for (j=info->ys; j<info->ys+info->ym; j++) {
            y = xyzmin[1] + j * hy;
            for (i=info->xs; i<info->xs+info->xm; i++) {
                x = xyzmin[0] + i * hx;
                au[k][j][i] = user->g_bdry(x,y,z,user);
            }
        }
    }
    ierr = DMDAVecRestoreArray(info->da, u, &au);CHKERRQ(ierr);
    return 0;
}

ext: c