text:
/* multilarge/test.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 <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multilarge.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>

static void test_random_matrix_orth(gsl_matrix *m, const gsl_rng *r);
static void test_random_matrix_ill(gsl_matrix *m, const gsl_rng *r);
static void test_random_vector(gsl_vector *v, const gsl_rng *r,
                               const double lower, const double upper);
static void test_random_matrix(gsl_matrix *m, const gsl_rng *r,
                               const double lower, const double upper);
static void test_random_vector_noise(const gsl_rng *r, gsl_vector *y);
static void test_compare_vectors(const double tol, const gsl_vector * a,
                                 const gsl_vector * b, const char * desc);
static void test_multifit_solve(const double lambda, const gsl_matrix * X,
                                const gsl_vector * y, const gsl_vector * wts,
                                const gsl_vector * diagL, const gsl_matrix * L,
                                double *rnorm, double *snorm, gsl_vector * c);
static void test_multilarge_solve(const gsl_multilarge_linear_type * T, const double lambda,
                                  const gsl_matrix * X, const gsl_vector * y, gsl_vector * wts,
                                  const gsl_vector * diagL, const gsl_matrix * L,
                                  double *rnorm, double *snorm, gsl_vector * c);

/* generate random square orthogonal matrix via QR decomposition */
static void
test_random_matrix_orth(gsl_matrix *m, const gsl_rng *r)
{
  const size_t M = m->size1;
  gsl_matrix *A = gsl_matrix_alloc(M, M);
  gsl_vector *tau = gsl_vector_alloc(M);
  gsl_matrix *R = gsl_matrix_alloc(M, M);

  test_random_matrix(A, r, -1.0, 1.0);
  gsl_linalg_QR_decomp(A, tau);
  gsl_linalg_QR_unpack(A, tau, m, R);

  gsl_matrix_free(A);
  gsl_matrix_free(R);
  gsl_vector_free(tau);
}

/* construct ill-conditioned matrix via SVD */
static void
test_random_matrix_ill(gsl_matrix *m, const gsl_rng *r)
{
  const size_t M = m->size1;
  const size_t N = m->size2;
  gsl_matrix *U = gsl_matrix_alloc(M, M);
  gsl_matrix *V = gsl_matrix_alloc(N, N);
  gsl_vector *S = gsl_vector_alloc(N);
  gsl_matrix_view Uv = gsl_matrix_submatrix(U, 0, 0, M, N);
  const double smin = 16.0 * GSL_DBL_EPSILON;
  const double smax = 10.0;
  const double ratio = pow(smin / smax, 1.0 / (N - 1.0));
  double s;
  size_t j;

  test_random_matrix_orth(U, r);
  test_random_matrix_orth(V, r);

  /* compute U * S */

  s = smax;
  for (j = 0; j < N; ++j)
    {
      gsl_vector_view uj = gsl_matrix_column(U, j);

      gsl_vector_scale(&uj.vector, s);
      s *= ratio;
    }

  /* compute m = (U * S) * V' */
  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, &Uv.matrix, V, 0.0, m);

  gsl_matrix_free(U);
  gsl_matrix_free(V);
  gsl_vector_free(S);
}

static void
test_random_vector(gsl_vector *v, const gsl_rng *r,
                   const double lower, const double upper)
{
  size_t i;
  size_t N = v->size;

  for (i = 0; i < N; ++i)
    {
      gsl_vector_set(v, i,
                     gsl_rng_uniform(r) * (upper - lower) + lower);
    }
}

static void
test_random_matrix(gsl_matrix *m, const gsl_rng *r,
                   const double lower, const double upper)
{
  size_t i, j;
  size_t M = m->size1;
  size_t N = m->size2;

  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
      {
        gsl_matrix_set(m, i, j,
                       gsl_rng_uniform(r) * (upper - lower) + lower);
      }
    }
}

/* generate Vandermonde matrix using equally spaced input points
 * on [0,1] */
static void
test_vander_matrix(gsl_matrix * m)
{
  const size_t M = m->size1;
  const size_t N = m->size2;
  const double dt = 1.0 / (M - 1.0);
  size_t i, j;

  for (i = 0; i < M; ++i)
    {
      double ti = i * dt;
      double mij = 1.0;

      for (j = 0; j < N; ++j)
        {
          gsl_matrix_set(m, i, j, mij);
          mij *= ti;
        }
    }
}

static void
test_random_vector_noise(const gsl_rng *r, gsl_vector *y)
{
  size_t i;

  for (i = 0; i < y->size; ++i)
    {
      double *ptr = gsl_vector_ptr(y, i);
      *ptr += 1.0e-3 * gsl_rng_uniform(r);
    }
}

static void
test_compare_vectors(const double tol, const gsl_vector * a,
                     const gsl_vector * b, const char * desc)
{
  size_t i;

  for (i = 0; i < a->size; ++i)
    {
      double ai = gsl_vector_get(a, i);
      double bi = gsl_vector_get(b, i);

      gsl_test_rel(bi, ai, tol, "%s i=%zu", desc, i);
    }
}

/* solve least squares system with multifit SVD */
static void
test_multifit_solve(const double lambda, const gsl_matrix * X,
                    const gsl_vector * y, const gsl_vector * wts,
                    const gsl_vector * diagL, const gsl_matrix * L,
                    double *rnorm, double *snorm, gsl_vector * c)
{
  const size_t n = X->size1;
  const size_t p = X->size2;
  gsl_multifit_linear_workspace *w =
    gsl_multifit_linear_alloc(n, p);
  gsl_matrix *Xs = gsl_matrix_alloc(n, p);
  gsl_vector *ys = gsl_vector_alloc(n);
  gsl_vector *cs = gsl_vector_alloc(p);
  gsl_matrix *LQR = NULL;
  gsl_vector *Ltau = NULL;
  gsl_matrix *M = NULL;

  /* convert to standard form */
  if (diagL)
    {
      gsl_multifit_linear_wstdform1(diagL, X, wts, y, Xs, ys, w);
    }
  else if (L)
    {
      const size_t m = L->size1;

      LQR = gsl_matrix_alloc(m, p);
      Ltau = gsl_vector_alloc(GSL_MIN(m, p));
      M = (m >= p) ? gsl_matrix_alloc(m, p) : gsl_matrix_alloc(n, p);

      gsl_matrix_memcpy(LQR, L);
      gsl_multifit_linear_L_decomp(LQR, Ltau);
      gsl_multifit_linear_wstdform2(LQR, Ltau, X, wts, y, Xs, ys, M, w);
    }
  else
    {
      gsl_matrix_memcpy(Xs, X);
      gsl_vector_memcpy(ys, y);
    }

  gsl_multifit_linear_svd(Xs, w);
  gsl_multifit_linear_solve(lambda, Xs, ys, cs, rnorm, snorm, w);

  /* convert to general form */
  if (diagL)
    gsl_multifit_linear_genform1(diagL, cs, c, w);
  else if (L)
    gsl_multifit_linear_wgenform2(LQR, Ltau, X, wts, y, cs, M, c, w);
  else
    gsl_vector_memcpy(c, cs);

  gsl_multifit_linear_free(w);
  gsl_matrix_free(Xs);
  gsl_vector_free(ys);
  gsl_vector_free(cs);

  if (LQR)
    gsl_matrix_free(LQR);
  if (Ltau)
    gsl_vector_free(Ltau);
  if (M)
    gsl_matrix_free(M);
}

/* solve least squares system with multilarge */
static void
test_multilarge_solve(const gsl_multilarge_linear_type * T, const double lambda,
                      const gsl_matrix * X, const gsl_vector * y, gsl_vector * wts,
                      const gsl_vector * diagL, const gsl_matrix * L,
                      double *rnorm, double *snorm, gsl_vector * c)
{
  const size_t n = X->size1;
  const size_t p = X->size2;
  const size_t nblock = 5;
  const size_t nrows = n / nblock; /* number of rows per block */
  gsl_multilarge_linear_workspace *w =
    gsl_multilarge_linear_alloc(T, p);
  gsl_matrix *Xs = gsl_matrix_alloc(nrows, p);
  gsl_vector *ys = gsl_vector_alloc(nrows);
  gsl_vector *cs = gsl_vector_alloc(p);
  gsl_matrix *LQR = NULL;
  gsl_vector *Ltau = NULL;
  size_t rowidx = 0;

  if (L)
    {
      const size_t m = L->size1;

      LQR = gsl_matrix_alloc(m, p);
      Ltau = gsl_vector_alloc(p);

      gsl_matrix_memcpy(LQR, L);
      gsl_multilarge_linear_L_decomp(LQR, Ltau);
    }

  while (rowidx < n)
    {
      size_t nleft = n - rowidx;
      size_t nr = GSL_MIN(nrows, nleft);
      gsl_matrix_const_view Xv = gsl_matrix_const_submatrix(X, rowidx, 0, nr, p);
      gsl_vector_const_view yv = gsl_vector_const_subvector(y, rowidx, nr);
      gsl_vector_view wv;
      gsl_matrix_view Xsv = gsl_matrix_submatrix(Xs, 0, 0, nr, p);
      gsl_vector_view ysv = gsl_vector_subvector(ys, 0, nr);

      if (wts)
        wv = gsl_vector_subvector(wts, rowidx, nr);

      /* convert to standard form */
      if (diagL)
        {
          gsl_multilarge_linear_wstdform1(diagL, &Xv.matrix, wts ? &wv.vector : NULL,
                                          &yv.vector, &Xsv.matrix, &ysv.vector, w);
        }
      else if (L)
        {
          gsl_multilarge_linear_wstdform2(LQR, Ltau, &Xv.matrix, wts ? &wv.vector : NULL,
                                          &yv.vector, &Xsv.matrix, &ysv.vector, w);
        }
      else
        {
          gsl_matrix_memcpy(&Xsv.matrix, &Xv.matrix);
          gsl_vector_memcpy(&ysv.vector, &yv.vector);
        }

      gsl_multilarge_linear_accumulate(&Xsv.matrix, &ysv.vector, w);

      rowidx += nr;
    }

  gsl_multilarge_linear_solve(lambda, cs, rnorm, snorm, w);

  if (diagL)
    gsl_multilarge_linear_genform1(diagL, cs, c, w);
  else if (L)
    gsl_multilarge_linear_genform2(LQR, Ltau, cs, c, w);
  else
    gsl_vector_memcpy(c, cs);

  gsl_multilarge_linear_free(w);
  gsl_matrix_free(Xs);
  gsl_vector_free(ys);
  gsl_vector_free(cs);

  if (LQR)
    gsl_matrix_free(LQR);
  if (Ltau)
    gsl_vector_free(Ltau);
}

