text:
/* bst/gsl_bst.h
 *
 * Copyright (C) 2018, 2019 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.
 */

#ifndef __GSL_BST_H__
#define __GSL_BST_H__

#include <gsl/gsl_math.h>
#include <gsl/gsl_bst_avl.h>
#include <gsl/gsl_bst_rb.h>
#include <gsl/gsl_bst_types.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

/* type of binary search tree */
typedef struct
{
  const char *name;
  const size_t node_size;
  int    (*init)(const gsl_bst_allocator * allocator,
                 gsl_bst_cmp_function * compare, void * params, void * vtable);
  size_t (*nodes) (const void * vtable);
  void * (*insert) (void * item, void * vtable);
  void * (*find) (const void *item, const void * vtable);
  void * (*remove) (const void *item, void * vtable);
  int    (*empty) (void * vtable);
  int    (*trav_init) (void * vtrav, const void * vtable);
  void * (*trav_first) (void * vtrav, const void * vtable);
  void * (*trav_last) (void * vtrav, const void * vtable);
  void * (*trav_find) (const void * item, void * vtrav, const void * vtable);
  void * (*trav_insert) (void * item, void * vtrav, void * vtable);
  void * (*trav_copy) (void * vtrav, const void * vsrc);
  void * (*trav_next) (void * vtrav);
  void * (*trav_prev) (void * vtrav);
  void * (*trav_cur) (const void * vtrav);
  void * (*trav_replace) (void * vtrav, void * new_item);
} gsl_bst_type;

typedef struct
{
  const gsl_bst_type * type;    /* binary search tree type */
  union
    {
      gsl_bst_avl_table avl_table;
      gsl_bst_rb_table rb_table;
    } table;
} gsl_bst_workspace;

typedef struct
{
  const gsl_bst_type * type;    /* binary search tree type */
  union
    {
      gsl_bst_avl_traverser avl_trav;
      gsl_bst_rb_traverser rb_trav;
    } trav_data;
} gsl_bst_trav;

/* tree types */
GSL_VAR const gsl_bst_type * gsl_bst_avl;
GSL_VAR const gsl_bst_type * gsl_bst_rb;

/*
 * Prototypes
 */

gsl_bst_workspace * gsl_bst_alloc(const gsl_bst_type * T, const gsl_bst_allocator * allocator,
                                  gsl_bst_cmp_function * compare, void * params);
void gsl_bst_free(gsl_bst_workspace * w);
int gsl_bst_empty(gsl_bst_workspace * w);
void * gsl_bst_insert(void * item, gsl_bst_workspace * w);
void * gsl_bst_find(const void * item, const gsl_bst_workspace * w);
void * gsl_bst_remove(const void * item, gsl_bst_workspace * w);
size_t gsl_bst_nodes(const gsl_bst_workspace * w);
size_t gsl_bst_node_size(const gsl_bst_workspace * w);
const char * gsl_bst_name(const gsl_bst_workspace * w);

int gsl_bst_trav_init(gsl_bst_trav * trav, const gsl_bst_workspace * w);
void * gsl_bst_trav_first(gsl_bst_trav * trav, const gsl_bst_workspace * w);
void * gsl_bst_trav_last (gsl_bst_trav * trav, const gsl_bst_workspace * w);
void * gsl_bst_trav_find (const void * item, gsl_bst_trav * trav, const gsl_bst_workspace * w);
void * gsl_bst_trav_insert (void * item, gsl_bst_trav * trav, gsl_bst_workspace * w);
void * gsl_bst_trav_copy(gsl_bst_trav * dest, const gsl_bst_trav * src);
void * gsl_bst_trav_next(gsl_bst_trav * trav);
void * gsl_bst_trav_prev(gsl_bst_trav * trav);
void * gsl_bst_trav_cur(const gsl_bst_trav * trav);
void * gsl_bst_trav_replace (gsl_bst_trav * trav, void * new_item);

__END_DECLS

#endif /* __GSL_BST_H__ */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
static char help[] =
"Coupled reaction-diffusion equations (Pearson 1993).  Option prefix -ptn_.\n"
"Demonstrates form  F(t,Y,dot Y) = G(t,Y)  where F() is IFunction and G() is\n"
"RHSFunction().  Implements IJacobian() and RHSJacobian().  Defaults to\n"
"ARKIMEX (= adaptive Runge-Kutta implicit-explicit) TS type.\n\n";

#include <petsc.h>

typedef struct {
  PetscReal u, v;
} Field;

typedef struct {
  PetscReal  L,     // domain side length
             Du,    // diffusion coefficient: u equation
             Dv,    //                        v equation
             phi,   // "dimensionless feed rate" (F in Pearson 1993)
             kappa; // "dimensionless rate constant" (k in Pearson 1993)
  PetscBool  IFcn_called, IJac_called, RHSFcn_called, RHSJac_called;
} PatternCtx;

extern PetscErrorCode InitialState(DM, Vec, PetscReal, PatternCtx*);
extern PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo*, PetscReal, Field**,
                                           Field**, PatternCtx*);
extern PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo*, PetscReal, Field**,
                                           Mat, Mat, PatternCtx*);
extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo*, PetscReal, Field**, Field**,
                                         Field **, PatternCtx*);
extern PetscErrorCode FormIJacobianLocal(DMDALocalInfo*, PetscReal, Field**, Field**,
                                         PetscReal, Mat, Mat, PatternCtx*);

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PatternCtx     user;
  TS             ts;
  Vec            x;
  DM             da;
  DMDALocalInfo  info;
  PetscReal      noiselevel = -1.0;  // negative value means no initial noise
  PetscBool      no_rhsjacobian = PETSC_FALSE,
                 no_ijacobian = PETSC_FALSE,
                 call_back_report = PETSC_FALSE;
  TSType         type;

  ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

  // parameter values from pages 21-22 in Hundsdorfer & Verwer (2003)
  user.L      = 2.5;
  user.Du     = 8.0e-5;
  user.Dv     = 4.0e-5;
  user.phi    = 0.024;
  user.kappa  = 0.06;
  user.IFcn_called   = PETSC_FALSE;
  user.IJac_called   = PETSC_FALSE;
  user.RHSFcn_called = PETSC_FALSE;
  user.RHSJac_called = PETSC_FALSE;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "ptn_", "options for patterns", ""); CHKERRQ(ierr);
  ierr = PetscOptionsBool("-call_back_report","report on which user-supplied call-backs were actually called",
           "pattern.c",call_back_report,&(call_back_report),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-Du","diffusion coefficient of first equation",
           "pattern.c",user.Du,&user.Du,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-Dv","diffusion coefficient of second equation",
           "pattern.c",user.Dv,&user.Dv,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-kappa","dimensionless rate constant (=k in (Pearson, 1993))",
           "pattern.c",user.kappa,&user.kappa,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-L","square domain side length; recommend L >= 0.5",
           "pattern.c",user.L,&user.L,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-no_ijacobian","do not set call-back DMDATSSetIJacobian()",
           "pattern.c",no_ijacobian,&(no_ijacobian),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-no_rhsjacobian","do not set call-back DMDATSSetRHSJacobian()",
           "pattern.c",no_rhsjacobian,&(no_rhsjacobian),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-noisy_init",
           "initialize u,v with this much random noise (e.g. 0.2) on top of usual initial values",
           "pattern.c",noiselevel,&noiselevel,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsReal("-phi","dimensionless feed rate (=F in (Pearson, 1993))",
           "pattern.c",user.phi,&user.phi,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
               DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC,
               DMDA_STENCIL_BOX,  // for 9-point stencil
               3,3,PETSC_DECIDE,PETSC_DECIDE,
               2, 1,              // degrees of freedom, stencil width
               NULL,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetFieldName(da,0,"u"); CHKERRQ(ierr);
  ierr = DMDASetFieldName(da,1,"v"); CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  if (info.mx != info.my) {
      SETERRQ(PETSC_COMM_SELF,1,"pattern.c requires mx == my");
  }
  ierr = DMDASetUniformCoordinates(da, 0.0, user.L, 0.0, user.L, -1.0, -1.0); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
           "running on %d x %d grid with square cells of side h = %.6f ...\n",
           info.mx,info.my,user.L/(PetscReal)(info.mx)); CHKERRQ(ierr);

//STARTTSSETUP
  ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(ierr);
  ierr = TSSetProblemType(ts,TS_NONLINEAR); CHKERRQ(ierr);
  ierr = TSSetDM(ts,da); CHKERRQ(ierr);
  ierr = TSSetApplicationContext(ts,&user); CHKERRQ(ierr);
  ierr = DMDATSSetRHSFunctionLocal(da,INSERT_VALUES,
           (DMDATSRHSFunctionLocal)FormRHSFunctionLocal,&user); CHKERRQ(ierr);
  if (!no_rhsjacobian) {
      ierr = DMDATSSetRHSJacobianLocal(da,
               (DMDATSRHSJacobianLocal)FormRHSJacobianLocal,&user); CHKERRQ(ierr);
  }
  ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,
           (DMDATSIFunctionLocal)FormIFunctionLocal,&user); CHKERRQ(ierr);
  if (!no_ijacobian) {
      ierr = DMDATSSetIJacobianLocal(da,
               (DMDATSIJacobianLocal)FormIJacobianLocal,&user); CHKERRQ(ierr);
  }
  ierr = TSSetType(ts,TSARKIMEX); CHKERRQ(ierr);
  ierr = TSSetTime(ts,0.0); CHKERRQ(ierr);
  ierr = TSSetMaxTime(ts,200.0); CHKERRQ(ierr);
  ierr = TSSetTimeStep(ts,5.0); CHKERRQ(ierr);
  ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
//ENDTSSETUP

  ierr = DMCreateGlobalVector(da,&x); CHKERRQ(ierr);
  ierr = InitialState(da,x,noiselevel,&user); CHKERRQ(ierr);
  ierr = TSSolve(ts,x); CHKERRQ(ierr);

  // optionally report on call-backs
  if (call_back_report) {
      ierr = TSGetType(ts,&type);CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,"CALL-BACK REPORT\n  solver type: %s\n",type); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,"  IFunction:   %D  | IJacobian:   %D\n",
                                          (int)user.IFcn_called,(int)user.IJac_called); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,"  RHSFunction: %D  | RHSJacobian: %D\n",
                                          (int)user.RHSFcn_called,(int)user.RHSJac_called); CHKERRQ(ierr);
  }

  VecDestroy(&x);  TSDestroy(&ts);  DMDestroy(&da);
  return PetscFinalize();
}

// Formulas from page 22 of Hundsdorfer & Verwer (2003).  Interpretation here is
// to always generate 0.5 x 0.5 non-trivial patch in (0,L) x (0,L) domain.
PetscErrorCode InitialState(DM da, Vec Y, PetscReal noiselevel, PatternCtx* user) {
  PetscErrorCode ierr;
  DMDALocalInfo    info;
  PetscInt         i,j;
  PetscReal        sx,sy;
  const PetscReal  ledge = (user->L - 0.5) / 2.0, // nontrivial initial values on
                   redge = user->L - ledge;       //   ledge < x,y < redge
  DMDACoor2d       **aC;
  Field            **aY;

  ierr = VecSet(Y,0.0); CHKERRQ(ierr);
  if (noiselevel > 0.0) {
      // noise added to usual initial condition is uniform on [0,noiselevel],
      //     independently for each location and component
      ierr = VecSetRandom(Y,NULL); CHKERRQ(ierr);
      ierr = VecScale(Y,noiselevel); CHKERRQ(ierr);
  }
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  ierr = DMDAGetCoordinateArray(da,&aC); CHKERRQ(ierr);
  ierr = DMDAVecGetArray(da,Y,&aY); CHKERRQ(ierr);
  for (j = info.ys; j < info.ys+info.ym; j++) {
    for (i = info.xs; i < info.xs+info.xm; i++) {
      if ((aC[j][i].x >= ledge) && (aC[j][i].x <= redge)
              && (aC[j][i].y >= ledge) && (aC[j][i].y <= redge)) {
          sx = PetscSinReal(4.0 * PETSC_PI * aC[j][i].x);
          sy = PetscSinReal(4.0 * PETSC_PI * aC[j][i].y);
          aY[j][i].v += 0.5 * sx * sx * sy * sy;
      }
      aY[j][i].u += 1.0 - 2.0 * aY[j][i].v;
    }
  }
  ierr = DMDAVecRestoreArray(da,Y,&aY); CHKERRQ(ierr);
  ierr = DMDARestoreCoordinateArray(da,&aC); CHKERRQ(ierr);
  return 0;
}

// in system form  F(t,Y,dot Y) = G(t,Y),  compute G():
//     G^u(t,u,v) = - u v^2 + phi (1 - u)
//     G^v(t,u,v) = + u v^2 - (phi + kappa) v
//STARTRHSFUNCTION
PetscErrorCode FormRHSFunctionLocal(DMDALocalInfo *info,
                   PetscReal t, Field **aY, Field **aG, PatternCtx *user) {
  PetscInt   i, j;
  PetscReal  uv2;

  user->RHSFcn_called = PETSC_TRUE;
  for (j = info->ys; j < info->ys + info->ym; j++) {
      for (i = info->xs; i < info->xs + info->xm; i++) {
          uv2 = aY[j][i].u * aY[j][i].v * aY[j][i].v;
          aG[j][i].u = - uv2 + user->phi * (1.0 - aY[j][i].u);
          aG[j][i].v = + uv2 - (user->phi + user->kappa) * aY[j][i].v;
      }
  }
  return 0;
}
//ENDRHSFUNCTION

PetscErrorCode FormRHSJacobianLocal(DMDALocalInfo *info,
                                    PetscReal t, Field **aY,
                                    Mat J, Mat P, PatternCtx *user) {
    PetscErrorCode ierr;
    PetscInt    i, j;
    PetscReal   v[2], uv, v2;
    MatStencil  col[2],row;

    user->RHSJac_called = PETSC_TRUE;
    for (j = info->ys; j < info->ys+info->ym; j++) {
        row.j = j;  col[0].j = j;  col[1].j = j;
        for (i = info->xs; i < info->xs+info->xm; i++) {
            row.i = i;  col[0].i = i;  col[1].i = i;
            uv = aY[j][i].u * aY[j][i].v;
            v2 = aY[j][i].v * aY[j][i].v;
            // u equation
            row.c = 0;  col[0].c = 0;  col[1].c = 1;
            v[0] = - v2 - user->phi;
            v[1] = - 2.0 * uv;
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
            // v equation
            row.c = 1;  col[0].c = 0;  col[1].c = 1;
            v[0] = v2;
            v[1] = 2.0 * uv - (user->phi + user->kappa);
            ierr = MatSetValuesStencil(P,1,&row,2,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}

// in system form  F(t,Y,dot Y) = G(t,Y),  compute F():
//     F^u(t,u,v,u_t,v_t) = u_t - D_u Laplacian u
//     F^v(t,u,v,u_t,v_t) = v_t - D_v Laplacian v
//STARTIFUNCTION
PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscReal t,
                                  Field **aY, Field **aYdot, Field **aF,
                                  PatternCtx *user) {
  PetscInt         i, j;
  const PetscReal  h = user->L / (PetscReal)(info->mx),
                   Cu = user->Du / (6.0 * h * h),
                   Cv = user->Dv / (6.0 * h * h);
  PetscReal        u, v, lapu, lapv;

  user->IFcn_called = PETSC_TRUE;
  for (j = info->ys; j < info->ys + info->ym; j++) {
      for (i = info->xs; i < info->xs + info->xm; i++) {
          u = aY[j][i].u;
          v = aY[j][i].v;
          lapu =     aY[j+1][i-1].u + 4.0*aY[j+1][i].u +   aY[j+1][i+1].u
                 + 4.0*aY[j][i-1].u -    20.0*u        + 4.0*aY[j][i+1].u
                 +   aY[j-1][i-1].u + 4.0*aY[j-1][i].u +   aY[j-1][i+1].u;
          lapv =     aY[j+1][i-1].v + 4.0*aY[j+1][i].v +   aY[j+1][i+1].v
                 + 4.0*aY[j][i-1].v -    20.0*v        + 4.0*aY[j][i+1].v
                 +   aY[j-1][i-1].v + 4.0*aY[j-1][i].v +   aY[j-1][i+1].v;
          aF[j][i].u = aYdot[j][i].u - Cu * lapu;
          aF[j][i].v = aYdot[j][i].v - Cv * lapv;
      }
  }
  return 0;
}
//ENDIFUNCTION

// in system form  F(t,Y,dot Y) = G(t,Y),  compute combined/shifted
// Jacobian of F():
//     J = (shift) dF/d(dot Y) + dF/dY
//STARTIJACOBIAN
PetscErrorCode FormIJacobianLocal(DMDALocalInfo *info,
                   PetscReal t, Field **aY, Field **aYdot,
                   PetscReal shift, Mat J, Mat P,
                   PatternCtx *user) {
    PetscErrorCode ierr;
    PetscInt         i, j, s, c;
    const PetscReal  h = user->L / (PetscReal)(info->mx),
                     Cu = user->Du / (6.0 * h * h),
                     Cv = user->Dv / (6.0 * h * h);
    PetscReal        val[9], CC;
    MatStencil       col[9], row;

    ierr = MatZeroEntries(P); CHKERRQ(ierr);  // workaround to address PETSc issue #734
    user->IJac_called = PETSC_TRUE;
    for (j = info->ys; j < info->ys + info->ym; j++) {
        row.j = j;
        for (i = info->xs; i < info->xs + info->xm; i++) {
            row.i = i;
            for (c = 0; c < 2; c++) { // u,v equations are c=0,1
                row.c = c;
                CC = (c == 0) ? Cu : Cv;
                for (s = 0; s < 9; s++)
                    col[s].c = c;
                col[0].i = i;   col[0].j = j;
                val[0] = shift + 20.0 * CC;
                col[1].i = i-1; col[1].j = j;    val[1] = - 4.0 * CC;
                col[2].i = i+1; col[2].j = j;    val[2] = - 4.0 * CC;
                col[3].i = i;   col[3].j = j-1;  val[3] = - 4.0 * CC;
                col[4].i = i;   col[4].j = j+1;  val[4] = - 4.0 * CC;
                col[5].i = i-1; col[5].j = j-1;  val[5] = - CC;
                col[6].i = i-1; col[6].j = j+1;  val[6] = - CC;
                col[7].i = i+1; col[7].j = j-1;  val[7] = - CC;
                col[8].i = i+1; col[8].j = j+1;  val[8] = - CC;
                ierr = MatSetValuesStencil(P,1,&row,9,col,val,INSERT_VALUES); CHKERRQ(ierr);
            }
        }
    }

    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}
//ENDIJACOBIAN

ext: c


text:
/* ode-initval/rk4imp.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Runge-Kutta 4, Gaussian implicit */

/* Author:  G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k1nu;
  double *k2nu;
  double *ytmp1;
  double *ytmp2;
}
rk4imp_state_t;

static void *
rk4imp_alloc (size_t dim)
{
  rk4imp_state_t *state = (rk4imp_state_t *) malloc (sizeof (rk4imp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk4imp_state",
                      GSL_ENOMEM);
    }

  state->k1nu = (double *) malloc (dim * sizeof (double));

  if (state->k1nu == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k1nu", GSL_ENOMEM);
    }

  state->k2nu = (double *) malloc (dim * sizeof (double));

  if (state->k2nu == 0)
    {
      free (state->k1nu);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k2nu", GSL_ENOMEM);
    }

  state->ytmp1 = (double *) malloc (dim * sizeof (double));

  if (state->ytmp1 == 0)
    {
      free (state->k2nu);
      free (state->k1nu);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp1", GSL_ENOMEM);
    }

  state->ytmp2 = (double *) malloc (dim * sizeof (double));

  if (state->ytmp2 == 0)
    {
      free (state->ytmp1);
      free (state->k2nu);
      free (state->k1nu);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp2", GSL_ENOMEM);
    }

  return state;
}

static int
rk4imp_apply (void *vstate,
              size_t dim,
              double t,
              double h,
              double y[],
              double yerr[],
              const double dydt_in[],
              double dydt_out[],
              const gsl_odeiv_system * sys)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  const double ir3 = 1.0 / M_SQRT3;
  const int iter_steps = 3;
  int status = 0;
  int nu;
  size_t i;

  double *const k1nu = state->k1nu;
  double *const k2nu = state->k2nu;
  double *const ytmp1 = state->ytmp1;
  double *const ytmp2 = state->ytmp2;

  /* initialization step */
  if (dydt_in != 0)
    {
      DBL_MEMCPY (k1nu, dydt_in, dim);
    }
  else
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t, y, k1nu);
      GSL_STATUS_UPDATE (&status, s);
    }

  DBL_MEMCPY (k2nu, k1nu, dim);

  /* iterative solution */
  for (nu = 0; nu < iter_steps; nu++)
    {
      for (i = 0; i < dim; i++)
        {
          ytmp1[i] =
            y[i] + h * (0.25 * k1nu[i] + 0.5 * (0.5 - ir3) * k2nu[i]);
          ytmp2[i] =
            y[i] + h * (0.25 * k2nu[i] + 0.5 * (0.5 + ir3) * k1nu[i]);
        }
      {
        int s =
          GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h * (1.0 - ir3), ytmp1, k1nu);
        GSL_STATUS_UPDATE (&status, s);
      }
      {
        int s =
          GSL_ODEIV_FN_EVAL (sys, t + 0.5 * h * (1.0 + ir3), ytmp2, k2nu);
        GSL_STATUS_UPDATE (&status, s);
      }
    }

  /* assignment */
  for (i = 0; i < dim; i++)
    {
      const double d_i = 0.5 * (k1nu[i] + k2nu[i]);
      if (dydt_out != NULL)
        dydt_out[i] = d_i;
      y[i] += h * d_i;
      yerr[i] = h * h * d_i;    /* FIXME: is this an overestimate ? */
    }

  return status;
}