static void
test_random(const gsl_multilarge_linear_type * T,
            const size_t n, const size_t p,
            const double tol,
            const gsl_rng * r)
{
  const double tol1 = 1.0e3 * tol;
  gsl_matrix *X = gsl_matrix_alloc(n, p);
  gsl_vector *y = gsl_vector_alloc(n);
  gsl_vector *c = gsl_vector_alloc(p);
  gsl_vector *w = gsl_vector_alloc(n);
  gsl_vector *diagL = gsl_vector_alloc(p);
  gsl_matrix *Lsqr = gsl_matrix_alloc(p, p);
  gsl_matrix *Ltall = gsl_matrix_alloc(5*p, p);
  gsl_vector *c0 = gsl_vector_alloc(p);
  gsl_vector *c1 = gsl_vector_alloc(p);
  double rnorm0, snorm0;
  double rnorm1, snorm1;
  char str[2048];
  size_t i;

  /* generate LS system */
  test_random_matrix_ill(X, r);
  /*test_random_matrix(X, r, -1.0, 1.0);*/
  test_random_vector(c, r, -1.0, 1.0);

  /* compute y = X c + noise */
  gsl_blas_dgemv(CblasNoTrans, 1.0, X, c, 0.0, y);
  test_random_vector_noise(r, y);

  /* random weights */
  test_random_vector(w, r, 0.0, 1.0);

  /* random diag(L) */
  test_random_vector(diagL, r, 1.0, 5.0);

  /* random square L */
  test_random_matrix(Lsqr, r, -5.0, 5.0);

  /* random tall L */
  test_random_matrix(Ltall, r, -10.0, 10.0);

  for (i = 0; i < 2; ++i)
    {
      double lambda = pow(10.0, -(double) i);

      /* unweighted with L = I */
      {
        test_multifit_solve(lambda, X, y, NULL, NULL, NULL, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, NULL, NULL, NULL, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s unweighted stdform n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* weighted, L = diag(L) */
      {
        test_multifit_solve(lambda, X, y, w, diagL, NULL, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, w, diagL, NULL, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s weighted diag(L) n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* unweighted, L = diag(L) */
      {
        test_multifit_solve(lambda, X, y, NULL, diagL, NULL, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, NULL, diagL, NULL, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s unweighted diag(L) n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* weighted, L = square */
      {
        test_multifit_solve(lambda, X, y, w, NULL, Lsqr, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, w, NULL, Lsqr, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s weighted Lsqr n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* unweighted, L = square */
      {
        test_multifit_solve(lambda, X, y, NULL, NULL, Lsqr, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, NULL, NULL, Lsqr, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s unweighted Lsqr n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* weighted, L = tall */
      {
        test_multifit_solve(lambda, X, y, w, NULL, Ltall, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, w, NULL, Ltall, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s weighted Ltall n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }

      /* unweighted, L = tall */
      {
        test_multifit_solve(lambda, X, y, NULL, NULL, Ltall, &rnorm0, &snorm0, c0);
        test_multilarge_solve(T, lambda, X, y, NULL, NULL, Ltall, &rnorm1, &snorm1, c1);

        sprintf(str, "random %s unweighted Ltall n=%zu p=%zu lambda=%g",
                T->name, n, p, lambda);
        test_compare_vectors(tol, c0, c1, str);

        gsl_test_rel(rnorm1, rnorm0, tol1, "rnorm %s", str);
        gsl_test_rel(snorm1, snorm0, tol, "snorm %s", str);
      }
    }

  gsl_matrix_free(X);
  gsl_vector_free(y);
  gsl_vector_free(c);
  gsl_vector_free(w);
  gsl_vector_free(diagL);
  gsl_matrix_free(Lsqr);
  gsl_matrix_free(Ltall);
  gsl_vector_free(c0);
  gsl_vector_free(c1);
}

int
main (void)
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  gsl_ieee_env_setup();

  {
    const double tol1 = 1.0e-8;
    const double tol2 = 1.0e-11;
    const size_t n_vals[] = { 40, 356, 501 };
    const size_t p_vals[] = { 40, 213, 345 };
    size_t i;

    for (i = 0; i < 2; ++i)
      {
        size_t n = n_vals[i];
        size_t p = p_vals[i];

        /* generate random ill-conditioned LS system and test */
        test_random(gsl_multilarge_linear_normal, n, p, tol1, r);
        test_random(gsl_multilarge_linear_tsqr, n, p, tol2, r);
      }
  }

  gsl_rng_free(r);

  exit (gsl_test_summary ());
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>

#include <StGermain/libStGermain/src/StGermain.h>
#include <StgDomain/libStgDomain/src/StgDomain.h>

#include "common-driver-utils.h"
#include "stokes_block_scaling.h"
#include "stokes_Kblock_scaling.h"
#include <StgFEM/libStgFEM/src/StgFEM.h>
#include <PICellerator/libPICellerator/src/PICellerator.h>
#include <Underworld/libUnderworld/src/Underworld.h>
#include "Solvers/SLE/src/SLE.h" /* to give the AugLagStokes_SLE type */
#include "Solvers/KSPSolvers/src/KSPSolvers.h" /* for __KSP_COMMON */

#include "BSSCR.h"

#undef __FUNCT__
#define __FUNCT__ "KSPUnscale_BSSCR"
PetscErrorCode KSPUnscale_BSSCR(KSP ksp)
{
    KSP_BSSCR  *bsscr = (KSP_BSSCR *)ksp->data;
    Mat            Amat,Pmat;
    Vec            B,X;
    Mat            ApproxS;
    MatStokesBlockScaling BA=bsscr->BA;
    PetscTruth     sym;
    MatStructure   pflag;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    X             = ksp->vec_sol;
    B             = ksp->vec_rhs;
    ApproxS       = bsscr->preconditioner->matrix;
    sym           = bsscr->DIsSym;
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
    if( BA->scaling_exists == PETSC_FALSE ){
	PetscPrintf( PETSC_COMM_WORLD, "SCALING does not EXIST in %s\n", __func__);
    }else{
	if( DEFAULT == bsscr->scaletype ){
	    PetscPrintf( PETSC_COMM_WORLD, "Operator scales (pre-Unscaling)\n");
	    BSSCR_MatStokesBlockReportOperatorScales(Amat, sym);
	    BSSCR_MatStokesBlockUnScaleSystem( BA, Amat, B, X, ApproxS, sym);
	    PetscPrintf( PETSC_COMM_WORLD, "Operator scales (post-Unscaling)\n");
	    BSSCR_MatStokesBlockReportOperatorScales(Amat, sym);
	}
	if( KONLY == bsscr->scaletype ){
	    PetscPrintf( PETSC_COMM_WORLD, "Operator scales KONLY (pre-Unscaling)\n");
	    BSSCR_MatStokesKBlockReportOperatorScales(Amat, sym);
	    BSSCR_MatStokesKBlockUnScaleSystem( BA, Amat, B, X, ApproxS, sym);
	    PetscPrintf( PETSC_COMM_WORLD, "Operator scales KONLY (post-Unscaling)\n");
	    BSSCR_MatStokesKBlockReportOperatorScales(Amat, sym);
	}
    }
    bsscr->scaled = PETSC_FALSE;
    PetscFunctionReturn(0);
}
/*
  KSPScale_BSSCR constructs scaling and applies scaling to system.
  Want to do both here when doing nonlinear iterations to keep the
  scaling up to date with each iteration.
 */
#undef __FUNCT__
#define __FUNCT__ "KSPScale_BSSCR"
PetscErrorCode KSPScale_BSSCR(KSP ksp)
{
    KSP_BSSCR  *bsscr = (KSP_BSSCR *)ksp->data;
    Mat            Amat,Pmat;
    Vec            B,X;
    Mat            ApproxS;
    MatStokesBlockScaling BA=bsscr->BA;
    PetscTruth     sym;
    MatStructure   pflag;
    PetscErrorCode ierr;

    PetscFunctionBegin;
    X             = ksp->vec_sol;
    B             = ksp->vec_rhs;
    ApproxS       = bsscr->preconditioner->matrix;
    sym           = bsscr->DIsSym;
    ierr=Stg_PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);CHKERRQ(ierr);

    if( DEFAULT == bsscr->scaletype ){
	if( BA->scaling_exists ){
	    BSSCR_MatStokesBlockDefaultBuildScaling( BA, Amat);/* rebuild scaling vectors on struct */
	    BA->scalings_have_been_inverted = PETSC_FALSE;
	}else{
	    BSSCR_MatBlock_ConstructScaling( BA, Amat, B, X );/* allocates scaling vectors then calls the above function */
	    BA->scalings_have_been_inverted = PETSC_FALSE;/* above function sets this but I am making it explicit here */
	}
	PetscPrintf( PETSC_COMM_WORLD, "Operator scales (pre-scaling)\n");
	BSSCR_MatStokesBlockReportOperatorScales(Amat, sym);
	BSSCR_MatStokesBlockScaleSystem( BA, Amat, B, X, ApproxS, sym );
	PetscPrintf( PETSC_COMM_WORLD, "Operator scales (post-scaling)\n");
	BSSCR_MatStokesBlockReportOperatorScales(Amat, sym);
	bsscr->scaled = PETSC_TRUE;
    }

    if( KONLY == bsscr->scaletype ){
	if( BA->scaling_exists ){
	    BSSCR_MatStokesKBlockDefaultBuildScaling( BA, Amat, B, X, sym);/* rebuild scaling vectors on struct */
	    BA->scalings_have_been_inverted = PETSC_FALSE;
	}else{
	    BSSCR_MatKBlock_ConstructScaling( BA, Amat, B, X, sym);/* allocates scaling vectors then calls the above function */
	    BA->scalings_have_been_inverted = PETSC_FALSE;/* above function sets this but I am making it explicit here */
	}
	PetscPrintf( PETSC_COMM_WORLD, "Operator scales KONLY (pre-scaling)\n");
	BSSCR_MatStokesKBlockReportOperatorScales(Amat, sym);
	BSSCR_MatStokesKBlockScaleSystem( BA, Amat, B, X, ApproxS, sym );
	PetscPrintf( PETSC_COMM_WORLD, "Operator scales KONLY (post-scaling)\n");
	BSSCR_MatStokesKBlockReportOperatorScales(Amat, sym);
	bsscr->scaled = PETSC_TRUE;
    }
    PetscFunctionReturn(0);
}

ext: c


text:
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include <stdbool.h>

#ifdef __APPLE__

#include <cblas.h>
#include <libiomp/omp.h>

#define set_num_threads(x) openblas_set_num_threads(x)
#define get_num_threads() openblas_get_num_threads()
#else
#include <mkl.h>
#include <omp.h>
#define set_num_threads(x) mkl_set_num_threads(x)
#define get_num_threads() mkl_get_num_threads()
#endif

double norm_inf_mat_Y11(double *X, int n);

int main(int argc, char **argv) {

    char *filename = NULL;
    int option = 0, omp_threads = 1;
    while ((option = getopt(argc, argv, "f:t:")) != -1) {
        switch (option) {
            case 'f':
                filename = optarg;
                break;
            case 't':
                omp_threads = atoi(optarg);
                break;
            default:
                printf("Usage: mpi_mexp -f string -t num_threads \n");
                return 0;
        }
    }
    if (filename == NULL) {
        printf("Usage: mpi_mexp -f string -t num_threads \n");
        return 0;
    }

    set_num_threads(omp_threads);
    omp_set_num_threads(omp_threads);

    int numTasks, rank;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numTasks);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int qProcs = (int) sqrt(numTasks);
    if (qProcs - sqrt(numTasks) != 0) {
        if (rank == 0)
            printf("np must be a square number\n");
        MPI_Finalize();
        return 0;
    }

    int q_root = qProcs - 1;
    int world_root = numTasks - 1;
    int n, dims[2] = {qProcs, qProcs}, periods[2] = {0, 0}, reorder = 0, coords[2], rc, i;
    MPI_Comm cartComm;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, &cartComm);
    MPI_Comm_rank(cartComm, &rank);
    MPI_Cart_coords(cartComm, rank, 2, coords);

    MPI_File matFile;
    rc = MPI_File_open(cartComm, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, &matFile);
    if (rc && (rank == 0)) {
        printf("Unable to open file %s\n", filename);
        fflush(stdout);
    }

    MPI_Status status;
    MPI_File_read(matFile, &n, 1, MPI_INT, &status);
    int blockSize = (n + qProcs) / qProcs;  /* number of rows in _block_ */
    int lastBlock = n + 1 - (qProcs - 1) * blockSize;

    double *LocalA = (double *) malloc(blockSize * blockSize * sizeof(double));

    MPI_Datatype readFileType;
    MPI_Type_vector(blockSize, blockSize, n, MPI_DOUBLE, &readFileType);
    MPI_Type_commit(&readFileType);
    MPI_Offset offset = (MPI_Offset) (1 * sizeof(int) +
                                      sizeof(double) * (blockSize * coords[1] + blockSize * n * coords[0]));
    MPI_File_set_view(matFile, offset, MPI_DOUBLE, readFileType,
                      "native", MPI_INFO_NULL);

    MPI_File_read_at(matFile, 0, LocalA, blockSize * blockSize, MPI_DOUBLE, &status);
    MPI_Type_free(&readFileType);
    if (coords[0] == q_root) {
        memset(LocalA + (lastBlock - 1) * blockSize, 0, (blockSize + 1 - lastBlock) * blockSize * sizeof(double));
    }
    if (coords[1] == q_root) {
        int kk = 0;
        for (kk = 0; kk < blockSize; ++kk) {
            memset(LocalA + kk * blockSize + lastBlock - 1, 0, (blockSize + 1 - lastBlock) * sizeof(double));
        }
    }

    double loc_d_var;
    double *loc_vec_var = (double *) malloc(blockSize * sizeof(double));
    memset(loc_vec_var, 0, blockSize * sizeof(double));
    offset = (MPI_Offset) (1 * sizeof(int) + sizeof(double) * n * n);
    MPI_File_set_view(matFile, offset, MPI_DOUBLE, MPI_DOUBLE,
                      "native", MPI_INFO_NULL);
    offset = coords[0] * blockSize;
    MPI_File_read_at_all(matFile, offset, loc_vec_var, coords[0] == qProcs - 1 ? lastBlock : blockSize, MPI_DOUBLE,
                         &status);

    MPI_File_close(&matFile);

    MPI_Comm rowComm, colComm;
    MPI_Comm_split(MPI_COMM_WORLD, coords[0], rank, &rowComm);
    MPI_Comm_split(MPI_COMM_WORLD, coords[1], rank, &colComm);

    int row_rank, col_rank;
    MPI_Comm_rank(colComm, &col_rank);
    MPI_Comm_rank(rowComm, &row_rank);

    if (coords[1] == qProcs - 1) {
        cblas_dcopy(blockSize, loc_vec_var, 1, LocalA + lastBlock - 1, blockSize);
    }

    double mpi_start = MPI_Wtime();
    //compute inf norm of matA
    for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
        loc_vec_var[i] = cblas_dasum((row_rank == q_root ? lastBlock - 1 : blockSize), LocalA + i * blockSize, 1);
    }

    if (row_rank == q_root) {
        MPI_Reduce(MPI_IN_PLACE, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
    } else {
        MPI_Reduce(loc_vec_var, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
    }

    if (row_rank == q_root) {
        loc_d_var = 0;
        for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
            if (loc_vec_var[i] > loc_d_var)
                loc_d_var = loc_vec_var[i];
        }
        if (col_rank == q_root) {
            MPI_Reduce(MPI_IN_PLACE, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
        } else {
            MPI_Reduce(&loc_d_var, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
        }
    }

    MPI_Bcast(&loc_d_var, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);

    //init Y martix

    for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
        cblas_dscal((row_rank == q_root ? lastBlock - 1 : blockSize), -1 / loc_d_var, LocalA + i * blockSize, 1);
    }
    if (row_rank == q_root) {
        cblas_dscal((col_rank == q_root ? lastBlock - 1 : blockSize), 1 / loc_d_var, LocalA + lastBlock - 1, blockSize);

    }
    if (row_rank == col_rank) {
        for (i = 0; i < (row_rank == q_root ? lastBlock-1 : blockSize); ++i) {
            LocalA[i * blockSize + i] += 1;
        }
    }
    if (rank == world_root) {
        LocalA[(lastBlock - 1) * (blockSize + 1)] = 1;
    }

    //loop
    double tolerance = 1e-8;
    int loopBreak = 0, iteration = 0;
    double *col_mat = (double *) malloc(qProcs * blockSize * blockSize * sizeof(double));
    memset(col_mat, 0, qProcs * blockSize * blockSize * sizeof(double));

    double *row_mat = (double *) malloc(qProcs * blockSize * blockSize * sizeof(double));
    memset(row_mat, 0, qProcs * blockSize * blockSize * sizeof(double));

    while (true) {
        //norm mat
        //compute inf norm of matA
        for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
            loc_vec_var[i] = cblas_dasum((row_rank == q_root ? lastBlock - 1 : blockSize), LocalA + i * blockSize, 1);
        }

        if (row_rank == q_root) {
            MPI_Reduce(MPI_IN_PLACE, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
        } else {
            MPI_Reduce(loc_vec_var, loc_vec_var, blockSize, MPI_DOUBLE, MPI_SUM, q_root, rowComm);
        }

        if (row_rank == q_root) {
            loc_d_var = 0;
            for (i = 0; i < (col_rank == q_root ? lastBlock - 1 : blockSize); ++i) {
                if (loc_vec_var[i] > loc_d_var)
                    loc_d_var = loc_vec_var[i];
            }
            if (col_rank == q_root) {
                MPI_Reduce(MPI_IN_PLACE, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
            } else {
                MPI_Reduce(&loc_d_var, &loc_d_var, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);
            }

            size_t index = cblas_idamax((col_rank == q_root ? lastBlock - 1 : blockSize), LocalA + lastBlock - 1,
                                        blockSize);
            double norm2b;
            MPI_Reduce(&LocalA[index * blockSize + lastBlock - 1], &norm2b, 1, MPI_DOUBLE, MPI_MAX, q_root, colComm);

            if (col_rank == q_root) {
                printf("normM:%.4f, normB:%.4f\n", loc_d_var, norm2b);
                loopBreak = (loc_d_var / norm2b < tolerance) ? 1 : 0;
            }

        }
        MPI_Bcast(&loopBreak, 1, MPI_INT, world_root, MPI_COMM_WORLD);
        if (loopBreak > 0)
            break;

        iteration++;

        //gather loc mat
        MPI_Allgather(LocalA, blockSize * blockSize, MPI_DOUBLE, col_mat, blockSize * blockSize, MPI_DOUBLE, colComm);
        MPI_Allgather(LocalA, blockSize * blockSize, MPI_DOUBLE, row_mat, blockSize * blockSize, MPI_DOUBLE, rowComm);

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, blockSize, blockSize, blockSize,
                    1, row_mat, blockSize, col_mat, blockSize, 0.0, LocalA, blockSize);
        for (i = 1; i < qProcs; ++i) {
            cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, blockSize, blockSize, blockSize,
                        1, row_mat + i * (blockSize * blockSize), blockSize, col_mat + i * (blockSize * blockSize),
                        blockSize, 1, LocalA, blockSize);
        }
    }

    if (rank == world_root) {
        printf("mpi_mexp_iter:%d \n", iteration);
        printf("mpi_mexp_time: %f\n", MPI_Wtime() - mpi_start);

    }

    if (row_rank == q_root) {
        for (i = 0; i < (row_rank == qProcs - 1 ? lastBlock : blockSize); ++i) {
            if (fabs(LocalA[i * blockSize + lastBlock - 1] - 1) > 1e-3) {
                printf("MPI_MEXP_ERR %f \n", LocalA[i * blockSize + lastBlock - 1]);
                break;
            }
        }
    }

    MPI_Comm_free(&cartComm);
    MPI_Comm_free(&colComm);
    MPI_Comm_free(&rowComm);

    MPI_Finalize();
}

max_stars_repo_name: baishuai/MEXP

max_issues_repo_name: baishuai/MEXP

max_forks_repo_name: baishuai/MEXP


text:
#ifndef STDAFX__H
#define STDAFX__H

#if defined(_WIN32)

#include <SDKDDKVer.h>

#if !defined(_STL_EXTRA_DISABLED_WARNINGS)
#define _STL_EXTRA_DISABLED_WARNINGS 4061 4324 4365 4514 4571 4582 4583 4623 4625 4626 4710 4774 4800 4820 4987 5026 5027 5039
#endif

#if !defined(_SCL_SECURE_NO_WARNINGS)
#define _SCL_SECURE_NO_WARNINGS 1
#endif

#if !defined(_CRT_SECURE_NO_WARNINGS)
#define _CRT_SECURE_NO_WARNINGS 1
#endif

#if !defined(_SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING)
#define _SILENCE_CXX17_ITERATOR_BASE_CLASS_DEPRECATION_WARNING 1
#endif

#if !defined(_SILENCE_CXX17_OLD_ALLOCATOR_MEMBERS_DEPRECATION_WARNING)
#define _SILENCE_CXX17_OLD_ALLOCATOR_MEMBERS_DEPRECATION_WARNING 1
#endif

#if !defined(_SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING)
#define _SILENCE_CXX17_CODECVT_HEADER_DEPRECATION_WARNING 1
#endif

#define STRICT
#define NOMINMAX

#pragma warning(disable: 4571) // warning C4571: Informational: catch(...) semantics changed since Visual C++ 7.1; structured exceptions (SEH) are no longer caught
#pragma warning(disable: 4668) // warning C4668: '%s' is not defined as a preprocessor macro, replacing with '0' for '#if/#elif'
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4711) // warning C4711: function '%s' selected for automatic inline expansion
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'
#pragma warning(disable: 5045) // warning C5045: Compiler will insert Spectre mitigation for memory load if /Qspectre switch specified

#include <Windows.h>

#pragma warning(push)
#pragma warning(disable: 4005) // warning C4005: '%s': macro redefinition
#include <Winternl.h>
#include <ntstatus.h>
#pragma warning(pop)

#else // defined(_WIN32)

#include <cpuid.h>

#endif

#if defined(_MSC_VER)

// currently broken (triggers on iterators and other objects that are usually unnamed)
#pragma warning(disable: 26444) // warning C26444: Avoid unnamed objects with custom construction and destruction (es.84: http://go.microsoft.com/fwlink/?linkid=862923).

// disable for everything
#pragma warning(disable: 4061) // warning C4061: enumerator '%s' in switch of enum '%s' is not explicitly handled by a case label
#pragma warning(disable: 4324) // warning C4234: structure was padded due to alignment specifier
#pragma warning(disable: 4514) // warning C4514: '%s': unreferenced inline function has been removed
#pragma warning(disable: 4623) // warning C4623: '%s': default constructor was implicitly defined as deleted
#pragma warning(disable: 4625) // warning C4625: '%s': copy constructor was implicitly defined as deleted
#pragma warning(disable: 4626) // warning C4626: '%s': assignment operator was implicitly defined as deleted
#pragma warning(disable: 4710) // warning C4710: '%s': function not inlined
#pragma warning(disable: 4820) // warning C4820: '%s': '%d' bytes padding added after data member '%s'
#pragma warning(disable: 5026) // warning C5026: '%s': move constructor was implicitly defined as deleted
#pragma warning(disable: 5027) // warning C5027: '%s': move assignment operator was implicitly defined as deleted

#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'no initialization'.
#pragma warning(disable: 26426) // warning C26426: Global initializer calls a non-constexpr function '%s' (i.22: http://go.microsoft.com/fwlink/?linkid=853919).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414).
#pragma warning(disable: 26485) // warning C26485: Expression '%s::`vbtable'': No array to pointer decay. (bounds.3: http://go.microsoft.com/fwlink/p/?LinkID=620415)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26499) // warning C26499: Could not find any lifetime tracking information for '%s'

// disable for standard headers
#pragma warning(push)
#pragma warning(disable: 26400) // warning C26400: Do not assign the result of an allocation or a function call with an owner<T> return value to a raw pointer, use owner<T> instead. (i.11 http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26401) // warning C26401: Do not delete a raw pointer that is not an owner<T>. (i.11: http://go.microsoft.com/fwlink/?linkid=845474)
#pragma warning(disable: 26408) // warning C26408: Avoid malloc() and free(), prefer the nothrow version of new with delete. (r.10 http://go.microsoft.com/fwlink/?linkid=845483)
#pragma warning(disable: 26409) // warning C26409: Avoid calling new and delete explicitly, use std::make_unique<T> instead. (r.11 http://go.microsoft.com/fwlink/?linkid=845485)
#pragma warning(disable: 26411) // warning C26411: The parameter '%s' is a reference to unique pointer and it is never reassigned or reset, use T* or T& instead. (r.33 http://go.microsoft.com/fwlink/?linkid=845479)
#pragma warning(disable: 26412) // warning C26412: Do not dereference an invalid pointer (lifetimes rule 1). 'return of %s' was invalidated at line %d by 'end of function scope (local lifetimes end)'.
#pragma warning(disable: 26413) // warning C26413: Do not dereference nullptr (lifetimes rule 2). 'nullptr' was pointed to nullptr at line %d.
#pragma warning(disable: 26423) // warning C26423: The allocation was not directly assigned to an owner.
#pragma warning(disable: 26424) // warning C26424: Failing to delete or assign ownership of allocation at line %d.
#pragma warning(disable: 26425) // warning C26425: Assigning '%s' to a static variable.
#pragma warning(disable: 26444) // warning C26444: Avoid unnamed objects with custom construction and destruction (es.84: http://go.microsoft.com/fwlink/?linkid=862923).
#pragma warning(disable: 26461) // warning C26461: The reference argument '%s' for function %s can be marked as const. (con.3: https://go.microsoft.com/fwlink/p/?LinkID=786684)
#pragma warning(disable: 26471) // warning C26471: Don't use reinterpret_cast. A cast from void* can use static_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417).
#pragma warning(disable: 26481) // warning C26481: Don't use pointer arithmetic. Use span instead. (bounds.1: http://go.microsoft.com/fwlink/p/?LinkID=620413)
#pragma warning(disable: 26482) // warning C26482: Only index into arrays using constant expressions. (bounds.2: http://go.microsoft.com/fwlink/p/?LinkID=620414)
#pragma warning(disable: 26490) // warning C26490: Don't use reinterpret_cast. (type.1: http://go.microsoft.com/fwlink/p/?LinkID=620417)
#pragma warning(disable: 26493) // warning C26493: Don't use C-style casts that would perform a static_cast downcast, const_cast, or reinterpret_cast. (type.4: http://go.microsoft.com/fwlink/p/?LinkID=620420)
#pragma warning(disable: 26494) // warning C26494: Variable '%s' is uninitialized. Always initialize an object. (type.5: http://go.microsoft.com/fwlink/p/?LinkID=620421)
#pragma warning(disable: 26495) // warning C26495: Variable '%s' is uninitialized. Always initialize a member variable. (type.6: http://go.microsoft.com/fwlink/p/?LinkID=620422)
#pragma warning(disable: 26496) // warning C26496: Variable '%s' is assigned only once, mark it as const. (con.4: https://go.microsoft.com/fwlink/p/?LinkID=784969)
#pragma warning(disable: 26497) // warning C26497: This function %s could be marked constexpr if compile-time evaluation is desired. (f.4: https://go.microsoft.com/fwlink/p/?LinkID=784970)

#else

// linux warnings go here

#endif

#include <cstddef>
#include <map>
#include <set>
#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <type_traits>
#include <utility>
#include <memory>
#include <tuple>
#include <thread>
#include <cstdlib>
#include <codecvt>

#if defined(_MSC_VER)

// disable additionally for third-party libraries
#pragma warning(push)
#pragma warning(disable: 4456) // warning C4456: declaration of '%s' hides previous local declaration
#pragma warning(disable: 4458) // warning C4458: declaration of '%s' hides class member
#pragma warning(disable: 4459) // warning C4459: declaration of '%s' hides global declaration
#pragma warning(disable: 4702) // warning C4702: unreacha3eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeble code

#include <CppCoreCheck\Warnings.h>
#pragma warning(disable: ALL_CPPCORECHECK_WARNINGS)

#endif

#ifndef BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#define BOOST_CONFIG_SUPPRESS_OUTDATED_MESSAGE
#endif

#include <boost/algorithm/string.hpp>
#include <boost/xpressive/xpressive.hpp>

#include <gsl/gsl>

#include <fmt/format.h>

#if defined(_MSC_VER)

#pragma warning(pop)
#pragma warning(pop)

#else

// linux warning restoration goes here

#endif

#endif

ext: h

max_stars_repo_name: DrPizza/cpuid

max_issues_repo_name: DrPizza/cpuid

max_forks_repo_name: DrPizza/cpuid


text:
/*
 * Simple example of NLopt usage from:
 * http://nlopt.readthedocs.io/en/latest/NLopt_Tutorial/
 */
#include <stdio.h>
#include <math.h>
#include <nlopt.h>

int count = 0;

//
// Objective function.
//
double myfunc(unsigned n, const double *x, double *grad, void *my_func_data)
{
    count++;
    if (grad) {
        grad[0] = 0.0;
        grad[1] = 0.5 / sqrt(x[1]);
    }
    return sqrt(x[1]);
}

typedef struct {
    double a, b;
} my_constraint_data;

//
// Constraint function.
//
double myconstraint(unsigned n, const double *x, double *grad, void *data)
{
    my_constraint_data *d = (my_constraint_data*) data;
    double a = d->a,
           b = d->b;
    if (grad) {
        grad[0] = 3 * a * (a*x[0] + b) * (a*x[0] + b);
        grad[1] = -1.0;
    }
    return ((a*x[0] + b) * (a*x[0] + b) * (a*x[0] + b) - x[1]);
}

int main()
{
    // Select algorithm.

    //
    // Global optimization
    //
//  nlopt_algorithm algorithm = NLOPT_GN_DIRECT;        // DIviding RECTangles algorithm for global optimization
//  nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L;      // locally biased variant of DIRECT
//  nlopt_algorithm algorithm = NLOPT_GN_DIRECT_L_RAND; // slightly randomized variant of DIRECT-L
//  nlopt_algorithm algorithm = NLOPT_GN_CRS2_LM;       // Controlled Random Search with local mutation
//  nlopt_algorithm algorithm = NLOPT_G_MLSL;           // Multi-Level Single-Linkage
//  nlopt_algorithm algorithm = NLOPT_G_MLSL_LDS;       // modification of MLSL with low-discrepancy sequence
//  nlopt_algorithm algorithm = NLOPT_GD_STOGO;         // StoGO
//  nlopt_algorithm algorithm = NLOPT_GD_STOGO_RAND;    // randomized variant of StoGO
//  nlopt_algorithm algorithm = NLOPT_GN_ISRES;         // Improved Stochastic Ranking Evolution Strategy
//  nlopt_algorithm algorithm = NLOPT_GN_ESCH;          // evolutionary algorithm

    //
    // Local derivative-free optimization
    //
//  nlopt_algorithm algorithm = NLOPT_LN_COBYLA;        // Constrained Optimization BY Linear Approximations
//  nlopt_algorithm algorithm = NLOPT_LN_BOBYQA;        // BOBYQA, iterative quadratic approximation
//  nlopt_algorithm algorithm = NLOPT_LN_NEWUOA;        // NEWUOA, quadratic approximation
//  nlopt_algorithm algorithm = NLOPT_LN_NEWUOA_BOUND;  // variant of NEWUOA with bound constraints
//  nlopt_algorithm algorithm = NLOPT_LN_PRAXIS;        // PRincipal AXIS
//  nlopt_algorithm algorithm = NLOPT_LN_NELDERMEAD;    // Nelder-Mead Simplex
//  nlopt_algorithm algorithm = NLOPT_LN_SBPLX;         // re-implementation of Subplex algorithm

    //
    // Local gradient-based optimization
    //
    nlopt_algorithm algorithm = NLOPT_LD_MMA;           // Method of Moving Asymptotes
//  nlopt_algorithm algorithm = NLOPT_LD_CCSAQ;         // conservative convex separable approximation, quadratic
//  nlopt_algorithm algorithm = NLOPT_LD_SLSQP;         // Sequential Least-Squares Quadratic Programming
//  nlopt_algorithm algorithm = NLOPT_LD_LBFGS;         // Low-storage BFGS
//  nlopt_algorithm algorithm = NLOPT_LD_TNEWTON;       // Truncated Newton
//  nlopt_algorithm algorithm = NLOPT_LD_VAR1;          // shifted limited-memory variable-metric algorithm, rank-1
//  nlopt_algorithm algorithm = NLOPT_LD_VAR2;          // shifted limited-memory variable-metric algorithm, rank-2

    // Create optimizer object.
    // Select dimensionality.
    nlopt_opt opt = nlopt_create(algorithm, 2);

    // Set objective function.
    nlopt_set_min_objective(opt, myfunc, NULL);

    // Specify bounds.
    double lb[2] = { -HUGE_VAL, 0 }; // lower bounds
    nlopt_set_lower_bounds(opt, lb);
    //TODO: nlopt_set_upper_bounds();

    // Add constraints.
    my_constraint_data data[2] = { {2, 0}, {-1, 1} };
    nlopt_add_inequality_constraint(opt, myconstraint, &data[0], 1e-8);
    nlopt_add_inequality_constraint(opt, myconstraint, &data[1], 1e-8);

    // Stopping criteria, or. a relative tolerance on the optimization parameters.
    nlopt_set_xtol_rel(opt, 1e-4);

    // Perform the optimization, starting with some initial guess.
    double x[2] = { 1.234, 5.678 }; // initial guess
    double minf;                    // the minimum objective value upon return

    if (nlopt_optimize(opt, x, &minf) < 0) {
        printf("nlopt failed!\n");
    } else {
        printf("found minimum after %d evaluations\n", count);
        printf("found minimum at f(%g, %g) = %0.10g\n", x[0], x[1], minf);
    }

    nlopt_destroy(opt);
    return 0;
}

max_stars_repo_name: sergev/vak-opensource

max_issues_repo_name: sergev/vak-opensource

max_forks_repo_name: sergev/vak-opensource


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Clemens Groepl $
// $Authors: $
// --------------------------------------------------------------------------

#ifndef OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H
#define OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H

#include <OpenMS/TRANSFORMATIONS/FEATUREFINDER/Fitter1D.h>

#include <gsl/gsl_rng.h> // gsl random number generators
#include <gsl/gsl_randist.h> // gsl random number distributions
#include <gsl/gsl_vector.h> // gsl vector and matrix definitions
#include <gsl/gsl_multifit_nlin.h> // gsl multidimensional fitting
#include <gsl/gsl_blas.h> // gsl linear algebra stuff

#include <OpenMS/DATASTRUCTURES/ListUtils.h>

namespace OpenMS
{

  /**
    @brief Abstract class for 1D-model fitter using Levenberg-Marquardt algorithm for parameter optimization.
      */
  class OPENMS_DLLAPI LevMarqFitter1D :
    public Fitter1D
  {

public:

    typedef std::vector<double> ContainerType;

    /// Default constructor
    LevMarqFitter1D() :
      Fitter1D()
    {
      this->defaults_.setValue("max_iteration", 500, "Maximum number of iterations using by Levenberg-Marquardt algorithm.", ListUtils::create<String>("advanced"));
      this->defaults_.setValue("deltaAbsError", 0.0001, "Absolute error used by the Levenberg-Marquardt algorithm.", ListUtils::create<String>("advanced"));
      this->defaults_.setValue("deltaRelError", 0.0001, "Relative error used by the Levenberg-Marquardt algorithm.", ListUtils::create<String>("advanced"));
    }

    /// copy constructor
    LevMarqFitter1D(const LevMarqFitter1D & source) :
      Fitter1D(source),
      max_iteration_(source.max_iteration_),
      abs_error_(source.abs_error_),
      rel_error_(source.rel_error_)
    {
    }

    /// destructor
    virtual ~LevMarqFitter1D()
    {
    }

    /// assignment operator
    virtual LevMarqFitter1D & operator=(const LevMarqFitter1D & source)
    {
      if (&source == this) return *this;

      Fitter1D::operator=(source);
      max_iteration_ = source.max_iteration_;
      abs_error_ = source.abs_error_;
      rel_error_ = source.rel_error_;

      return *this;
    }

protected:

    /// GSL status
    Int gsl_status_;
    /// Parameter indicates symmetric peaks
    bool symmetric_;
    /// Maximum number of iterations
    Int max_iteration_;
    /** Test for the convergence of the sequence by comparing the last iteration step dx with the absolute error epsabs and relative error epsrel to the current position x */
    /// Absolute error
    CoordinateType abs_error_;
    /// Relative error
    CoordinateType rel_error_;

    /** Display the intermediate state of the solution. The solver state contains
        the vector s->x which is the current position, and the vector s->f with
        corresponding function values */
    virtual void printState_(Int iter, gsl_multifit_fdfsolver * s) = 0;

    /// Return GSL status as string
    const String getGslStatus_()
    {
      return gsl_strerror(gsl_status_);
    }

    /**
        @brief Optimize start parameter

        @exception Exception::UnableToFit is thrown if fitting cannot be performed
    */
    void optimize_(const RawDataArrayType & set, Int num_params, CoordinateType x_init[],
                   Int (* residual)(const gsl_vector * x, void * params, gsl_vector * f),
                   Int (* jacobian)(const gsl_vector * x, void * params, gsl_matrix * J),
                   Int (* evaluate)(const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix * J),
                   void * advanced_params
                   )
    {

      const gsl_multifit_fdfsolver_type * T;
      gsl_multifit_fdfsolver * s;

      Int status;
      Int iter = 0;
      const UInt n = (UInt)set.size();

      // number of parameters to be optimized
      UInt p = num_params;

      // gsl always expects N>=p or default gsl error handler invoked,
      // cause Jacobian be rectangular M x N with M>=N
      if (n < p) throw Exception::UnableToFit(__FILE__, __LINE__, __PRETTY_FUNCTION__, "UnableToFit-FinalSet", "Skipping feature, gsl always expects N>=p");

      // allocate space for a covariance matrix of size p by p
      gsl_matrix * covar = gsl_matrix_alloc(p, p);
      gsl_multifit_function_fdf f;

      gsl_vector_view x = gsl_vector_view_array(x_init, p);

      gsl_rng_env_setup();

      // set up the function to be fit
      f.f = (residual);     // the function of residuals
      f.df = (jacobian);     // the gradient of this function
      f.fdf = (evaluate);     // combined function and gradient
      f.n = set.size();     // number of points in the data set
      f.p = p;     // number of parameters in the fit function
      f.params = advanced_params;     // // structure with the data and error bars

      T = gsl_multifit_fdfsolver_lmsder;
      s = gsl_multifit_fdfsolver_alloc(T, n, p);
      gsl_multifit_fdfsolver_set(s, &f, &x.vector);

#ifdef DEBUG_FEATUREFINDER
      printState_(iter, s);
#endif

      // this is the loop for fitting
      do
      {
        iter++;

        // perform a single iteration of the fitting routine
        status = gsl_multifit_fdfsolver_iterate(s);

#ifdef DEBUG_FEATUREFINDER
        // customized routine to print out current parameters
        printState_(iter, s);
#endif

        /* check if solver is stuck */
        if (status) break;

        // test for convergence with an absolute and relative error
        status = gsl_multifit_test_delta(s->dx, s->x, abs_error_, rel_error_);
      }
      while (status == GSL_CONTINUE && iter < max_iteration_);

      // This function uses Jacobian matrix J to compute the covariance matrix of the best-fit parameters, covar.
      // The parameter epsrel (0.0) is used to remove linear-dependent columns when J is rank deficient.
      gsl_multifit_covar(s->J, 0.0, covar);

#ifdef DEBUG_FEATUREFINDER
      gsl_matrix_fprintf(stdout, covar, "covar %g");
#endif

#define FIT(i) gsl_vector_get(s->x, i)
#define ERR(i) sqrt(gsl_matrix_get(covar, i, i))

      // Set GSl status
      gsl_status_ = status;

#ifdef DEBUG_FEATUREFINDER
      {
        // chi-squared value
        DoubleReal chi = gsl_blas_dnrm2(s->f);
        DoubleReal dof = n - p;
        DoubleReal c = GSL_MAX_DBL(1, chi / sqrt(dof));

        printf("chisq/dof = %g\n", pow(chi, 2.0) / dof);

        for (Size i = 0; i < p; ++i)
        {
          std::cout << i;
          printf(".Parameter = %.5f +/- %.5f\n", FIT(i), c * ERR(i));
        }
      }
#endif

      // set optimized parameters
      for (Size i = 0; i < p; ++i)
      {
        x_init[i] = FIT(i);
      }

      gsl_multifit_fdfsolver_free(s);
      gsl_matrix_free(covar);

    }

    void updateMembers_()
    {
      Fitter1D::updateMembers_();
      max_iteration_ = this->param_.getValue("max_iteration");
      abs_error_ = this->param_.getValue("deltaAbsError");
      rel_error_ = this->param_.getValue("deltaRelError");
    }

  };
}

#endif // OPENMS_TRANSFORMATIONS_FEATUREFINDER_LEVMARQFITTER1D_H

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: kreinert/OpenMS

max_forks_repo_name: kreinert/OpenMS


text:
/*
 *  siegert_neuron.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef SIEGERT_NEURON_H
#define SIEGERT_NEURON_H

#include "config.h"

#ifdef HAVE_GSL

// C includes:
#include <gsl/gsl_integration.h>

// Includes from nestkernel:
#include "archiving_node.h"
#include "connection.h"
#include "event.h"
#include "nest_types.h"
#include "node.h"
#include "normal_randomdev.h"
#include "poisson_randomdev.h"
#include "ring_buffer.h"
#include "recordables_map.h"
#include "universal_data_logger.h"

namespace nest
{

/* BeginDocumentation
Name: siegert_neuron

Description:

 siegert_neuron is an implementation of a rate model with the
 non-linearity given by the gain function of the
 leaky-integrate-and-fire neuron with delta or exponentially decaying
 synapses [2] and [3, their eq. 25]. The model can be used for a
 mean-field analysis of spiking networks.

 The model supports connections to other rate models with zero
 delay, and uses the secondary_event concept introduced with the
 gap-junction framework.

Parameters:

 The following parameters can be set in the status dictionary.

 rate                double - Rate (1/s)
 tau                 double - Time constant in ms.
 mean                double - Additional constant input

 The following parameters can be set in the status directory and are
 used in the evaluation of the gain function. Parameters as in
 iaf_psc_exp/delta.

 tau_m               double - Membrane time constant in ms.
 tau_syn             double - Time constant of postsynaptic currents in ms.
 t_ref               double - Duration of refractory period in ms.
 theta               double - Threshold relative to resting potential in mV.
 V_reset             double - Reset relative to resting membrane potential in
                              mV.

Notes:

References:

 [1] Hahne, J., Dahmen, D., Schuecker, J., Frommer, A.,
 Bolten, M., Helias, M. and Diesmann, M. (2017).
 Integration of Continuous-Time Dynamics in a
 Spiking Neural Network Simulator.
 Front. Neuroinform. 11:34. doi: 10.3389/fninf.2017.00034

 [2] Fourcaud, N and Brunel, N. (2002). Dynamics of the firing
 probability of noisy integrate-and-fire neurons, Neural computation,
 14:9, pp 2057--2110

 [3] Schuecker, J., Diesmann, M. and Helias, M. (2015).
 Modulated escape from a metastable state driven by colored noise.
 Physical Review E 92:052119

 [4] Hahne, J., Helias, M., Kunkel, S., Igarashi, J.,
 Bolten, M., Frommer, A. and Diesmann, M. (2015).
 A unified framework for spiking and gap-junction interactions
 in distributed neuronal network simulations.
 Front. Neuroinform. 9:22. doi: 10.3389/fninf.2015.00022

Sends: DiffusionConnectionEvent

Receives: DiffusionConnectionEvent, DataLoggingRequest

Author: Jannis Schuecker, David Dahmen, Jan Hahne
SeeAlso: diffusion_connection
*/
class siegert_neuron : public Archiving_Node
{

public:
  typedef Node base;

  siegert_neuron();
  siegert_neuron( const siegert_neuron& );

  /**
   * Import sets of overloaded virtual functions.
   * @see Technical Issues / Virtual Functions: Overriding, Overloading, and
   * Hiding
   */
  using Node::handle;
  using Node::sends_secondary_event;

  void handle( DiffusionConnectionEvent& );
  void handle( DataLoggingRequest& );

  port handles_test_event( DiffusionConnectionEvent&, rport );
  port handles_test_event( DataLoggingRequest&, rport );

  void
  sends_secondary_event( DiffusionConnectionEvent& )
  {
  }

  void get_status( DictionaryDatum& ) const;
  void set_status( const DictionaryDatum& );

private:
  void init_state_( const Node& proto );
  void init_buffers_();
  void calibrate();

  /** This is the actual update function. The additional boolean parameter
   * determines if the function is called by update (false) or wfr_update (true)
   */
  bool update_( Time const&, const long, const long, const bool );

  void update( Time const&, const long, const long );
  bool wfr_update( Time const&, const long, const long );

  // siegert functions and helpers
  double siegert( double, double );
  double siegert1( double, double, double, double );
  double siegert2( double, double, double, double );

  // The next two classes need to be friends to access the State_ class/member
  friend class RecordablesMap< siegert_neuron >;
  friend class UniversalDataLogger< siegert_neuron >;

  // ----------------------------------------------------------------

  /**
   * Independent parameters of the model.
   */
  struct Parameters_
  {
    /** Time constant in ms. */
    double tau_;

    /** Membrane time constant in ms. */
    double tau_m_;

    /** Synaptic time constant in ms. */
    double tau_syn_;

    /** Refractory period in ms. */
    double t_ref_;

    /** Constant input in Hz. */
    double mean_;

    /** Threshold in mV. */
    double theta_;

    /** reset value in mV. */
    double V_reset_;

    Parameters_(); //!< Sets default parameter values

    void get( DictionaryDatum& ) const; //!< Store current values in dictionary

    void set( const DictionaryDatum& );
  };

  // ----------------------------------------------------------------

  /**
   * State variables of the model.
   */
  struct State_
  {
    double r_; //!< Rate

    State_(); //!< Default initialization

    void get( DictionaryDatum& ) const;
    void set( const DictionaryDatum& );
  };

  // ----------------------------------------------------------------

  /**
   * Buffers of the model.
   */
  struct Buffers_
  {
    Buffers_( siegert_neuron& );
    Buffers_( const Buffers_&, siegert_neuron& );

    std::vector< double >
      drift_input_; //!< buffer for drift term received by DiffusionConnection
    std::vector< double > diffusion_input_; //!< buffer for diffusion term
    // received by DiffusionConnection
    std::vector< double >
      last_y_values; //!< remembers y_values from last wfr_update
    UniversalDataLogger< siegert_neuron >
      logger_; //!< Logger for all analog data
  };

  // ----------------------------------------------------------------

  /**
   * Internal variables of the model.
   */
  struct Variables_
  {

    // propagators
    double P1_;
    double P2_;
  };

  //! Read out the rate
  double
  get_rate_() const
  {
    return S_.r_;
  }

  // ----------------------------------------------------------------

  Parameters_ P_;
  State_ S_;
  Variables_ V_;
  Buffers_ B_;

  //! Mapping of recordables names to access functions
  static RecordablesMap< siegert_neuron > recordablesMap_;
};

inline void
siegert_neuron::update( Time const& origin, const long from, const long to )
{
  update_( origin, from, to, false );
}

inline bool
siegert_neuron::wfr_update( Time const& origin, const long from, const long to )
{
  State_ old_state = S_; // save state before wfr update
  const bool wfr_tol_exceeded = update_( origin, from, to, true );
  S_ = old_state; // restore old state

  return not wfr_tol_exceeded;
}

inline port
siegert_neuron::handles_test_event( DiffusionConnectionEvent&,
  rport receptor_type )
{
  if ( receptor_type == 0 )
    return 0;
  else if ( receptor_type == 1 )
    return 1;
  else
    throw UnknownReceptorType( receptor_type, get_name() );
}

inline port
siegert_neuron::handles_test_event( DataLoggingRequest& dlr,
  rport receptor_type )
{
  if ( receptor_type != 0 )
    throw UnknownReceptorType( receptor_type, get_name() );
  return B_.logger_.connect_logging_device( dlr, recordablesMap_ );
}

inline void
siegert_neuron::get_status( DictionaryDatum& d ) const
{
  P_.get( d );
  S_.get( d );
  Archiving_Node::get_status( d );
  ( *d )[ names::recordables ] = recordablesMap_.get_list();
}

inline void
siegert_neuron::set_status( const DictionaryDatum& d )
{
  Parameters_ ptmp = P_; // temporary copy in case of errors
  ptmp.set( d );         // throws if BadProperty
  State_ stmp = S_;      // temporary copy in case of errors
  stmp.set( d );         // throws if BadProperty

  // We now know that (ptmp, stmp) are consistent. We do not
  // write them back to (P_, S_) before we are also sure that
  // the properties to be set in the parent class are internally
  // consistent.
  Archiving_Node::set_status( d );

  // if we get here, temporaries contain consistent set of properties
  P_ = ptmp;
  S_ = stmp;
}

} // namespace

#endif // HAVE_GSL
#endif // SIEGERT_NEURON_H

max_stars_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster

max_issues_repo_name: zlchai/SNN-simulator-on-PYNQcluster

max_forks_repo_name: OpenHEC/SNN-simulator-on-PYNQcluster


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>
#include <petscsnes.h>

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/vecimpl.h"
     #include "petsc/private/matimpl.h"
     #include "petsc/private/pcimpl.h"
     #include "petsc/private/kspimpl.h"
     #include "petsc/private/snesimpl.h"
  #else
     #include "petsc-private/vecimpl.h"
     #include "petsc-private/matimpl.h"
     #include "petsc-private/pcimpl.h"
     #include "petsc-private/kspimpl.h"
     #include "petsc-private/snesimpl.h"
  #endif
#else
  #include "private/vecimpl.h"
  #include "private/matimpl.h"
  #include "private/pcimpl.h"
  #include "private/kspimpl.h"
  #include "private/snesimpl.h"
#endif

#include "common-driver-utils.h"

#define BSSCR_report_op(op,v,i,func_name) \
  if((op)) { PetscViewerASCIIPrintf((v), "[%.2d] %s: valid operation \n", i, func_name ); } \
  else { PetscViewerASCIIPrintf((v),     "[%.2d] %s: OPERATION NOT DEFINED \n", i, func_name ); } \

PetscErrorCode BSSCR_MatListOperations( Mat A, PetscViewer v )
{
	MatOps ops = A->ops;
	MatType type;
	PetscInt i;

	MatGetType( A, &type );
	PetscViewerASCIIPrintf(v, "Operations available for MatType: %s \n", type );

	PetscViewerASCIIPushTab( v );

	i=0;
	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->setvalues, v,  i++, "MatSetValues" );
	BSSCR_report_op( ops->getrow, v,     i++, "MatSetValues" );
	BSSCR_report_op( ops->restorerow, v, i++, "MatRestoreRow" );
	BSSCR_report_op( ops->mult, v,       i++, "MatMult" );
	BSSCR_report_op( ops->multadd, v,    i++, "MatMultAdd" );

	BSSCR_report_op( ops->multtranspose, v,    i++, "MatMultTranspose" );
	BSSCR_report_op( ops->multtransposeadd, v, i++, "MatMultTransposeAdd" );
	BSSCR_report_op( ops->solve, v,            i++, "MatSolve" );
	BSSCR_report_op( ops->solveadd, v,         i++, "MatSolveAdd" );
	BSSCR_report_op( ops->solvetranspose, v,   i++, "MatSolveTranspose" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->solvetransposeadd, v,  i++, "MatSolveTransposeAdd" );
	BSSCR_report_op( ops->lufactor, v,           i++, "MatLUFactor" );
	BSSCR_report_op( ops->choleskyfactor, v,     i++, "MatCholeskyFactor" );
#if( ( PETSC_VERSION_MAJOR >= 3 ) && ( PETSC_VERSION_MINOR > 0 ) )
	BSSCR_report_op( ops->sor, v,                i++, "MatSor" );
#else
	BSSCR_report_op( ops->relax, v,              i++, "MatRelax" );
#endif
	BSSCR_report_op( ops->transpose, v,          i++, "MatTranspose" );

	BSSCR_report_op( ops->getinfo, v,       i++, "MatGetInfo" );
	BSSCR_report_op( ops->equal, v,         i++, "MatEqual" );
	BSSCR_report_op( ops->getdiagonal, v,   i++, "MatGetDiagonal" );
	BSSCR_report_op( ops->diagonalscale, v, i++, "MatDiagonalScale" );
	BSSCR_report_op( ops->norm, v,          i++, "MatNorm" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->assemblybegin, v, i++, "MatAssemblyBegin" );
	BSSCR_report_op( ops->assemblyend, v,   i++, "MatAssemblyEnd" );
#if( ( ( PETSC_VERSION_MAJOR == 3 ) && ( PETSC_VERSION_MINOR < 1 ) ) || ( PETSC_VERSION_MAJOR <= 2 ) )
	BSSCR_report_op( ops->compress, v,      i++, "MatCompress" );
#endif
	BSSCR_report_op( ops->setoption, v,     i++, "MatSetOption" );
	BSSCR_report_op( ops->zeroentries, v,   i++, "MatZeroEntries" );

	BSSCR_report_op( ops->zerorows, v,               i++, "Stg_MatZeroRows" );
	BSSCR_report_op( ops->lufactorsymbolic, v,       i++, "MatLUFactorSymbolic" );
	BSSCR_report_op( ops->lufactornumeric, v,        i++, "MatLUFactorNumeric" );
	BSSCR_report_op( ops->choleskyfactorsymbolic, v, i++, "MatCholeskyFactorSymbolic" );
	BSSCR_report_op( ops->choleskyfactornumeric, v,  i++, "MatCholeskyFactorNumeric" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");

	BSSCR_report_op( ops->ilufactorsymbolic, v,  i++, "MatILUFactorSymbolic" );
	BSSCR_report_op( ops->iccfactorsymbolic, v,  i++, "MatICCFactorSymbolic" );

	BSSCR_report_op( ops->duplicate, v,     i++, "MatDuplicate" );
	BSSCR_report_op( ops->forwardsolve, v,  i++, "MatForwardSolve" );
	BSSCR_report_op( ops->backwardsolve, v, i++, "MatBackwardSolve" );
	BSSCR_report_op( ops->ilufactor, v,     i++, "MatILUFactor" );
	BSSCR_report_op( ops->iccfactor, v,     i++, "MatICCFactor" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->axpy, v,            i++, "MatAXPY" );

#if( (PETSC_VERSION_MAJOR == 3) &&  (PETSC_VERSION_MINOR > 7) )
    BSSCR_report_op( ops->createsubmatrices, v,  i++, "MatCreateSubMatrices" );
#else
    BSSCR_report_op( ops->getsubmatrices, v,  i++, "MatGetSubMatrices" );
#endif

    BSSCR_report_op( ops->increaseoverlap, v, i++, "MatIncreaseOverlap" );
	BSSCR_report_op( ops->getvalues, v,       i++, "MatGetValues" );
	BSSCR_report_op( ops->copy, v,            i++, "MatCopy" );

	BSSCR_report_op( ops->getrowmax, v,   i++, "MatGetRowMax" );
	BSSCR_report_op( ops->scale, v,       i++, "MatScale" );
	BSSCR_report_op( ops->shift, v,       i++, "MatShift" );
	BSSCR_report_op( ops->diagonalset, v, i++, "MatDiagonalSet" );
#if( ( ( PETSC_VERSION_MAJOR == 3 ) && ( PETSC_VERSION_MINOR < 1 ) ) || ( PETSC_VERSION_MAJOR <= 2 ) )
	BSSCR_report_op( ops->iludtfactor, v, i++, "MatILUDTFactor" );
#endif

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");

	BSSCR_report_op( ops->getrowij, v,        i++, "MatGetRowIJ" );
	BSSCR_report_op( ops->restorerowij, v,    i++, "MatRestoreRowIJ" );
	BSSCR_report_op( ops->getcolumnij, v,     i++, "MatGetColumnIJ" );
	BSSCR_report_op( ops->restorecolumnij, v, i++, "MatRestoreColumnIJ" );

	BSSCR_report_op( ops->fdcoloringcreate, v, i++, "MatFDColoringCreate" );
	BSSCR_report_op( ops->coloringpatch, v,    i++, "MatColoringPatch" );
	BSSCR_report_op( ops->setunfactored, v,    i++, "MatSetUnFactored" );
	BSSCR_report_op( ops->permute, v,          i++, "MatPermute" );
	BSSCR_report_op( ops->setvaluesblocked, v, i++, "MatSetValuesBlocked" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");

#if( (PETSC_VERSION_MAJOR == 3) &&  (PETSC_VERSION_MINOR > 7) )
    BSSCR_report_op( ops->createsubmatrix, v,  i++, "MatCreateSubMatrix" );
#else
    BSSCR_report_op( ops->getsubmatrix, v,  i++, "MatGetSubMatrix" );
#endif

	BSSCR_report_op( ops->destroy, v,       i++, "MatDestroy" );
	BSSCR_report_op( ops->view, v,          i++, "MatView" );
	BSSCR_report_op( ops->convertfrom, v,   i++, "MatConvertFrom" );

	BSSCR_report_op( ops->setlocaltoglobalmapping, v, i++, "MatSetLocalToGlobalMapping" );
	BSSCR_report_op( ops->setvalueslocal, v,          i++, "MatSetValuesLocal" );
	BSSCR_report_op( ops->zerorowslocal, v,           i++, "Stg_MatZeroRowsLocal" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->getrowmaxabs, v,    i++, "MatGetRowMaxAbs" );
	BSSCR_report_op( ops->convert, v,         i++, "MatConvert" );
//    BSSCR_report_op( ops->setcoloring, v,     i++, "MatSetColoring" );  NOT SURE IF THIS IS AVAIL IN 3.8

	BSSCR_report_op( ops->setvaluesadifor, v, i++, "MatSetValuesAdicfor" );

	BSSCR_report_op( ops->fdcoloringapply, v,              i++, "MatFDCloringApply" );
	BSSCR_report_op( ops->setfromoptions, v,               i++, "MatSetFromOptions" );
	BSSCR_report_op( ops->multconstrained, v,              i++, "MatMultConstrained" );
	BSSCR_report_op( ops->multtransposeconstrained, v,     i++, "MatMultTransposeConstrained" );
	/* data member depreciated from PETSc v. 3.0.0 */
#if( PETSC_VERSION_MAJOR <= 2 )
	BSSCR_report_op( ops->ilufactorsymbolicconstrained, v, i++, "MatILUFactorSymbolicConstrained" );
#endif
	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
#if( ( ( PETSC_VERSION_MAJOR == 3 ) && ( PETSC_VERSION_MINOR < 2 ) ) )
	BSSCR_report_op( ops->permutesparsify, v, i++, "MatPermuteSparsify" );
#endif
	BSSCR_report_op( ops->mults, v,           i++, "MatMults" );
	BSSCR_report_op( ops->solves, v,          i++, "MatSolves" );
	BSSCR_report_op( ops->getinertia, v,      i++, "MatGetInertia" );
	BSSCR_report_op( ops->load, v,            i++, "MatLoad" );

	BSSCR_report_op( ops->issymmetric, v,             i++, "MatIsSymmetric" );
	BSSCR_report_op( ops->ishermitian, v,             i++, "MatIsHermitian" );
	BSSCR_report_op( ops->isstructurallysymmetric, v, i++, "MatIsStructurallySymmetric" );
#if( ( ( PETSC_VERSION_MAJOR == 3 ) && ( PETSC_VERSION_MINOR < 1 ) ) || ( PETSC_VERSION_MAJOR <= 2 ) )
	BSSCR_report_op( ops->pbrelax, v,                 i++, "MatPBRelax" );
#endif
	BSSCR_report_op( ops->getvecs, v,                 i++, "MatGetVec" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->matmult, v,         i++, "MatMatMult" );
	BSSCR_report_op( ops->matmultsymbolic, v, i++, "MatMatMultSymbolic" );
	BSSCR_report_op( ops->matmultnumeric, v,  i++, "MatMatMultNumeric" );
	BSSCR_report_op( ops->ptap, v,            i++, "MatPtAP" );
	BSSCR_report_op( ops->ptapsymbolic, v,    i++, "MatPtAPSymbolic" );

	BSSCR_report_op( ops->ptapnumeric, v,              i++, "MatPtAPNumeric" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");

	BSSCR_report_op( ops->conjugate, v,           i++, "MatConjugate" );
	//BSSCR_report_op( ops->setsizes, v,            i++, "MatSetSizes" );

	BSSCR_report_op( ops->setvaluesrow, v,              i++, "MatSetValuesRow" );
	BSSCR_report_op( ops->realpart, v,                  i++, "MatRealPart" );
	BSSCR_report_op( ops->imaginarypart, v,             i++, "MatImaginaryPart" );
	BSSCR_report_op( ops->getrowuppertriangular, v,     i++, "MatGetRowUpperTriangular" );
	BSSCR_report_op( ops->restorerowuppertriangular, v, i++, "MatRestoreRowUpperTriangular" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->matsolve, v,           i++, "MatMatSolve" );
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR <6) )
	BSSCR_report_op( ops->getredundantmatrix, v, i++, "MatGetRedundantMatrix" );
#endif
	BSSCR_report_op( ops->getrowmin, v,          i++, "MatGetRowMin" );
	BSSCR_report_op( ops->getcolumnvector, v,    i++, "MatGetColumnVector" );

	PetscViewerASCIIPopTab( v );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_VecListOperations( Vec x, PetscViewer v )
{
	VecOps ops = x->ops;
	VecType type;
	PetscInt i;

	VecGetType( x, &type );
	PetscViewerASCIIPrintf(v, "Operations available for VecType: %s \n", type );

	PetscViewerASCIIPushTab( v );

	PetscViewerASCIIPopTab( v );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_KSPListOperations( KSP ksp, PetscViewer v )
{
	KSPOps ops = ksp->ops;
	KSPType type;
	PetscInt i;

	KSPGetType( ksp, &type );
	PetscViewerASCIIPrintf(v, "Operations available for KSPType: %s \n", type );

	PetscViewerASCIIPushTab( v );

	i=0;
	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->buildsolution, v,  i++, "KSPBuildSolution" );
	BSSCR_report_op( ops->buildresidual, v,  i++, "KSPBuildResidual" );
	BSSCR_report_op( ops->solve, v,          i++, "KSPSolve" );
	BSSCR_report_op( ops->setup, v,          i++, "KSPSetUp" );
	BSSCR_report_op( ops->setfromoptions, v, i++, "KSPSetFromOptions" );

	BSSCR_report_op( ops->publishoptions, v,               i++, "KSPPublishOptions" );
	BSSCR_report_op( ops->computeextremesingularvalues, v, i++, "KSPComputeExtremeSingularValues" );
	BSSCR_report_op( ops->computeeigenvalues, v,           i++, "KSPComputeEigenvalues" );
	BSSCR_report_op( ops->destroy, v, i++, "BSSCR_KSPDestroy" );
	BSSCR_report_op( ops->view, v,    i++, "KSPView" );

	PetscViewerASCIIPopTab( v );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_PCListOperations( PC pc, PetscViewer v )
{
	PCOps ops = pc->ops;
	PCType type;
	PetscInt i;

	PCGetType( pc, &type );
	PetscViewerASCIIPrintf(v, "Operations available for PCType: %s \n", type );

	PetscViewerASCIIPushTab( v );

	i=0;
	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->setup, v,           i++, "PCSetUp" );
	BSSCR_report_op( ops->apply, v,           i++, "PCApply" );
	BSSCR_report_op( ops->applyrichardson, v, i++, "PCApplyRichardson" );
	BSSCR_report_op( ops->applyBA, v,         i++, "PCApplyBA" );
	BSSCR_report_op( ops->applytranspose, v,  i++, "PCApplyTranspose" );

	BSSCR_report_op( ops->applyBAtranspose, v,  i++, "PCApplyBATranspose" );
	BSSCR_report_op( ops->setfromoptions, v,    i++, "PCSetFromOptions" );
	BSSCR_report_op( ops->presolve, v,          i++, "PCPreSolve" );
	BSSCR_report_op( ops->postsolve, v,         i++, "PCPostSolve" );
	BSSCR_report_op( ops->getfactoredmatrix, v, i++, "PCGetFactoredMatrix" );

	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");
	BSSCR_report_op( ops->applysymmetricleft, v,  i++, "PCApplySymmetricLeft" );
	BSSCR_report_op( ops->applysymmetricright, v, i++, "PCApplySymmetricRight" );
	BSSCR_report_op( ops->setuponblocks, v,       i++, "PCSetupOnBlocks" );
	BSSCR_report_op( ops->destroy, v,             i++, "PCDestroy" );
	BSSCR_report_op( ops->view, v,                i++, "PCView" );

	PetscViewerASCIIPopTab( v );

	PetscFunctionReturn(0);
}

PetscErrorCode BSSCR_SNESListOperations( SNES snes, PetscViewer v )
{
	SNESOps ops = snes->ops;
	SNESType type;
	PetscInt i;

	SNESGetType( snes, &type );
	PetscViewerASCIIPrintf(v, "Operations available for SNESType: %s \n", type );

	PetscViewerASCIIPushTab( v );

	i=0;
	PetscViewerASCIIPrintf(v, "------------------------------------------------\n");

	BSSCR_report_op( ops->computescaling, v,  i++, "SNESComputeScaling" );
	BSSCR_report_op( ops->update, v,          i++, "SNESUpdate" );
	BSSCR_report_op( ops->converged, v,       i++, "SNESConverged" );

	BSSCR_report_op( ops->setup, v,           i++, "SNESSetUp" );
	BSSCR_report_op( ops->solve, v,           i++, "SNESSolve" );
	BSSCR_report_op( ops->view, v,            i++, "SNESView" );
	BSSCR_report_op( ops->setfromoptions, v,  i++, "SNESSetFromOptions" );
	BSSCR_report_op( ops->destroy, v,         i++, "SNESDestroy" );

	PetscViewerASCIIPopTab( v );

	PetscFunctionReturn(0);
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <gsl/gsl_interp.h>
#include "field.h"

#define StartTest -4E-7
#define StopTest 4E-7

signal_type *
formatExpImpResp(int numPnts, double *timeValues, double *voltageValues,
	int samplingFrequencyHz, int verbose)
{
double *newTimeValues;
double *newVoltageValues;
double timeAtMaxIntensity;
double minTime = DBL_MIN;
double maxTime = -DBL_MIN;
double maxVoltage = -DBL_MIN;
double stepSize;
int i, j, maxVoltageIndex, numSteps;
int startIndex, stopIndex;
signal_type *impulseResponse;
int pntsSampled;

/* find the max voltage */

	for (i = 0; i < numPnts; i++) {
		if (maxVoltage < voltageValues[i]) {
			maxVoltage = voltageValues[i];
			maxVoltageIndex = i;
			}
		}

	if (verbose >= 1) fprintf(stderr, "index %d max voltage %e\n",
		maxVoltageIndex, maxVoltage);

/* normalize the voltages */

	for (i = 0; i < numPnts; i++)
		voltageValues[i] /= maxVoltage;

/*
 * center the time axis around the max intensity. the matlab code recalculated
 * the index of the maximum value, but the normalization shouldn't change
 * that, so I skipped it. note that I have to save the time associated with
 * the max intensity, because in the time array, it goes to 0 as I do the
 * centering.
 */

	timeAtMaxIntensity = timeValues[maxVoltageIndex];

	for (i = 0; i < numPnts; i++)
		timeValues[i] = timeValues[i] - timeAtMaxIntensity;

/* now find the min and max times */

	for (i = 0; i < numPnts; i++) {
		if (minTime > timeValues[i]) minTime = timeValues[i];

		if (maxTime < timeValues[i]) maxTime = timeValues[i];
		}

	if (verbose >= 1) fprintf(stderr, "min time %e max time %e\n", minTime,
		maxTime);

/* re-sample the data to match the Field II sampling frequency */

	if (verbose >= 1) fprintf(stderr, "sampling %d\n", samplingFrequencyHz);
	stepSize = 1.0 / samplingFrequencyHz;

	if (verbose >= 1) fprintf(stderr, "diff %e\n", maxTime - minTime);
	numSteps = (int )ceil(((maxTime - minTime) * samplingFrequencyHz));

	if (verbose >= 1) fprintf(stderr, "step size %e numSteps %d\n", stepSize, numSteps);

	if ((newTimeValues = (double *)malloc(sizeof(double) * numSteps)) == NULL) {
		fprintf(stderr, "in readExpData, couldn't allocate space for new times\n");
		return(0);
		}

	for (i = 0; i < numSteps; i++) {
		newTimeValues[i] = minTime + i * stepSize;
		}

/*
 * now we find the voltage values at each of the new times by interpolating
 * from the old times and voltages.
 */

/* initialize and allocate the gsl objects */

	gsl_interp *interpolation = gsl_interp_alloc (gsl_interp_linear, numPnts);

	gsl_interp_init(interpolation, timeValues, voltageValues, numPnts);

	gsl_interp_accel * accelerator =  gsl_interp_accel_alloc();

/* get interpolation the new times */

	if ((newVoltageValues = (double *)malloc(sizeof(double) * numSteps)) == NULL) {
		fprintf(stderr, "in readExpData, couldn't allocate space for new voltage\n");
		return(0);
		}

	for (i = 0; i < numSteps; i++) {
		newVoltageValues[i] = gsl_interp_eval(interpolation, timeValues,
			voltageValues, newTimeValues[i], accelerator);

		if (verbose >= 3) fprintf(stderr, "\n%e", newVoltageValues[i]);
		}

/* find the indices of NewTime to use for the Field II impulse response */

	for (i = 0; i < numSteps; i++) {
		if (newTimeValues[i] > StartTest) {
			startIndex = i;
			break;
			}
		}

	for (i = 0; i < numSteps; i++)
		if (newTimeValues[i] > StopTest) {
			stopIndex = i;
			break;
			}

	pntsSampled = stopIndex - startIndex;

	if (verbose >= 1) fprintf(stderr,
		"\nstartIndex %d, stopIndex %d pntsSampled %d\n", startIndex,
		stopIndex, pntsSampled);

	impulseResponse = alloc_signal(pntsSampled, 0);

	for (i = startIndex, j = 0; i < stopIndex; i++, j++) {
		impulseResponse->data[j] = newVoltageValues[i];
		if (verbose >= 3) fprintf(stderr, "%e\n", newVoltageValues[i]);
		}

	return(impulseResponse);
}

max_stars_repo_name: suyashkumar/fem

max_issues_repo_name: suyashkumar/fem

max_forks_repo_name: suyashkumar/fem


text:
// File: nlopt_wrapper.c
// Author: Jannis Harder
// Description: Wrapper for nlopt using context switching

#include <stdlib.h>
#include <ucontext.h>
#include <string.h>
#include <stdio.h> // For debugging TODO: remove
#include <nlopt.h>

#define STACKSIZE (8 * 1024 * 1024)

enum status {
	ST_IDLE = 0,
	ST_RUNNING,
	ST_VALUE,
	ST_GRAD,
	ST_DONE
};

struct function {
	struct wrapper * w;
	struct function * next;
	int id;
};

struct wrapper {
	struct function f;
	nlopt_opt opt;
	int n;
	ucontext_t julia_ctx;
	ucontext_t nlopt_ctx;
	char *nlopt_stack;
	double final_f;
	double *opt_x;
	double *eval_x;
	double *eval_grad;
	double eval_f;
	int force_stop;
	int result;
	int function_id;
	enum status status;
};

double nlopt_wrapper_f(unsigned n, const double *x, double *grad, void *f_data);
void nlopt_wrapper_optimize_thread(struct wrapper* w);

void nlopt_wrapper_version(int *version)
{
	nlopt_version(&version[0], &version[1], &version[2]);
	version[3] = 2;
	version[4] = 2;
	version[5] = 4;
	version[6] = 0;
}

struct wrapper *nlopt_wrapper_create(int type_id, int dimensions)
{
	nlopt_opt opt = nlopt_create(type_id, dimensions);
	if (!opt)
		return NULL;
	struct wrapper *w = calloc(1, sizeof(struct wrapper));
	w->n = dimensions;
	w->opt = opt;
	w->f.w = w;
	w->f.id = 1;
	return w;
}

void nlopt_wrapper_free(struct wrapper *w)
{
	nlopt_destroy(w->opt);
	free(w->opt_x);
	struct function *f = &w->f;
	while (f) {
		struct function *next = f->next;
		free(f);
		f = next;
	}
}

int nlopt_wrapper_objective(struct wrapper *w, int max)
{
	if (max)
		return nlopt_set_max_objective(w->opt, &nlopt_wrapper_f, &w->f);
	else
		return nlopt_set_min_objective(w->opt, &nlopt_wrapper_f, &w->f);
}

double nlopt_wrapper_f(unsigned n, const double *x, double *grad, void *f_data)
{
	struct function *f = f_data;
	struct wrapper *w = f->w;

	w->function_id = f->id;

	w->status = grad ? ST_GRAD : ST_VALUE;
	memcpy(w->eval_x, x, n * sizeof(double));
	swapcontext(&w->nlopt_ctx, &w->julia_ctx);
	if (w->force_stop)
		nlopt_force_stop(w->opt);
	if (grad)
		memcpy(grad, w->eval_grad, n * sizeof(double));
	return w->eval_f;
}

void nlopt_wrapper_optimize_start(struct wrapper *w, double *x)
{
	int i;
	getcontext(&w->nlopt_ctx);
	w->nlopt_stack = calloc(1, STACKSIZE);
	w->nlopt_ctx.uc_stack.ss_sp = w->nlopt_stack;
	w->nlopt_ctx.uc_stack.ss_size = STACKSIZE;
	w->nlopt_ctx.uc_link = &w->julia_ctx;
	makecontext(&w->nlopt_ctx, (void(*)())&nlopt_wrapper_optimize_thread, 1, w);

	w->opt_x = calloc(w->n, sizeof(double));
	memcpy(w->opt_x, x, w->n * sizeof(double));
}

int nlopt_wrapper_optimize_callback(struct wrapper *w, double *x, double *grad, double f, int force_stop, int *function_id)
{
	w->eval_f = f;
	w->eval_x = x;
	w->eval_grad = grad;
	w->force_stop = force_stop;
	swapcontext(&w->julia_ctx, &w->nlopt_ctx);
	*function_id = w->function_id;
	return w->status;
}

int nlopt_wrapper_optimize_finalize(struct wrapper *w, double *x, double *f)
{
	memcpy(x, w->opt_x, w->n * sizeof(double));
	free(w->opt_x);
	w->opt_x = NULL;
	*f = w->final_f;
	free(w->nlopt_stack);
	w->nlopt_stack = NULL;
	return w->result;
}

void nlopt_wrapper_optimize_thread(struct wrapper *w)
{
	w->result = nlopt_optimize(w->opt, w->opt_x, &w->final_f);
	w->status = ST_DONE;
}

void nlopt_wrapper_dimopt(struct wrapper *w, double *x, int i)
{
	switch (i) {
	case 0: nlopt_set_lower_bounds(w->opt, x); break;
	case 1: nlopt_set_upper_bounds(w->opt, x); break;
	case 2: nlopt_set_xtol_abs(w->opt, x); break;
	}
}

void nlopt_wrapper_doubleopt(struct wrapper *w, double v, int i)
{
	switch (i) {
	case 0: nlopt_set_stopval(w->opt, v); break;
	case 1: nlopt_set_ftol_rel(w->opt, v); break;
	case 2: nlopt_set_ftol_abs(w->opt, v); break;
	case 3: nlopt_set_xtol_rel(w->opt, v); break;
	case 4: nlopt_set_maxtime(w->opt, v); break;
	}
}

void nlopt_wrapper_intopt(struct wrapper *w, long int v, int i)
{
	switch (i) {
	case 0: nlopt_set_maxeval(w->opt, v); break;
	case 1: nlopt_set_population(w->opt, v < 0 ? 0 : v); break;
	case 2: nlopt_srand(v); break;
	}
}

void nlopt_wrapper_add_constraint(struct wrapper *w, int id, double tolerance, int equality)
{
	struct function *f = calloc(1, sizeof(struct function));

	f->w = w;
	f->next = w->f.next;
	f->id = id;
	w->f.next = f;
	if (equality)
		nlopt_add_equality_constraint(w->opt, &nlopt_wrapper_f, f, tolerance);
	else
		nlopt_add_inequality_constraint(w->opt, &nlopt_wrapper_f, f, tolerance);
}

void nlopt_wrapper_local_optimizer(struct wrapper *w, struct wrapper *local_w)
{
	nlopt_set_local_optimizer(w->opt, local_w->opt);
}

max_stars_repo_name: MetalNinjas/julia-nlopt

max_issues_repo_name: MetalNinjas/julia-nlopt

max_forks_repo_name: MetalNinjas/julia-nlopt