static int
rk4imp_reset (void *vstate, size_t dim)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->k1nu, dim);
  DBL_ZERO_MEMSET (state->k2nu, dim);
  DBL_ZERO_MEMSET (state->ytmp1, dim);
  DBL_ZERO_MEMSET (state->ytmp2, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk4imp_order (void *vstate)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 4;
}

static void
rk4imp_free (void *vstate)
{
  rk4imp_state_t *state = (rk4imp_state_t *) vstate;
  free (state->k1nu);
  free (state->k2nu);
  free (state->ytmp1);
  free (state->ytmp2);
  free (state);
}

static const gsl_odeiv_step_type rk4imp_type = { "rk4imp",      /* name */
  1,                            /* can use dydt_in */
  0,                            /* gives exact dydt_out */
  &rk4imp_alloc,
  &rk4imp_apply,
  &rk4imp_reset,
  &rk4imp_order,
  &rk4imp_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_rk4imp = &rk4imp_type;


text:
/*
 * This file is part of the GROMACS molecular simulation package.
 *
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
 * Copyright (c) 2001-2004, The GROMACS development team,
 * check out http://www.gromacs.org for more information.
 * Copyright (c) 2012,2013, by the GROMACS development team, led by
 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
 * others, as listed in the AUTHORS file in the top-level source
 * directory and at http://www.gromacs.org.
 *
 * GROMACS is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 *
 * GROMACS 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
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with GROMACS; if not, see
 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 *
 * If you want to redistribute modifications to GROMACS, please
 * consider that scientific software is very special. Version
 * control is crucial - bugs must be traceable. We will be happy to
 * consider code for inclusion in the official distribution, but
 * derived work must not be called official GROMACS. Details are found
 * in the README & COPYING files - if they are missing, get the
 * official version at http://www.gromacs.org.
 *
 * To help us fund GROMACS development, we humbly ask that you cite
 * the research papers on the package. Check out http://www.gromacs.org.
 */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#ifdef HAVE_LIBGSL
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_version.h>
#endif

#include <math.h>
#include "typedefs.h"
#include "smalloc.h"
#include "vec.h"
#include "geminate.h"
#include "gmx_omp.h"

/* The first few sections of this file contain functions that were adopted,
 * and to some extent modified, by Erik Marklund (erikm[aT]xray.bmc.uu.se,
 * http://folding.bmc.uu.se) from code written by Omer Markovitch (email, url).
 * This is also the case with the function eq10v2().
 *
 * The parts menetioned in the previous paragraph were contributed under the BSD license.
 */

/* This first part is from complex.c which I recieved from Omer Markowitch.
 * - Erik Marklund
 *
 * ------------- from complex.c ------------- */

/* Complexation of a paired number (r,i)                                     */
static gem_complex gem_cmplx(double r, double i)
{
    gem_complex value;
    value.r = r;
    value.i = i;
    return value;
}

/* Complexation of a real number, x */
static gem_complex gem_c(double x)
{
    gem_complex value;
    value.r = x;
    value.i = 0;
    return value;
}

/* Magnitude of a complex number z                                           */
static double gem_cx_abs(gem_complex z)
{
    return (sqrt(z.r*z.r+z.i*z.i));
}

/* Addition of two complex numbers z1 and z2                                 */
static gem_complex gem_cxadd(gem_complex z1, gem_complex z2)
{
    gem_complex value;
    value.r = z1.r+z2.r;
    value.i = z1.i+z2.i;
    return value;
}

/* Addition of a complex number z1 and a real number r */
static gem_complex gem_cxradd(gem_complex z, double r)
{
    gem_complex value;
    value.r = z.r + r;
    value.i = z.i;
    return value;
}

/* Subtraction of two complex numbers z1 and z2                              */
static gem_complex gem_cxsub(gem_complex z1, gem_complex z2)
{
    gem_complex value;
    value.r = z1.r-z2.r;
    value.i = z1.i-z2.i;
    return value;
}

/* Multiplication of two complex numbers z1 and z2                           */
static gem_complex gem_cxmul(gem_complex z1, gem_complex z2)
{
    gem_complex value;
    value.r = z1.r*z2.r-z1.i*z2.i;
    value.i = z1.r*z2.i+z1.i*z2.r;
    return value;
}

/* Square of a complex number z                                              */
static gem_complex gem_cxsq(gem_complex z)
{
    gem_complex value;
    value.r = z.r*z.r-z.i*z.i;
    value.i = z.r*z.i*2.;
    return value;
}

/* multiplication of a complex number z and a real number r */
static gem_complex gem_cxrmul(gem_complex z, double r)
{
    gem_complex value;
    value.r = z.r*r;
    value.i = z.i*r;
    return value;
}

/* Division of two complex numbers z1 and z2                                 */
static gem_complex gem_cxdiv(gem_complex z1, gem_complex z2)
{
    gem_complex value;
    double      num;
    num = z2.r*z2.r+z2.i*z2.i;
    if (num == 0.)
    {
        fprintf(stderr, "ERROR in gem_cxdiv function\n");
    }
    value.r = (z1.r*z2.r+z1.i*z2.i)/num; value.i = (z1.i*z2.r-z1.r*z2.i)/num;
    return value;
}

/* division of a complex z number by a real number x */
static gem_complex gem_cxrdiv(gem_complex z, double r)
{
    gem_complex value;
    value.r = z.r/r;
    value.i = z.i/r;
    return value;
}

/* division of a real number r by a complex number x */
static gem_complex gem_rxcdiv(double r, gem_complex z)
{
    gem_complex value;
    double      f;
    f       = r/(z.r*z.r+z.i*z.i);
    value.r = f*z.r;
    value.i = -f*z.i;
    return value;
}

/* Exponential of a complex number-- exp (z)=|exp(z.r)|*{cos(z.i)+I*sin(z.i)}*/
static gem_complex gem_cxdexp(gem_complex z)
{
    gem_complex value;
    double      exp_z_r;
    exp_z_r = exp(z.r);
    value.r = exp_z_r*cos(z.i);
    value.i = exp_z_r*sin(z.i);
    return value;
}

/* Logarithm of a complex number -- log(z)=log|z|+I*Arg(z),                  */
/*                                  where -PI < Arg(z) < PI                  */
static gem_complex gem_cxlog(gem_complex z)
{
    gem_complex value;
    double      mag2;
    mag2 = z.r*z.r+z.i*z.i;
    if (mag2 < 0.)
    {
        fprintf(stderr, "ERROR in gem_cxlog func\n");
    }
    value.r = log(sqrt(mag2));
    if (z.r == 0.)
    {
        value.i = PI/2.;
        if (z.i < 0.)
        {
            value.i = -value.i;
        }
    }
    else
    {
        value.i = atan2(z.i, z.r);
    }
    return value;
}

/* Square root of a complex number z = |z| exp(I*the) -- z^(1/2)             */
/*                               z^(1/2)=|z|^(1/2)*[cos(the/2)+I*sin(the/2)] */
/*                               where 0 < the < 2*PI                        */
static gem_complex gem_cxdsqrt(gem_complex z)
{
    gem_complex value;
    double      sq;
    sq      = gem_cx_abs(z);
    value.r = sqrt(fabs((sq+z.r)*0.5)); /* z'.r={|z|*[1+cos(the)]/2}^(1/2) */
    value.i = sqrt(fabs((sq-z.r)*0.5)); /* z'.i={|z|*[1-cos(the)]/2}^(1/2) */
    if (z.i < 0.)
    {
        value.r = -value.r;
    }
    return value;
}

/* Complex power of a complex number z1^z2                                   */
static gem_complex gem_cxdpow(gem_complex z1, gem_complex z2)
{
    gem_complex value;
    value = gem_cxdexp(gem_cxmul(gem_cxlog(z1), z2));
    return value;
}

/* ------------ end of complex.c ------------ */

/* This next part was derived from cubic.c, also received from Omer Markovitch.
 * ------------- from cubic.c ------------- */

/* Solver for a cubic equation: x^3-a*x^2+b*x-c=0                            */
static void gem_solve(gem_complex *al, gem_complex *be, gem_complex *gam,
                      double a, double b, double c)
{
    double      t1, t2, two_3, temp;
    gem_complex ctemp, ct3;

    two_3 = pow(2., 1./3.); t1 = -a*a+3.*b; t2 = 2.*a*a*a-9.*a*b+27.*c;
    temp  = 4.*t1*t1*t1+t2*t2;

    ctemp = gem_cmplx(temp, 0.);   ctemp = gem_cxadd(gem_cmplx(t2, 0.), gem_cxdsqrt(ctemp));
    ct3   = gem_cxdpow(ctemp, gem_cmplx(1./3., 0.));

    ctemp = gem_rxcdiv(-two_3*t1/3., ct3);
    ctemp = gem_cxadd(ctemp, gem_cxrdiv(ct3, 3.*two_3));

    *gam = gem_cxadd(gem_cmplx(a/3., 0.), ctemp);

    ctemp = gem_cxmul(gem_cxsq(*gam), gem_cxsq(gem_cxsub(*gam, gem_cmplx(a, 0.))));
    ctemp = gem_cxadd(ctemp, gem_cxmul(gem_cmplx(-4.*c, 0.), *gam));
    ctemp = gem_cxdiv(gem_cxdsqrt(ctemp), *gam);
    *al   = gem_cxrmul(gem_cxsub(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
    *be   = gem_cxrmul(gem_cxadd(gem_cxsub(gem_cmplx(a, 0.), *gam), ctemp), 0.5);
}

/* ------------ end of cubic.c ------------ */

/* This next part was derived from cerror.c and rerror.c, also received from Omer Markovitch.
 * ------------- from [cr]error.c ------------- */

/************************************************************/
/* Real valued error function and related functions         */
/* x, y     : real variables                                */
/* erf(x)   : error function                                */
/* erfc(x)  : complementary error function                  */
/* omega(x) : exp(x*x)*erfc(x)                              */
/* W(x,y)   : exp(-x*x)*omega(x+y)=exp(2*x*y+y^2)*erfc(x+y) */
/************************************************************/

/*---------------------------------------------------------------------------*/
/* Utilzed the series approximation of erf(x)                                */
/* Relative error=|err(x)|/erf(x)<EPS                                        */
/* Handbook of Mathematical functions, Abramowitz, p 297                     */
/* Note: When x>=6 series sum deteriorates -> Asymptotic series used instead */
/*---------------------------------------------------------------------------*/
/*  This gives erfc(z) correctly only upto >10-15 */

static double gem_erf(double x)
{
    double n, sum, temp, exp2, xm, x2, x4, x6, x8, x10, x12;
    temp = x;
    sum  = temp;
    xm   = 26.;
    x2   = x*x;
    x4   = x2*x2;
    x6   = x4*x2;
    x8   = x6*x2;
    x10  = x8*x2;
    x12  = x10*x2;
    exp2 = exp(-x2);
    if (x <= xm)
    {
        for (n = 1.; n <= 2000.; n += 1.)
        {
            temp *= 2.*x2/(2.*n+1.);
            sum  += temp;
            if (fabs(temp/sum) < 1.E-16)
            {
                break;
            }
        }

        if (n >= 2000.)
        {
            fprintf(stderr, "In Erf calc - iteration exceeds %lg\n", n);
        }
        sum *= 2./sPI*exp2;
    }
    else
    {
        /* from the asymptotic expansion of experfc(x) */
        sum = (1. - 0.5/x2 + 0.75/x4
               - 1.875/x6 + 6.5625/x8
               - 29.53125/x10 + 162.421875/x12)
            / sPI/x;
        sum *= exp2; /* now sum is erfc(x) */
        sum  = -sum+1.;
    }
    return x >= 0.0 ? sum : -sum;
}

/* Result --> Alex's code for erfc and experfc behaves better than this      */
/* complementray error function.                Returns 1.-erf(x)            */
static double gem_erfc(double x)
{
    double t, z, ans;
    z = fabs(x);
    t = 1.0/(1.0+0.5*z);

    ans = t * exp(-z*z - 1.26551223 +
                  t*(1.00002368 +
                     t*(0.37409196 +
                        t*(0.09678418 +
                           t*(-0.18628806 +
                              t*(0.27886807 +
                                 t*(-1.13520398 +
                                    t*(1.48851587 +
                                       t*(-0.82215223 +
                                          t*0.17087277)))))))));

    return x >= 0.0 ? ans : 2.0-ans;
}

/* omega(x)=exp(x^2)erfc(x)                                                  */
static double gem_omega(double x)
{
    double xm, ans, xx, x4, x6, x8, x10, x12;
    xm  = 26;
    xx  = x*x;
    x4  = xx*xx;
    x6  = x4*xx;
    x8  = x6*xx;
    x10 = x8*xx;
    x12 = x10*xx;

    if (x <= xm)
    {
        ans = exp(xx)*gem_erfc(x);
    }
    else
    {
        /* Asymptotic expansion */
        ans = (1. - 0.5/xx + 0.75/x4 - 1.875/x6 + 6.5625/x8 - 29.53125/x10 + 162.421875/x12) / sPI/x;
    }
    return ans;
}

/*---------------------------------------------------------------------------*/
/* Utilzed the series approximation of erf(z=x+iy)                           */
/* Relative error=|err(z)|/|erf(z)|<EPS                                      */
/* Handbook of Mathematical functions, Abramowitz, p 299                     */
/* comega(z=x+iy)=cexp(z^2)*cerfc(z)                                         */
/*---------------------------------------------------------------------------*/
static gem_complex gem_comega(gem_complex z)
{
    gem_complex value;
    double      x, y;
    double      sumr, sumi, n, n2, f, temp, temp1;
    double      x2, y2, cos_2xy, sin_2xy, cosh_2xy, sinh_2xy, cosh_ny, sinh_ny, exp_y2;

    x        = z.r;
    y        = z.i;
    x2       = x*x;
    y2       = y*y;
    sumr     = 0.;
    sumi     = 0.;
    cos_2xy  = cos(2.*x*y);
    sin_2xy  = sin(2.*x*y);
    cosh_2xy = cosh(2.*x*y);
    sinh_2xy = sinh(2.*x*y);
    exp_y2   = exp(-y2);

    for (n = 1.0, temp = 0.; n <= 2000.; n += 1.0)
    {
        n2      = n*n;
        cosh_ny = cosh(n*y);
        sinh_ny = sinh(n*y);
        f       = exp(-n2/4.)/(n2+4.*x2);
        /* if(f<1.E-200) break; */
        sumr += (2.*x - 2.*x*cosh_ny*cos_2xy + n*sinh_ny*sin_2xy)*f;
        sumi += (2.*x*cosh_ny*sin_2xy + n*sinh_ny*cos_2xy)*f;
        temp1 = sqrt(sumr*sumr+sumi*sumi);
        if (fabs((temp1-temp)/temp1) < 1.E-16)
        {
            break;
        }
        temp = temp1;
    }
    if (n == 2000.)
    {
        fprintf(stderr, "iteration exceeds %lg\n", n);
    }
    sumr *= 2./PI;
    sumi *= 2./PI;

    if (x != 0.)
    {
        f = 1./2./PI/x;
    }
    else
    {
        f = 0.;
    }
    value.r = gem_omega(x)-(f*(1.-cos_2xy)+sumr);
    value.i = -(f*sin_2xy+sumi);
    value   = gem_cxmul(value, gem_cmplx(exp_y2*cos_2xy, exp_y2*sin_2xy));
    return (value);
}

/* ------------ end of [cr]error.c ------------ */

/*_ REVERSIBLE GEMINATE RECOMBINATION
 *
 * Here are the functions for reversible geminate recombination. */

/* Changes the unit from square cm per s to square Ångström per ps,
 * since Omers code uses the latter units while g_mds outputs the former.
 * g_hbond expects a diffusion coefficent given in square cm per s. */
static double sqcm_per_s_to_sqA_per_ps (real D)
{
    fprintf(stdout, "Diffusion coefficient is %f A^2/ps\n", D*1e4);
    return (double)(D*1e4);
}

static double eq10v2(double theoryCt[], double time[], int manytimes,
                     double ka, double kd, t_gemParams *params)
{
    /* Finding the 3 roots */
    int    i;
    double kD, D, r, a, b, c, tsqrt, sumimaginary;
    gem_complex
           alpha, beta, gamma,
        c1, c2, c3, c4,
        oma, omb, omc,
        part1, part2, part3, part4;

    kD = params->kD;
    D  = params->D;
    r  = params->sigma;
    a  = (1.0 + ka/kD) * sqrt(D)/r;
    b  = kd;
    c  = kd * sqrt(D)/r;

    gem_solve(&alpha, &beta, &gamma, a, b, c);
    /* Finding the 3 roots */

    sumimaginary = 0;
    part1        = gem_cxmul(alpha, gem_cxmul(gem_cxadd(beta,  gamma), gem_cxsub(beta,  gamma)));                 /* 1(2+3)(2-3) */
    part2        = gem_cxmul(beta,  gem_cxmul(gem_cxadd(gamma, alpha), gem_cxsub(gamma, alpha)));                 /* 2(3+1)(3-1) */
    part3        = gem_cxmul(gamma, gem_cxmul(gem_cxadd(alpha, beta), gem_cxsub(alpha, beta)));                   /* 3(1+2)(1-2) */
    part4        = gem_cxmul(gem_cxsub(gamma, alpha), gem_cxmul(gem_cxsub(alpha, beta), gem_cxsub(beta, gamma))); /* (3-1)(1-2)(2-3) */

#pragma omp parallel for \
    private(i, tsqrt, oma, omb, omc, c1, c2, c3, c4) \
    reduction(+:sumimaginary) \
    default(shared) \
    schedule(guided)
    for (i = 0; i < manytimes; i++)
    {
        tsqrt         = sqrt(time[i]);
        oma           = gem_comega(gem_cxrmul(alpha, tsqrt));
        c1            = gem_cxmul(oma, gem_cxdiv(part1, part4));
        omb           = gem_comega(gem_cxrmul(beta, tsqrt));
        c2            = gem_cxmul(omb, gem_cxdiv(part2, part4));
        omc           = gem_comega(gem_cxrmul(gamma, tsqrt));
        c3            = gem_cxmul(omc, gem_cxdiv(part3, part4));
        c4.r          = c1.r+c2.r+c3.r;
        c4.i          = c1.i+c2.i+c3.i;
        theoryCt[i]   = c4.r;
        sumimaginary += c4.i * c4.i;
    }

    return sumimaginary;

} /* eq10v2 */

/* This returns the real-valued index(!) to an ACF, equidistant on a log scale. */
static double getLogIndex(const int i, const t_gemParams *params)
{
    return (exp(((double)(i)) * params->logQuota) -1);
}

extern t_gemParams *init_gemParams(const double sigma, const double D,
                                   const real *t, const int len, const int nFitPoints,
                                   const real begFit, const real endFit,
                                   const real ballistic, const int nBalExp, const gmx_bool bDt)
{
    double       tDelta;
    t_gemParams *p;

    snew(p, 1);

    /* A few hardcoded things here. For now anyway. */
/*   p->ka_min   = 0; */
/*   p->ka_max   = 100; */
/*   p->dka      = 10; */
/*   p->kd_min   = 0; */
/*   p->kd_max   = 2; */
/*   p->dkd      = 0.2; */
    p->ka       = 0;
    p->kd       = 0;
/*   p->lsq      = -1; */
/*   p->lifetime = 0; */
    p->sigma    = sigma*10; /* Omer uses Å, not nm */
/*   p->lsq_old  = 99999; */
    p->D        = sqcm_per_s_to_sqA_per_ps(D);
    p->kD       = 4*acos(-1.0)*sigma*p->D;

    /* Parameters used by calcsquare(). Better to calculate them
     * here than in calcsquare every time it's called. */
    p->len = len;
/*   p->logAfterTime = logAfterTime; */
    tDelta       = (t[len-1]-t[0]) / len;
    if (tDelta <= 0)
    {
        gmx_fatal(FARGS, "Time between frames is non-positive!");
    }
    else
    {
        p->tDelta = tDelta;
    }

    p->nExpFit      = nBalExp;
/*   p->nLin         = logAfterTime / tDelta; */
    p->nFitPoints   = nFitPoints;
    p->begFit       = begFit;
    p->endFit       = endFit;
    p->logQuota     = (double)(log(p->len))/(p->nFitPoints-1);
/*   if (p->nLin <= 0) { */
/*     fprintf(stderr, "Number of data points in the linear regime is non-positive!\n"); */
/*     sfree(p); */
/*     return NULL; */
/*   } */
/* We want the same number of data points in the log regime. Not crucial, but seems like a good idea. */
/* p->logDelta = log(((float)len)/p->nFitPoints) / p->nFitPoints;/\* log(((float)len)/p->nLin) / p->nLin; *\/ */
/*   p->logPF    = p->nFitPoints*p->nFitPoints/(float)len; /\* p->nLin*p->nLin/(float)len; *\/ */
/* logPF and logDelta are stitched together with the macro GETLOGINDEX defined in geminate.h */

/*   p->logMult      = pow((float)len, 1.0/nLin);/\* pow(t[len-1]-t[0], 1.0/p->nLin); *\/ */
    p->ballistic    =  ballistic;
    return p;
}

/* There was a misunderstanding regarding the fitting. From our
 * recent correspondence it appears that Omer's code require
 * the ACF data on a log-scale and does not operate on the raw data.
 * This needs to be redone in gemFunc_residual() as well as in the
 * t_gemParams structure. */
#ifdef HAVE_LIBGSL
static double gemFunc_residual2(const gsl_vector *p, void *data)
{
    gemFitData *GD;
    int         i, iLog, nLin, nFitPoints, nData;
    double      r, residual2, *ctTheory, *y;

    GD         = (gemFitData *)data;
    nLin       = GD->params->nLin;
    nFitPoints = GD->params->nFitPoints;
    nData      = GD->nData;
    residual2  = 0;
    ctTheory   = GD->ctTheory;
    y          = GD->y;

    /* Now, we'd save a lot of time by not calculating eq10v2 for all
     * time points, but only those that we sample to calculate the mse.
     */

    eq10v2(GD->ctTheory, GD->doubleLogTime /* GD->time */, nFitPoints /* GD->nData */,
           gsl_vector_get(p, 0), gsl_vector_get(p, 1),
           GD->params);

    fixGemACF(GD->ctTheory, nFitPoints);

    /* Removing a bunch of points from the log-part. */
#pragma omp parallel for schedule(dynamic) \
    firstprivate(nData, ctTheory, y, nFitPoints) \
    private (i, iLog, r) \
    reduction(+:residual2) \
    default(shared)
    for (i = 0; i < nFitPoints; i++)
    {
        iLog       = GD->logtime[i];
        r          = log(ctTheory[i /* iLog */]);
        residual2 += sqr(r-log(y[iLog]));
    }
    residual2 /= nFitPoints; /* Not really necessary for the fitting, but gives more meaning to the returned data. */
    /* printf("ka=%3.5f\tkd=%3.5f\trmse=%3.5f\n", gsl_vector_get(p,0), gsl_vector_get(p,1), residual2); */
    return residual2;

/*   for (i=0; i<nLin; i++) { */
/*     /\* Linear part ----------*\/ */
/*     r = ctTheory[i]; */
/*     residual2 += sqr(r-y[i]); */
/*     /\* Log part -------------*\/ */
/*     iLog = GETLOGINDEX(i, GD->params); */
/*     if (iLog >= nData) */
/*       gmx_fatal(FARGS, "log index out of bounds: %i", iLog); */
/*     r = ctTheory[iLog]; */
/*     residual2 += sqr(r-y[iLog]); */

/*   } */
/*   residual2 /= GD->n; */
/*   return residual2; */
}
#endif

static real* d2r(const double *d, const int nn);

extern real fitGemRecomb(double *ct, double *time, double **ctFit,
                         const int nData, t_gemParams *params)
{

    int         nThreads, i, iter, status, maxiter;
    real        size, d2, tol, *dumpdata;
    size_t      p, n;
    gemFitData *GD;
    char       *dumpstr, dumpname[128];

    /* nmsimplex2 had convergence problems prior to gsl v1.14,
     * but it's O(N) instead of O(N) operations, so let's use it if v >= 1.14 */
#ifdef HAVE_LIBGSL
    gsl_multimin_fminimizer *s;
    gsl_vector              *x, *dx; /* parameters and initial step size */
    gsl_multimin_function    fitFunc;
#ifdef GSL_MAJOR_VERSION
#ifdef GSL_MINOR_VERSION
#if ((GSL_MAJOR_VERSION == 1 && GSL_MINOR_VERSION >= 14) || \
    (GSL_MAJOR_VERSION > 1))
    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2;
#else
    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
#endif /* #if ... */
#endif /* GSL_MINOR_VERSION */
#else
    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
#endif /* GSL_MAJOR_VERSION */
    fprintf(stdout, "Will fit ka and kd to the ACF according to the reversible geminate recombination model.\n");
#else  /* HAVE_LIBGSL */
    fprintf(stderr, "Sorry, can't do reversible geminate recombination without gsl. "
            "Recompile using --with-gsl.\n");
    return -1;
#endif /* HAVE_LIBGSL */

#ifdef HAVE_LIBGSL
#ifdef GMX_OPENMP
    nThreads = gmx_omp_get_max_threads();
    gmx_omp_set_num_threads(nThreads);
    fprintf(stdout, "We will be using %i threads.\n", nThreads);
#endif

    iter    = 0;
    status  = 0;
    maxiter = 100;
    tol     = 1e-10;

    p = 2;                                        /* Number of parameters to fit. ka and kd.  */
    n = params->nFitPoints; /* params->nLin*2 */; /* Number of points in the reduced dataset  */

    if (params->D <= 0)
    {
        fprintf(stderr, "Fitting of D is not implemented yet. It must be provided on the command line.\n");
        return -1;
    }

/*   if (nData<n) { */
/*     fprintf(stderr, "Reduced data set larger than the complete data set!\n"); */
/*     n=nData; */
/*   } */
    snew(dumpdata, nData);
    snew(GD, 1);

    GD->n        = n;
    GD->y        = ct;
    GD->ctTheory = NULL;
    snew(GD->ctTheory, nData);
    GD->LinLog = NULL;
    snew(GD->LinLog, n);
    GD->time   = time;
    GD->ka     = 0;
    GD->kd     = 0;
    GD->tDelta = time[1]-time[0];
    GD->nData  = nData;
    GD->params = params;
    snew(GD->logtime, params->nFitPoints);
    snew(GD->doubleLogTime, params->nFitPoints);

    for (i = 0; i < params->nFitPoints; i++)
    {
        GD->doubleLogTime[i]  = (double)(getLogIndex(i, params));
        GD->logtime[i]        = (int)(GD->doubleLogTime[i]);
        GD->doubleLogTime[i] *= GD->tDelta;

        if (GD->logtime[i] >= nData)
        {
            fprintf(stderr, "Ayay. It seems we're indexing out of bounds.\n");
            params->nFitPoints = i;
        }
    }

    fitFunc.f      = &gemFunc_residual2;
    fitFunc.n      = 2;
    fitFunc.params = (void*)GD;

    x  = gsl_vector_alloc (fitFunc.n);
    dx = gsl_vector_alloc (fitFunc.n);
    gsl_vector_set (x,  0, 25);
    gsl_vector_set (x,  1, 0.5);
    gsl_vector_set (dx, 0, 0.1);
    gsl_vector_set (dx, 1, 0.01);

    s = gsl_multimin_fminimizer_alloc (T, fitFunc.n);
    gsl_multimin_fminimizer_set (s, &fitFunc, x, dx);
    gsl_vector_free (x);
    gsl_vector_free (dx);

    do
    {
        iter++;
        status = gsl_multimin_fminimizer_iterate (s);

        if (status != 0)
        {
            gmx_fatal(FARGS, "Something went wrong in the iteration in minimizer %s:\n \"%s\"\n",
                      gsl_multimin_fminimizer_name(s), gsl_strerror(status));
        }

        d2         = gsl_multimin_fminimizer_minimum(s);
        size       = gsl_multimin_fminimizer_size(s);
        params->ka = gsl_vector_get (s->x, 0);
        params->kd = gsl_vector_get (s->x, 1);

        if (status)
        {
            fprintf(stderr, "%s\n", gsl_strerror(status));
            break;
        }

        status = gsl_multimin_test_size(size, tol);

        if (status == GSL_SUCCESS)
        {
            fprintf(stdout, "Converged to minimum at\n");
        }

        printf ("iter %5d: ka = %2.5f  kd = %2.5f  f() = %7.3f  size = %.3f  chi2 = %2.5f\n",
                iter,
                params->ka,
                params->kd,
                s->fval, size, d2);

        if (iter%1 == 0)
        {
            eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params);
            /* fixGemACF(GD->ctTheory, nFitPoints); */
            sprintf(dumpname, "Iter_%i.xvg", iter);
            for (i = 0; i < GD->nData; i++)
            {
                dumpdata[i] = (real)(GD->ctTheory[i]);
                if (!gmx_isfinite(dumpdata[i]))
                {
                    gmx_fatal(FARGS, "Non-finite value in acf.");
                }
            }
            dumpN(dumpdata, GD->nData, dumpname);
        }
    }
    while ((status == GSL_CONTINUE) && (iter < maxiter));

    /*   /\* Calculate the theoretical ACF from the parameters one last time. *\/ */
    eq10v2(GD->ctTheory, time, nData, params->ka, params->kd, params);
    *ctFit = GD->ctTheory;

    sfree(GD);
    gsl_multimin_fminimizer_free (s);

    return d2;

#endif /* HAVE_LIBGSL */
}

#ifdef HAVE_LIBGSL
static int balFunc_f(const gsl_vector *x, void *data, gsl_vector *f)
{
    /* C + sum{ A_i * exp(-B_i * t) }*/

    balData *BD;
    int      n, i, j, nexp;
    double  *y, *A, *B, C,     /* There are the parameters to be optimized. */
             t, ct;

    BD   = (balData *)data;
    n    = BD->n;
    nexp = BD->nexp;
    y    = BD->y,
    snew(A, nexp);
    snew(B, nexp);

    for (i = 0; i < nexp; i++)
    {
        A[i] = gsl_vector_get(x, i*2);
        B[i] = gsl_vector_get(x, i*2+1);
    }
    C = gsl_vector_get(x, nexp*2);

    for (i = 0; i < n; i++)
    {
        t  = i*BD->tDelta;
        ct = 0;
        for (j = 0; j < nexp; j++)
        {
            ct += A[j] * exp(B[j] * t);
        }
        ct += C;
        gsl_vector_set (f, i, ct - y[i]);
    }
    return GSL_SUCCESS;
}

/* The derivative stored in jacobian form (J)*/
static int balFunc_df(const gsl_vector *params, void *data, gsl_matrix *J)
{
    balData *BD;
    size_t   n, i, j;
    double  *y, *A, *B, C, /* There are the parameters. */
             t;
    int      nexp;

    BD   = (balData*)data;
    n    = BD->n;
    y    = BD->y;
    nexp = BD->nexp;

    snew(A, nexp);
    snew(B, nexp);

    for (i = 0; i < nexp; i++)
    {
        A[i] = gsl_vector_get(params, i*2);
        B[i] = gsl_vector_get(params, i*2+1);
    }
    C = gsl_vector_get(params, nexp*2);
    for (i = 0; i < n; i++)
    {
        t = i*BD->tDelta;
        for (j = 0; j < nexp; j++)
        {
            gsl_matrix_set (J, i, j*2,   exp(B[j]*t));        /* df(t)/dA_j */
            gsl_matrix_set (J, i, j*2+1, A[j]*t*exp(B[j]*t)); /* df(t)/dB_j */
        }
        gsl_matrix_set (J, i, nexp*2, 1);                     /* df(t)/dC */
    }
    return GSL_SUCCESS;
}

/* Calculation of the function and its derivative */
static int balFunc_fdf(const gsl_vector *params, void *data,
                       gsl_vector *f, gsl_matrix *J)
{
    balFunc_f(params, data, f);
    balFunc_df(params, data, J);
    return GSL_SUCCESS;
}
#endif /* HAVE_LIBGSL */

/* Removes the ballistic term from the beginning of the ACF,
 * just like in Omer's paper.
 */
extern void takeAwayBallistic(double *ct, double *t, int len, real tMax, int nexp, gmx_bool bDerivative)
{

    /* Use nonlinear regression with GSL instead.
     * Fit with 4 exponentials and one constant term,
     * subtract the fatest exponential. */

    int      nData, i, status, iter;
    balData *BD;
    double  *guess,           /* Initial guess. */
    *A,                       /* The fitted parameters. (A1, B1, A2, B2,... C) */
             a[2],
             ddt[2];
    gmx_bool sorted;
    size_t   n;
    size_t   p;

    nData = 0;
    do
    {
        nData++;
    }
    while (t[nData] < tMax+t[0] && nData < len);

    p = nexp*2+1;            /* Number of parameters. */

#ifdef HAVE_LIBGSL
    const gsl_multifit_fdfsolver_type *T
        = gsl_multifit_fdfsolver_lmsder;

    gsl_multifit_fdfsolver   *s;           /* The solver itself. */
    gsl_multifit_function_fdf fitFunction; /* The function to be fitted. */
    gsl_matrix               *covar;       /* Covariance matrix for the parameters.
                                            * We'll not use the result, though. */
    gsl_vector_view           theParams;

    nData = 0;
    do
    {
        nData++;
    }
    while (t[nData] < tMax+t[0] && nData < len);

    guess = NULL;
    n     = nData;

    snew(guess, p);
    snew(A, p);
    covar = gsl_matrix_alloc (p, p);

    /* Set up an initial gess for the parameters.
     * The solver is somewhat sensitive to the initial guess,
     * but this worked fine for a TIP3P box with -geminate dd
     * EDIT: In fact, this seems like a good starting pont for other watermodels too. */
    for (i = 0; i < nexp; i++)
    {
        guess[i*2]   = 0.1;
        guess[i*2+1] = -0.5 + (((double)i)/nexp - 0.5)*0.3;
    }
    guess[nexp * 2] = 0.01;

    theParams = gsl_vector_view_array(guess, p);

    snew(BD, 1);
    BD->n      = n;
    BD->y      = ct;
    BD->tDelta = t[1]-t[0];
    BD->nexp   = nexp;

    fitFunction.f      =  &balFunc_f;
    fitFunction.df     =  &balFunc_df;
    fitFunction.fdf    =  &balFunc_fdf;
    fitFunction.n      =  nData;
    fitFunction.p      =  p;
    fitFunction.params =  BD;

    s = gsl_multifit_fdfsolver_alloc (T, nData, p);
    if (s == NULL)
    {
        gmx_fatal(FARGS, "Could not set up the nonlinear solver.");
    }

    gsl_multifit_fdfsolver_set(s, &fitFunction, &theParams.vector);

    /* \=============================================/ */

    iter = 0;
    do
    {
        iter++;
        status = gsl_multifit_fdfsolver_iterate (s);

        if (status)
        {
            break;
        }
        status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
    }
    while (iter < 5000 && status == GSL_CONTINUE);

    if (iter == 5000)
    {
        fprintf(stderr, "The non-linear fitting did not converge in 5000 steps.\n"
                "Check the quality of the fit!\n");
    }
    else
    {
        fprintf(stderr, "Non-linear fitting of ballistic term converged in %d steps.\n\n", (int)iter);
    }
    for (i = 0; i < nexp; i++)
    {
        fprintf(stdout, "%c * exp(%c * t) + ", 'A'+(char)i*2, 'B'+(char)i*2);
    }

    fprintf(stdout, "%c\n", 'A'+(char)nexp*2);
    fprintf(stdout, "Here are the actual numbers for A-%c:\n", 'A'+nexp*2);

    for (i = 0; i < nexp; i++)
    {
        A[i*2]   = gsl_vector_get(s->x, i*2);
        A[i*2+1] = gsl_vector_get(s->x, i*2+1);
        fprintf(stdout, " %g*exp(%g * x) +", A[i*2], A[i*2+1]);
    }
    A[i*2] = gsl_vector_get(s->x, i*2);        /* The last and constant term */
    fprintf(stdout, " %g\n", A[i*2]);

    fflush(stdout);

    /* Implement some check for parameter quality */
    for (i = 0; i < nexp; i++)
    {
        if (A[i*2] < 0 || A[i*2] > 1)
        {
            fprintf(stderr, "WARNING: ----------------------------------\n"
                    " | A coefficient does not lie within [0,1].\n"
                    " | This may or may not be a problem.\n"
                    " | Double check the quality of the fit!\n");
        }
        if (A[i*2+1] > 0)
        {
            fprintf(stderr, "WARNING: ----------------------------------\n"
                    " | One factor in the exponent is positive.\n"
                    " | This could be a problem if the coefficient\n"
                    " | is large. Double check the quality of the fit!\n");
        }
    }
    if (A[i*2] < 0 || A[i*2] > 1)
    {
        fprintf(stderr, "WARNING: ----------------------------------\n"
                " | The constant term does not lie within [0,1].\n"
                " | This may or may not be a problem.\n"
                " | Double check the quality of the fit!\n");
    }

    /* Sort the terms */
    sorted = (nexp > 1) ?  FALSE : TRUE;
    while (!sorted)
    {
        sorted = TRUE;
        for (i = 0; i < nexp-1; i++)
        {
            ddt[0] = A[i*2] * A[i*2+1];
            ddt[1] = A[i*2+2] * A[i*2+3];

            if ((bDerivative && (ddt[0] < 0 && ddt[1] < 0 && ddt[0] > ddt[1])) || /* Compare derivative at t=0... */
                (!bDerivative && (A[i*2+1] > A[i*2+3])))                          /* Or just the coefficient in the exponent */
            {
                sorted = FALSE;
                a[0]   = A[i*2];   /* coefficient */
                a[1]   = A[i*2+1]; /* parameter in the exponent */

                A[i*2]   = A[i*2+2];
                A[i*2+1] = A[i*2+3];

                A[i*2+2] = a[0];
                A[i*2+3] = a[1];
            }
        }
    }

    /* Subtract the fastest component */
    fprintf(stdout, "Fastest component is %g * exp(%g * t)\n"
            "Subtracting fastest component from ACF.\n", A[0], A[1]);

    for (i = 0; i < len; i++)
    {
        ct[i] = (ct[i] - A[0] * exp(A[1] * i*BD->tDelta)) / (1-A[0]);
    }

    sfree(guess);
    sfree(A);

    gsl_multifit_fdfsolver_free(s);
    gsl_matrix_free(covar);
    fflush(stdout);

#else
    /* We have no gsl. */
    fprintf(stderr, "Sorry, can't take away ballistic component without gsl. "
            "Recompile using --with-gsl.\n");
    return;
#endif /* HAVE_LIBGSL */

}

extern void dumpN(const real *e, const int nn, char *fn)
{
    /* For debugging only */
    int   i;
    FILE *f;
    char  standardName[] = "Nt.xvg";
    if (fn == NULL)
    {
        fn = standardName;
    }

    f = fopen(fn, "w");
    fprintf(f,
            "@ type XY\n"
            "@ xaxis label \"Frame\"\n"
            "@ yaxis label \"N\"\n"
            "@ s0 line type 3\n");

    for (i = 0; i < nn; i++)
    {
        fprintf(f, "%-10i %-g\n", i, e[i]);
    }

    fclose(f);
}

static real* d2r(const double *d, const int nn)
{
    real *r;
    int   i;

    snew(r, nn);
    for (i = 0; i < nn; i++)
    {
        r[i] = (real)d[i];
    }

    return r;
}

static void _patchBad(double *ct, int n, double dy)
{
    /* Just do lin. interpolation for now. */
    int i;

    for (i = 1; i < n; i++)
    {
        ct[i] = ct[0]+i*dy;
    }
}

static void patchBadPart(double *ct, int n)
{
    _patchBad(ct, n, (ct[n] - ct[0])/n);
}

static void patchBadTail(double *ct, int n)
{
    _patchBad(ct+1, n-1, ct[1]-ct[0]);

}

extern void fixGemACF(double *ct, int len)
{
    int      i, j, b, e;
    gmx_bool bBad;

    /* Let's separate two things:
     * - identification of bad parts
     * - patching of bad parts.
     */

    b    = 0; /* Start of a bad stretch */
    e    = 0; /* End of a bad stretch */
    bBad = FALSE;

    /* An acf of binary data must be one at t=0. */
    if (abs(ct[0]-1.0) > 1e-6)
    {
        ct[0] = 1.0;
        fprintf(stderr, "|ct[0]-1.0| = %1.6d. Setting ct[0] to 1.0.\n", abs(ct[0]-1.0));
    }

    for (i = 0; i < len; i++)
    {

#ifdef HAS_ISFINITE
        if (isfinite(ct[i]))
#elif defined(HAS__ISFINITE)
        if (_isfinite(ct[i]))
#else
        if (1)
#endif
        {
            if (!bBad)
            {
                /* Still on a good stretch. Proceed.*/
                continue;
            }

            /* Patch up preceding bad stretch. */
            if (i == (len-1))
            {
                /* It's the tail */
                if (b <= 1)
                {
                    gmx_fatal(FARGS, "The ACF is mostly NaN or Inf. Aborting.");
                }
                patchBadTail(&(ct[b-2]), (len-b)+1);
            }

            e = i;
            patchBadPart(&(ct[b-1]), (e-b)+1);
            bBad = FALSE;
        }
        else
        {
            if (!bBad)
            {
                b = i;

                bBad = TRUE;
            }
        }
    }
}


text:
/* matrix/gsl_matrix_float.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 *
 * 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_MATRIX_FLOAT_H__
#define __GSL_MATRIX_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.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

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  float * data;
  gsl_block_float * block;
  int owner;
} gsl_matrix_float;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_view;

typedef _gsl_matrix_float_view gsl_matrix_float_view;

typedef struct
{
  gsl_matrix_float matrix;
} _gsl_matrix_float_const_view;

typedef const _gsl_matrix_float_const_view gsl_matrix_float_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_float *
gsl_matrix_float_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_float *
gsl_matrix_float_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_float *
gsl_matrix_float_alloc_from_block (gsl_block_float * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_float *
gsl_matrix_float_alloc_from_matrix (gsl_matrix_float * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_float *
gsl_vector_float_alloc_row_from_matrix (gsl_matrix_float * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_float *
gsl_vector_float_alloc_col_from_matrix (gsl_matrix_float * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_float_free (gsl_matrix_float * m);

/* Views */

GSL_EXPORT
_gsl_matrix_float_view
gsl_matrix_float_submatrix (gsl_matrix_float * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_float_view
gsl_matrix_float_row (gsl_matrix_float * m, const size_t i);

GSL_EXPORT
_gsl_vector_float_view
gsl_matrix_float_column (gsl_matrix_float * m, const size_t j);

GSL_EXPORT
_gsl_vector_float_view
gsl_matrix_float_diagonal (gsl_matrix_float * m);

GSL_EXPORT
_gsl_vector_float_view
gsl_matrix_float_subdiagonal (gsl_matrix_float * m, const size_t k);

GSL_EXPORT
_gsl_vector_float_view
gsl_matrix_float_superdiagonal (gsl_matrix_float * m, const size_t k);

GSL_EXPORT
_gsl_matrix_float_view
gsl_matrix_float_view_array (float * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_float_view
gsl_matrix_float_view_array_with_tda (float * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_float_view
gsl_matrix_float_view_vector (gsl_vector_float * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_float_view
gsl_matrix_float_view_vector_with_tda (gsl_vector_float * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_float_const_view
gsl_matrix_float_const_submatrix (const gsl_matrix_float * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_float_const_view
gsl_matrix_float_const_row (const gsl_matrix_float * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_float_const_view
gsl_matrix_float_const_column (const gsl_matrix_float * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_float_const_view
gsl_matrix_float_const_diagonal (const gsl_matrix_float * m);

GSL_EXPORT
_gsl_vector_float_const_view
gsl_matrix_float_const_subdiagonal (const gsl_matrix_float * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_float_const_view
gsl_matrix_float_const_superdiagonal (const gsl_matrix_float * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_float_const_view
gsl_matrix_float_const_view_array (const float * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_float_const_view
gsl_matrix_float_const_view_array_with_tda (const float * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector (const gsl_vector_float * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_float_const_view
gsl_matrix_float_const_view_vector_with_tda (const gsl_vector_float * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT float   gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x);

GSL_EXPORT float * gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_EXPORT const float * gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_float_set_zero (gsl_matrix_float * m);
GSL_EXPORT void gsl_matrix_float_set_identity (gsl_matrix_float * m);
GSL_EXPORT void gsl_matrix_float_set_all (gsl_matrix_float * m, float x);

GSL_EXPORT int gsl_matrix_float_fread (FILE * stream, gsl_matrix_float * m) ;
GSL_EXPORT int gsl_matrix_float_fwrite (FILE * stream, const gsl_matrix_float * m) ;
GSL_EXPORT int gsl_matrix_float_fscanf (FILE * stream, gsl_matrix_float * m);
GSL_EXPORT int gsl_matrix_float_fprintf (FILE * stream, const gsl_matrix_float * m, const char * format);

GSL_EXPORT int gsl_matrix_float_memcpy(gsl_matrix_float * dest, const gsl_matrix_float * src);
GSL_EXPORT int gsl_matrix_float_swap(gsl_matrix_float * m1, gsl_matrix_float * m2);

GSL_EXPORT int gsl_matrix_float_swap_rows(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_float_swap_columns(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_float_swap_rowcol(gsl_matrix_float * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_float_transpose (gsl_matrix_float * m);
GSL_EXPORT int gsl_matrix_float_transpose_memcpy (gsl_matrix_float * dest, const gsl_matrix_float * src);

GSL_EXPORT float gsl_matrix_float_max (const gsl_matrix_float * m);
GSL_EXPORT float gsl_matrix_float_min (const gsl_matrix_float * m);
GSL_EXPORT void gsl_matrix_float_minmax (const gsl_matrix_float * m, float * min_out, float * max_out);

GSL_EXPORT void gsl_matrix_float_max_index (const gsl_matrix_float * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_float_min_index (const gsl_matrix_float * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_float_minmax_index (const gsl_matrix_float * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_float_isnull (const gsl_matrix_float * m);

GSL_EXPORT int gsl_matrix_float_add (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_EXPORT int gsl_matrix_float_sub (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_EXPORT int gsl_matrix_float_mul_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_EXPORT int gsl_matrix_float_div_elements (gsl_matrix_float * a, const gsl_matrix_float * b);
GSL_EXPORT int gsl_matrix_float_scale (gsl_matrix_float * a, const double x);
GSL_EXPORT int gsl_matrix_float_add_constant (gsl_matrix_float * a, const double x);
GSL_EXPORT int gsl_matrix_float_add_diagonal (gsl_matrix_float * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_float_get_row(gsl_vector_float * v, const gsl_matrix_float * m, const size_t i);
GSL_EXPORT int gsl_matrix_float_get_col(gsl_vector_float * v, const gsl_matrix_float * m, const size_t j);
GSL_EXPORT int gsl_matrix_float_set_row(gsl_matrix_float * m, const size_t i, const gsl_vector_float * v);
GSL_EXPORT int gsl_matrix_float_set_col(gsl_matrix_float * m, const size_t j, const gsl_vector_float * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
float
gsl_matrix_float_get(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_float_set(gsl_matrix_float * m, const size_t i, const size_t j, const float x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
float *
gsl_matrix_float_ptr(gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (float *) (m->data + (i * m->tda + j)) ;
}

extern inline
const float *
gsl_matrix_float_const_ptr(const gsl_matrix_float * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const float *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_FLOAT_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/* vector/gsl_vector_short.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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.
 */

#ifndef __GSL_VECTOR_SHORT_H__
#define __GSL_VECTOR_SHORT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_short.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

typedef struct
{
  size_t size;
  size_t stride;
  short *data;
  gsl_block_short *block;
  int owner;
}
gsl_vector_short;

typedef struct
{
  gsl_vector_short vector;
} _gsl_vector_short_view;

typedef _gsl_vector_short_view gsl_vector_short_view;

typedef struct
{
  gsl_vector_short vector;
} _gsl_vector_short_const_view;

typedef const _gsl_vector_short_const_view gsl_vector_short_const_view;

/* Allocation */

GSL_FUN gsl_vector_short *gsl_vector_short_alloc (const size_t n);
GSL_FUN gsl_vector_short *gsl_vector_short_calloc (const size_t n);

GSL_FUN gsl_vector_short *gsl_vector_short_alloc_from_block (gsl_block_short * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_short *gsl_vector_short_alloc_from_vector (gsl_vector_short * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_short_free (gsl_vector_short * v);

/* Views */

GSL_FUN _gsl_vector_short_view
gsl_vector_short_view_array (short *v, size_t n);

GSL_FUN _gsl_vector_short_view
gsl_vector_short_view_array_with_stride (short *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_short_const_view
gsl_vector_short_const_view_array (const short *v, size_t n);

GSL_FUN _gsl_vector_short_const_view
gsl_vector_short_const_view_array_with_stride (const short *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_short_view
gsl_vector_short_subvector (gsl_vector_short *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_short_view
gsl_vector_short_subvector_with_stride (gsl_vector_short *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_short_const_view
gsl_vector_short_const_subvector (const gsl_vector_short *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_short_const_view
gsl_vector_short_const_subvector_with_stride (const gsl_vector_short *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_short_set_zero (gsl_vector_short * v);
GSL_FUN void gsl_vector_short_set_all (gsl_vector_short * v, short x);
GSL_FUN int gsl_vector_short_set_basis (gsl_vector_short * v, size_t i);

GSL_FUN int gsl_vector_short_fread (FILE * stream, gsl_vector_short * v);
GSL_FUN int gsl_vector_short_fwrite (FILE * stream, const gsl_vector_short * v);
GSL_FUN int gsl_vector_short_fscanf (FILE * stream, gsl_vector_short * v);
GSL_FUN int gsl_vector_short_fprintf (FILE * stream, const gsl_vector_short * v,
                              const char *format);

GSL_FUN int gsl_vector_short_memcpy (gsl_vector_short * dest, const gsl_vector_short * src);

GSL_FUN int gsl_vector_short_reverse (gsl_vector_short * v);

GSL_FUN int gsl_vector_short_swap (gsl_vector_short * v, gsl_vector_short * w);
GSL_FUN int gsl_vector_short_swap_elements (gsl_vector_short * v, const size_t i, const size_t j);

GSL_FUN short gsl_vector_short_max (const gsl_vector_short * v);
GSL_FUN short gsl_vector_short_min (const gsl_vector_short * v);
GSL_FUN void gsl_vector_short_minmax (const gsl_vector_short * v, short * min_out, short * max_out);

GSL_FUN size_t gsl_vector_short_max_index (const gsl_vector_short * v);
GSL_FUN size_t gsl_vector_short_min_index (const gsl_vector_short * v);
GSL_FUN void gsl_vector_short_minmax_index (const gsl_vector_short * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_short_add (gsl_vector_short * a, const gsl_vector_short * b);
GSL_FUN int gsl_vector_short_sub (gsl_vector_short * a, const gsl_vector_short * b);
GSL_FUN int gsl_vector_short_mul (gsl_vector_short * a, const gsl_vector_short * b);
GSL_FUN int gsl_vector_short_div (gsl_vector_short * a, const gsl_vector_short * b);
GSL_FUN int gsl_vector_short_scale (gsl_vector_short * a, const double x);
GSL_FUN int gsl_vector_short_add_constant (gsl_vector_short * a, const double x);

GSL_FUN int gsl_vector_short_isnull (const gsl_vector_short * v);
GSL_FUN int gsl_vector_short_ispos (const gsl_vector_short * v);
GSL_FUN int gsl_vector_short_isneg (const gsl_vector_short * v);
GSL_FUN int gsl_vector_short_isnonneg (const gsl_vector_short * v);

GSL_FUN INLINE_DECL short gsl_vector_short_get (const gsl_vector_short * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_short_set (gsl_vector_short * v, const size_t i, short x);
GSL_FUN INLINE_DECL short * gsl_vector_short_ptr (gsl_vector_short * v, const size_t i);
GSL_FUN INLINE_DECL const short * gsl_vector_short_const_ptr (const gsl_vector_short * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
short
gsl_vector_short_get (const gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_short_set (gsl_vector_short * v, const size_t i, short x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
short *
gsl_vector_short_ptr (gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (short *) (v->data + i * v->stride);
}

INLINE_FUN
const short *
gsl_vector_short_const_ptr (const gsl_vector_short * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const short *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_SHORT_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#ifndef mooney_rivlin_h
#define mooney_rivlin_h

#include <petsc.h>
#include "../include/structs.h"

#ifndef PHYSICS_STRUCT_MR
#define PHYSICS_STRUCT_MR
typedef struct Physics_MR_ *Physics_MR;
struct Physics_MR_ {
  // Material properties for MR
  CeedScalar mu_1;
  CeedScalar mu_2;
  CeedScalar lambda;
};
#endif // PHYSICS_STRUCT_MR

// Create context object
PetscErrorCode PhysicsContext_MR(MPI_Comm comm, Ceed ceed, Units *units,
                                 CeedQFunctionContext *ctx);
PetscErrorCode PhysicsSmootherContext_MR(MPI_Comm comm, Ceed ceed,
    CeedQFunctionContext ctx, CeedQFunctionContext *ctx_smoother);

// Process physics options - Mooney-Rivlin
PetscErrorCode ProcessPhysics_MR(MPI_Comm comm, Physics_MR phys, Units units);

#endif // mooney_rivlin_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
#ifndef SVGP_H
#define SVGP_H

#define M_PI 3.14159265358979323846

#include "Vector.h"
#include "Matrix.h"
#include "GaussianProcess.h"
#include <gsl/gsl_math.h>

class SVGP //: public GaussianProcess
	{
	protected:
		Matrix X; ///< Samples (N x M)
		Vector Y; ///< Output (N x 1)
		int N; ///< Number of samples
		int d;

		// Kernel parameters
		real noise_var;
		real sig_var;
		Vector scale_length;

		/// SVI parameters
		int num_inducing; ///< inducing inputs (J)
		Matrix Z; ///< Hidden variables (N x J)
		Vector u; ///< Global variables (N x 1) (targets)

		// variational distribution parametrization
		Matrix S;
		Vector m;

		real l; //step length
		int samples; //not currently used.. but should be

		Vector p_mean;
		Matrix p_var;

		Matrix Kmm;
		Matrix Knm;
		Matrix Kmn;
		Matrix Knn;
		Matrix invKmm;
		Matrix K_tilde;

		real Beta;
		real KL;
		real L;

        Matrix currentSample;
        Vector currentObservation;

        //preferably these could be done in minibatches but only with single examples for now
        //the example is repeated <subsamples> times as in Hoffman et al [2013]
        virtual void local_update(); //updating Z (local/latent variables)
        virtual void global_update(const Matrix& X_samples, const Vector& Y_samples); //updating m, S (which parametrizes q(u) and in turn gives global param u)
		//virtual Vector optimize_Z(int max_iters);
        virtual void init();
        //virtual real LogLikelihood(double* data); //computes the likelihood given a Z*, to use for optimizing the position of Z
        virtual void getSamples(Matrix& X_, Vector& Y_);
	public:
		SVGP(Matrix& X, Vector& Y, Matrix& Z, real noise_var, real sig_var, Vector scale_length, int samples);
        //initializes Z through k-means
		SVGP(Matrix& X, Vector& Y, int num_inducing, real noise_var, real sig_var, Vector scale_length, int samples);
		virtual Matrix Kernel(const Matrix& A, const Matrix& B);
		virtual void Prediction(const Vector& x, real& mean, real& var);
		virtual void UpdateGaussianProcess(); //update
		virtual void FullUpdateGaussianProcess();
		virtual real LogLikelihood();
        //virtual real LogLikelihood(const gsl_vector *v);
		virtual void AddObservation(const Vector& x, const real& y);
		virtual void AddObservation(const std::vector<Vector>& x, const std::vector<real>& y);
		virtual void Clear();

};
#endif