text:
/* linalg/householder.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004, 2007, 2010 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.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

#include <gsl/gsl_linalg.h>

double
gsl_linalg_householder_transform (gsl_vector * v)
{
  /* replace v[0:n-1] with a householder vector (v[0:n-1]) and
     coefficient tau that annihilate v[1:n-1] */

  const size_t n = v->size ;

  if (n == 1)
    {
      return 0.0; /* tau = 0 */
    }
  else
    {
      double alpha, beta, tau ;

      gsl_vector_view x = gsl_vector_subvector (v, 1, n - 1) ;

      double xnorm = gsl_blas_dnrm2 (&x.vector);

      if (xnorm == 0)
        {
          return 0.0; /* tau = 0 */
        }

      alpha = gsl_vector_get (v, 0) ;
      beta = - (alpha >= 0.0 ? +1.0 : -1.0) * hypot(alpha, xnorm) ;
      tau = (beta - alpha) / beta ;

      {
        double s = (alpha - beta);

        if (fabs(s) > GSL_DBL_MIN)
          {
            gsl_blas_dscal (1.0 / s, &x.vector);
            gsl_vector_set (v, 0, beta) ;
          }
        else
          {
            gsl_blas_dscal (GSL_DBL_EPSILON / s, &x.vector);
            gsl_blas_dscal (1.0 / GSL_DBL_EPSILON, &x.vector);
            gsl_vector_set (v, 0, beta) ;
          }
      }

      return tau;
    }
}

int
gsl_linalg_householder_hm (double tau, const gsl_vector * v, gsl_matrix * A)
{
  /* applies a householder transformation v,tau to matrix m */

  if (tau == 0.0)
    {
      return GSL_SUCCESS;
    }

#ifdef USE_BLAS
  {
    gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
    gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
    size_t j;

    for (j = 0; j < A->size2; j++)
      {
        double wj = 0.0;
        gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
        gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);
        wj += gsl_matrix_get(A,0,j);

        {
          double A0j = gsl_matrix_get (A, 0, j);
          gsl_matrix_set (A, 0, j, A0j - tau *  wj);
        }

        gsl_blas_daxpy (-tau * wj, &v1.vector, &A1j.vector);
      }
  }
#else
  {
    size_t i, j;

    for (j = 0; j < A->size2; j++)
      {
        /* Compute wj = Akj vk */

        double wj = gsl_matrix_get(A,0,j);

        for (i = 1; i < A->size1; i++)  /* note, computed for v(0) = 1 above */
          {
            wj += gsl_matrix_get(A,i,j) * gsl_vector_get(v,i);
          }

        /* Aij = Aij - tau vi wj */

        /* i = 0 */
        {
          double A0j = gsl_matrix_get (A, 0, j);
          gsl_matrix_set (A, 0, j, A0j - tau *  wj);
        }

        /* i = 1 .. M-1 */

        for (i = 1; i < A->size1; i++)
          {
            double Aij = gsl_matrix_get (A, i, j);
            double vi = gsl_vector_get (v, i);
            gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
          }
      }
  }
#endif

  return GSL_SUCCESS;
}

int
gsl_linalg_householder_mh (double tau, const gsl_vector * v, gsl_matrix * A)
{
  /* applies a householder transformation v,tau to matrix m from the
     right hand side in order to zero out rows */

  if (tau == 0)
    return GSL_SUCCESS;

  /* A = A - tau w v' */

#ifdef USE_BLAS
  {
    gsl_vector_const_view v1 = gsl_vector_const_subvector (v, 1, v->size - 1);
    gsl_matrix_view A1 = gsl_matrix_submatrix (A, 0, 1, A->size1, A->size2-1);
    size_t i;

    for (i = 0; i < A->size1; i++)
      {
        double wi = 0.0;
        gsl_vector_view A1i = gsl_matrix_row(&A1.matrix, i);
        gsl_blas_ddot (&A1i.vector, &v1.vector, &wi);
        wi += gsl_matrix_get(A,i,0);

        {
          double Ai0 = gsl_matrix_get (A, i, 0);
          gsl_matrix_set (A, i, 0, Ai0 - tau *  wi);
        }

        gsl_blas_daxpy(-tau * wi, &v1.vector, &A1i.vector);
      }
  }
#else
  {
    size_t i, j;

    for (i = 0; i < A->size1; i++)
      {
        double wi = gsl_matrix_get(A,i,0);

        for (j = 1; j < A->size2; j++)  /* note, computed for v(0) = 1 above */
          {
            wi += gsl_matrix_get(A,i,j) * gsl_vector_get(v,j);
          }

        /* j = 0 */

        {
          double Ai0 = gsl_matrix_get (A, i, 0);
          gsl_matrix_set (A, i, 0, Ai0 - tau *  wi);
        }

        /* j = 1 .. N-1 */

        for (j = 1; j < A->size2; j++)
          {
            double vj = gsl_vector_get (v, j);
            double Aij = gsl_matrix_get (A, i, j);
            gsl_matrix_set (A, i, j, Aij - tau * wi * vj);
          }
      }
  }
#endif

  return GSL_SUCCESS;
}

int
gsl_linalg_householder_hv (double tau, const gsl_vector * v, gsl_vector * w)
{
  /* applies a householder transformation v to vector w */
  const size_t N = v->size;

  if (tau == 0)
    return GSL_SUCCESS ;

  {
    /* compute d = v'w */

    double d0 = gsl_vector_get(w,0);
    double d1, d;

    gsl_vector_const_view v1 = gsl_vector_const_subvector(v, 1, N-1);
    gsl_vector_view w1 = gsl_vector_subvector(w, 1, N-1);

    gsl_blas_ddot (&v1.vector, &w1.vector, &d1);

    d = d0 + d1;

    /* compute w = w - tau (v) (v'w) */

    {
      double w0 = gsl_vector_get (w,0);
      gsl_vector_set (w, 0, w0 - tau * d);
    }

    gsl_blas_daxpy (-tau * d, &v1.vector, &w1.vector);
  }

  return GSL_SUCCESS;
}

int
gsl_linalg_householder_hm1 (double tau, gsl_matrix * A)
{
  /* applies a householder transformation v,tau to a matrix being
     build up from the identity matrix, using the first column of A as
     a householder vector */

  if (tau == 0)
    {
      size_t i,j;

      gsl_matrix_set (A, 0, 0, 1.0);

      for (j = 1; j < A->size2; j++)
        {
          gsl_matrix_set (A, 0, j, 0.0);
        }

      for (i = 1; i < A->size1; i++)
        {
          gsl_matrix_set (A, i, 0, 0.0);
        }

      return GSL_SUCCESS;
    }

  /* w = A' v */

#ifdef USE_BLAS
  {
    gsl_matrix_view A1 = gsl_matrix_submatrix (A, 1, 0, A->size1 - 1, A->size2);
    gsl_vector_view v1 = gsl_matrix_column (&A1.matrix, 0);
    size_t j;

    for (j = 1; j < A->size2; j++)
      {
        double wj = 0.0;   /* A0j * v0 */

        gsl_vector_view A1j = gsl_matrix_column(&A1.matrix, j);
        gsl_blas_ddot (&A1j.vector, &v1.vector, &wj);

        /* A = A - tau v w' */

        gsl_matrix_set (A, 0, j, - tau *  wj);

        gsl_blas_daxpy(-tau*wj, &v1.vector, &A1j.vector);
      }

    gsl_blas_dscal(-tau, &v1.vector);

    gsl_matrix_set (A, 0, 0, 1.0 - tau);
  }
#else
  {
    size_t i, j;

    for (j = 1; j < A->size2; j++)
      {
        double wj = 0.0;   /* A0j * v0 */

        for (i = 1; i < A->size1; i++)
          {
            double vi = gsl_matrix_get(A, i, 0);
            wj += gsl_matrix_get(A,i,j) * vi;
          }

        /* A = A - tau v w' */

        gsl_matrix_set (A, 0, j, - tau *  wj);

        for (i = 1; i < A->size1; i++)
          {
            double vi = gsl_matrix_get (A, i, 0);
            double Aij = gsl_matrix_get (A, i, j);
            gsl_matrix_set (A, i, j, Aij - tau * vi * wj);
          }
      }

    for (i = 1; i < A->size1; i++)
      {
        double vi = gsl_matrix_get(A, i, 0);
        gsl_matrix_set(A, i, 0, -tau * vi);
      }

    gsl_matrix_set (A, 0, 0, 1.0 - tau);
  }
#endif

  return GSL_SUCCESS;
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2017 Couchbase, Inc
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */

#pragma once

#include "monotonic.h"

#include <memcached/dockey.h>
#include <memcached/types.h>
#include <nlohmann/json_fwd.hpp>
#include <platform/sized_buffer.h>
#include <gsl/gsl>

#include <unordered_map>
#include <vector>

namespace Collections {

// The reserved name of the system owned, default collection.
const char* const _DefaultCollectionIdentifier = "_default";
static cb::const_char_buffer DefaultCollectionIdentifier(
        _DefaultCollectionIdentifier);

const char* const _DefaultScopeIdentifier = "_default";
static cb::const_char_buffer DefaultScopeIdentifier(_DefaultScopeIdentifier);

// SystemEvent keys or parts which will be made into keys
const char* const SystemSeparator = ":"; // Note this never changes
const char* const CollectionEventPrefixWithSeparator = "_collection:";
const char* const ScopeEventPrefixWithSeparator = "_scope:";

// Couchstore private file name for manifest data
const char CouchstoreManifest[] = "_local/collections_manifest";

// Length of the string excluding the zero terminator (i.e. strlen)
const size_t CouchstoreManifestLen = sizeof(CouchstoreManifest) - 1;

using ManifestUid = WeaklyMonotonic<uint64_t>;

// Map used in summary stats
using Summary = std::unordered_map<CollectionID, uint64_t>;

struct ManifestUidNetworkOrder {
    ManifestUidNetworkOrder(ManifestUid uid) : uid(htonll(uid)) {
    }
    ManifestUid to_host() const {
        return ntohll(uid);
    }
    ManifestUid::value_type uid;
};
static_assert(sizeof(ManifestUidNetworkOrder) == 8,
              "ManifestUidNetworkOrder must have fixed size of 8 bytes as "
              "written to disk.");

/**
 * Return a ManifestUid from a C-string.
 * A valid ManifestUid is a C-string where each character satisfies
 * std::isxdigit and can be converted to a ManifestUid by std::strtoull.
 *
 * @param uid C-string uid
 * @param len a length for validation purposes
 * @throws std::invalid_argument if uid is invalid
 */
ManifestUid makeUid(const char* uid, size_t len = 16);

/**
 * Return a ManifestUid from a std::string
 * A valid ManifestUid is a std::string where each character satisfies
 * std::isxdigit and can be converted to a ManifestUid by std::strtoull.
 *
 * @param uid std::string
 * @throws std::invalid_argument if uid is invalid
 */
static inline ManifestUid makeUid(const std::string& uid) {
    return makeUid(uid.c_str());
}

/**
 * Return a CollectionID from a C-string.
 * A valid CollectionID is a C-string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul.
 *
 * @param uid C-string uid
 * @throws std::invalid_argument if uid is invalid
 */
static inline CollectionID makeCollectionID(const char* uid) {
    // CollectionID is 8 characters max and smaller than a ManifestUid
    return gsl::narrow_cast<CollectionID>(makeUid(uid, 8));
}

/**
 * Return a CollectionID from a std::string
 * A valid CollectionID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 *
 * @param uid std::string
 * @throws std::invalid_argument if uid is invalid
 */
static inline CollectionID makeCollectionID(const std::string& uid) {
    return makeCollectionID(uid.c_str());
}

/**
 * Return a ScopeID from a C-string.
 * A valid CollectionID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 * @param uid C-string uid
 * @return std::invalid_argument if uid is invalid
 */
static inline ScopeID makeScopeID(const char* uid) {
    // ScopeId is 8 characters max and smaller than a ManifestUid
    return gsl::narrow_cast<ScopeID>(makeUid(uid, 8));
}

/**
 * Return a ScopeID from a std::string
 * A valid ScopeID is a std::string where each character satisfies
 * std::isxdigit and can be converted to a CollectionID by std::strtoul
 * @param uid std::string
 * @return std::invalid_argument if uid is invalid
 */
static inline ScopeID makeScopeID(const std::string& uid) {
    return makeScopeID(uid.c_str());
}

/**
 * The metadata of a single collection
 *
 * Default construction yields the default collection
 */
struct CollectionMetaData {
    ScopeID sid{ScopeID::Default}; // The scope that the collection belongs to
    CollectionID cid{CollectionID::Default}; // The collection's ID
    std::string name{_DefaultCollectionIdentifier}; // The collection's name
    cb::ExpiryLimit maxTtl; // The collection's maxTTL

    bool operator==(const CollectionMetaData& other) const {
        return sid == other.sid && cid == other.cid && name == other.name &&
               maxTtl == other.maxTtl;
    }
};

/**
 * All of the data a system event needs
 */
struct CreateEventData {
    ManifestUid manifestUid; // The Manifest which generated the event
    CollectionMetaData metaData; // The data of the new collection
};

struct DropEventData {
    ManifestUid manifestUid; // The Manifest which generated the event
    ScopeID sid; // The scope that the collection belonged to
    CollectionID cid; // The collection the event belongs to
};

struct CreateScopeEventData {
    ManifestUid manifestUid; // The Manifest which generated the event
    ScopeID sid; // The scope id
    std::string name; // The scope name
};

struct DropScopeEventData {
    ManifestUid manifestUid; // The Manifest which generated the event
    ScopeID sid; // The scope the event belongs to
};

/**
 * All of the data a DCP create event message will transmit in the value of the
 * message. This is the layout to be used on the wire and is in the correct
 * byte order
 */
struct CreateEventDcpData {
    CreateEventDcpData(const CreateEventData& ev)
        : manifestUid(ev.manifestUid),
          sid(ev.metaData.sid),
          cid(ev.metaData.cid) {
    }
    /// The manifest uid stored in network byte order ready for sending
    ManifestUidNetworkOrder manifestUid;
    /// The scope id stored in network byte order ready for sending
    ScopeIDNetworkOrder sid;
    /// The collection id stored in network byte order ready for sending
    CollectionIDNetworkOrder cid;
    // The size is sizeof(manifestUid) + sizeof(cid) + sizeof(sid)
    // (msvc won't allow that expression)
    constexpr static size_t size{16};
};

/**
 * All of the data a DCP create event message will transmit in the value of a
 * DCP system event message (when the collection is created with a TTL). This is
 * the layout to be used on the wire and is in the correct byte order
 */
struct CreateWithMaxTtlEventDcpData {
    CreateWithMaxTtlEventDcpData(const CreateEventData& ev)
        : manifestUid(ev.manifestUid),
          sid(ev.metaData.sid),
          cid(ev.metaData.cid),
          maxTtl(htonl(gsl::narrow_cast<uint32_t>(
                  ev.metaData.maxTtl.get().count()))) {
    }
    /// The manifest uid stored in network byte order ready for sending
    ManifestUidNetworkOrder manifestUid;
    /// The scope id stored in network byte order ready for sending
    ScopeIDNetworkOrder sid;
    /// The collection id stored in network byte order ready for sending
    CollectionIDNetworkOrder cid;
    /// The collection's maxTTL value (in network byte order)
    uint32_t maxTtl;
    // The size is sizeof(manifestUid) + sizeof(cid) + sizeof(sid) +
    //             sizeof(maxTTL) (msvc won't allow that expression)
    constexpr static size_t size{20};
};

/**
 * All of the data a DCP drop event message will transmit in the value of the
 * message. This is the layout to be used on the wire and is in the correct
 * byte order
 */
struct DropEventDcpData {
    DropEventDcpData(const DropEventData& data)
        : manifestUid(data.manifestUid), sid(data.sid), cid(data.cid) {
    }

    /// The manifest uid stored in network byte order ready for sending
    ManifestUidNetworkOrder manifestUid;
    /// The scope id stored in network byte order ready for sending
    ScopeIDNetworkOrder sid;
    /// The collection id stored in network byte order ready for sending
    CollectionIDNetworkOrder cid;
    // The size is sizeof(manifestUid) + sizeof(cid) (msvc won't allow that
    // expression)
    constexpr static size_t size{16};
};

/**
 * All of the data a DCP create scope event message will transmit in the value
 * of the message. This is the layout to be used on the wire and is in the
 * correct byte order
 */
struct CreateScopeEventDcpData {
    CreateScopeEventDcpData(const CreateScopeEventData& data)
        : manifestUid(data.manifestUid), sid(data.sid) {
    }
    /// The manifest uid stored in network byte order ready for sending
    ManifestUidNetworkOrder manifestUid;
    /// The scope id stored in network byte order ready for sending
    ScopeIDNetworkOrder sid;
    constexpr static size_t size{12};
};

/**
 * All of the data a DCP drop scope event message will transmit in the value of
 * the message. This is the layout to be used on the wire and is in the correct
 * byte order
 */
struct DropScopeEventDcpData {
    DropScopeEventDcpData(const DropScopeEventData& data)
        : manifestUid(data.manifestUid), sid(data.sid) {
    }

    /// The manifest uid stored in network byte order ready for sending
    ManifestUidNetworkOrder manifestUid;
    /// The collection id stored in network byte order ready for sending
    ScopeIDNetworkOrder sid;
    constexpr static size_t size{12};
};

/**
 * All unknown_collection errors should be accompanied by a error context
 * value which includes the manifest ID which was what the collection lookup
 * failed against.
 * @return json for use in setErrorJsonExtras
 */
nlohmann::json getUnknownCollectionErrorContext(uint64_t manifestUid);

/**
 * For creation of collection SystemEvents - The SystemEventFactory
 * glues the CollectionID into the event key (so create of x doesn't
 * collide with create of y). This method yields the 'keyExtra' parameter
 *
 * @param collection The value to turn into a string
 * @return the keyExtra parameter to be passed to SystemEventFactory
 */
std::string makeCollectionIdIntoString(CollectionID collection);

/**
 * For creation of scope SystemEvents - The SystemEventFactory
 * glues the ScopeID into the event key (so create of x doesn't
 * collide with create of y). This method yields the 'keyExtra' parameter
 *
 * @param sid The ScopeId to turn into a string
 * @return the keyExtra parameter to be passed to SystemEventFactory
 */
std::string makeScopeIdIntoString(ScopeID sid);

/**
 * For creation of collection SystemEvents - The SystemEventFactory
 * glues the CollectionID into the event key (so create of x doesn't
 * collide with create of y). This method basically reverses
 * makeCollectionIdIntoString so we can get a CollectionID from a
 * SystemEvent key
 *
 * @param key DocKey from a SystemEvent
 * @param separator the separator between the SystemEvent prefix and the
 *        CollectionID
 * @return the ID which was in the event
 */
CollectionID getCollectionIDFromKey(
        const DocKey& key,
        const char* separator = Collections::SystemSeparator);

/// Same as getCollectionIDFromKey but for events changing scopes
ScopeID getScopeIDFromKey(const DocKey& key,
                          const char* separator = Collections::SystemSeparator);

/**
 * Callback function for processing against dropped collections in an ephemeral
 * vb, returns true if the key at seqno should be dropped
 */
using IsDroppedEphemeralCb = std::function<bool(const DocKey&, int64_t)>;

} // end namespace Collections

ext: h


text:
#include <math.h>
#include <string.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_sf_hyperg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>

#include <fastpm/libfastpm.h>
#include <fastpm/prof.h>
#include <fastpm/lightcone.h>
#include <fastpm/logging.h>
#include <gsl/gsl_linalg.h>

#include "pmpfft.h"
#include "pmghosts.h"

static void
_matrix_invert(double (*matrix)[4],
                    double (*matrix_inv)[4],
                    int matrix_size)
                                                                        //square matrix..
                                                                        //matrix_size fixed to 4
{
    double temp[16];

    memcpy(temp, matrix, sizeof(double) * 16);

    gsl_matrix_view m = gsl_matrix_view_array(temp, matrix_size, matrix_size);
    gsl_matrix_view inv = gsl_matrix_view_array(&matrix_inv[0][0],matrix_size, matrix_size);
    gsl_permutation * p = gsl_permutation_alloc (matrix_size);

    int signum;
    gsl_linalg_LU_decomp (&m.matrix, p, &signum);
    gsl_linalg_LU_invert (&m.matrix, p, &inv.matrix);

    gsl_permutation_free (p);
}

void
fastpm_lc_init(FastPMLightCone * lc)
{
    gsl_set_error_handler_off(); // Turn off GSL error handler

    /* Allocation */
    lc->horizon = malloc(sizeof(FastPMHorizon));
    fastpm_horizon_init(lc->horizon, lc->cosmology);
    _matrix_invert(lc->glmatrix,lc->glmatrix_inv,4);

}

void
fastpm_lc_destroy(FastPMLightCone * lc)
{
    fastpm_horizon_destroy(lc->horizon);
    free(lc->horizon);
}

void
fastpm_usmesh_init(FastPMUSMesh * mesh, FastPMLightCone * lc,
            double target_volume,
            FastPMStore * source,
            size_t np_upper,
            double (*tileshifts)[3],
            int ntiles,
            double amin,
            double amax)
{

    mesh->target_volume = target_volume;
    mesh->source = source;
    mesh->amin = amin;
    mesh->amax = amax;
    mesh->lc = lc;
    mesh->tileshifts = malloc(sizeof(tileshifts[0]) * ntiles);
    mesh->ntiles = ntiles;
    mesh->np_before = 0;

    memcpy(mesh->tileshifts, tileshifts, sizeof(tileshifts[0]) * ntiles);

    mesh->event_handlers = NULL;
    mesh->p = malloc(sizeof(FastPMStore));
    /* for saving the density with particles */
    fastpm_store_init(mesh->p, source->name, np_upper,
                  COLUMN_ID | COLUMN_POS | COLUMN_VEL | COLUMN_MASK
                | COLUMN_AEMIT,
                FASTPM_MEMORY_HEAP
    );

    mesh->p->meta.M0 = source->meta.M0;       // FIXME: change this for ncdm mass defn?
}

void fastpm_usmesh_destroy(FastPMUSMesh * mesh)
{
    fastpm_destroy_event_handlers(&mesh->event_handlers);
    fastpm_store_destroy(mesh->p);
    free(mesh->tileshifts);
    free(mesh->p);
}

struct funct_params {
    FastPMLightCone *lc;
    FastPMStore * p;
    FastPMDriftFactor * drift;
    ptrdiff_t i;
    double tileshift[4];
    double a1;
    double a2;
    void * context;
};

static double
funct(double a, void *params)
{
    struct funct_params *Fp = (struct funct_params *) params;

    FastPMLightCone *lc = Fp->lc;
    FastPMStore * p = Fp->p;
    FastPMDriftFactor * drift = Fp->drift;
    ptrdiff_t i = Fp->i;
    int d;
    double xi[4];
    double xo[4];

    xi[3] = 1;
    if(p->v) {
        fastpm_drift_one(drift, p, i, xi, a);
    } else {
        for(d = 0; d < 3; d ++) {
            xi[d] = p->x[i][d];
        }
    }
    for(d = 0; d < 4; d ++) {
        xi[d] += Fp->tileshift[d];
    }
    /* transform the coordinate */
    fastpm_gldot(lc->glmatrix, xi, xo);

    double distance = fastpm_lc_distance(lc, xo);

    /* XXX: may need to worry about periodic boundary */
    return distance - lc->speedfactor * HorizonDistance(a, lc->horizon);
}

double
fastpm_lc_distance(FastPMLightCone * lc, double x[3])
{
    double distance;
    if (lc->fov <= 0) {
        distance = x[2];
    } else {
        distance = 0;
        int d;
        for (d = 0; d < 3; d ++) {
            distance += x[d] * x[d];
        }
        distance = sqrt(distance);
    }
    return distance;
}

static int
_fastpm_usmesh_intersect_one(FastPMUSMesh * mesh,
        struct funct_params * params,
        ptrdiff_t i,
        double * solution)
{
    params->i = i;

    return fastpm_horizon_solve(mesh->lc->horizon,
        params->context,
        solution,
        params->a1, params->a2,
        funct, params);
}

static double
zangle(double * x) {
    double dxy = 0;
    double dz = x[2];
    dxy = x[0] * x[0] + x[1] * x[1];

    double rt = atan2(sqrt(dxy), dz) / M_PI * 180.;
    if (rt < 0) rt += 360.;
    return rt;
}

/* is a vector in the i-th octant.
 * The ordering of octants is in C ordering, with
 * the Z-axis the fast changing direction.
 *
 * Thus [0, 4) are x > 0.
 *
 * The tolerance is in the same unit as the vector.
 * */
static int
_in_octant(int i, double vec[3], double tol)
{
    static const int signs[][3] = {
        {1, 1, 1},
        {1, 1, -1},
        {1, -1, 1},
        {1, -1, -1},
        {-1, 1, 1},
        {-1, 1, -1},
        {-1, -1, 1},
        {-1, -1, -1},
    };
    int d;
    for(d = 0; d < 3; d ++) {
        double s = vec[d] * signs[i][d];
        if(s < -tol) {
            return 0;
        }
    }
    return 1;
}
int
fastpm_lc_inside(FastPMLightCone * lc, double vec[3])
{
    if(lc->fov > 0) {
        int r;
        if(lc->fov < 360) {
            r = zangle(vec) <= lc->fov * 0.5;
        } else {
            r = 1;
        }
        if(r) {
            double norm = 0;
            int d;
            for(d = 0; d < 3; d ++) {
                norm += vec[d] * vec[d];
            }
            norm = sqrt(norm);
            int i;
            for(i = 0; i < 8; i ++) {
                if(lc->octants[i] && _in_octant(i, vec, lc->tol * norm)) return 1;
            }
            return 0;
        } else {
            return 0;
        }
    } else {
        /* smesh takes care of this */
        return 1;
    }
}

/* Compute an approximate AABB of all particles
 * based on position at two times, assuming linear motion.
 * */
static void
fastpm_compute_bbox(FastPMStore * p,
        FastPMDriftFactor * drift,
        double a1,
        double a2,
        double padding,
        double xmin[3],
        double xmax[3],
        MPI_Comm comm
) {
    ptrdiff_t i;
    for(int d = 0; d < 3; d ++) {
        xmin[d] = 1e20;
        xmax[d] = -1e20;
    }
    for(i = 0; i < p->np; i ++) {
        int d;
        for(d = 0; d < 3; d ++) {
            double xo[3];
            fastpm_drift_one(drift, p, i, xo, a1);
            if(xo[d] < xmin[d]) xmin[d] = xo[d];
            if(xo[d] > xmax[d]) xmax[d] = xo[d];
            fastpm_drift_one(drift, p, i, xo, a2);
            if(xo[d] < xmin[d]) xmin[d] = xo[d];
            if(xo[d] > xmax[d]) xmax[d] = xo[d];
        }
    }
    // It is sufficient not doing a reduce; each rank can cull its own
    // volume.
    // MPI_Allreduce(MPI_IN_PLACE, xmin, 3, MPI_DOUBLE, MPI_MIN, comm);
    // MPI_Allreduce(MPI_IN_PLACE, xmax, 3, MPI_DOUBLE, MPI_MAX, comm);
    for(int d = 0; d < 3; d ++) {
        xmin[d] -= padding;
        xmax[d] += padding;
    }
}

#include "spherebox.h"
void
rotate_vec3(vec3 v3, double glmatrix[4][4], int translate) {
    double xi[4];
    double xo[4];
    for(int d = 0; d < 3; d ++) {
        xi[d] = v3.v[d];
    }
    if(translate) xi[3] = 1;
    else xi[3] = 0;
    fastpm_gldot(glmatrix, xi, xo);
    for(int d = 0; d < 3; d ++) {
        v3.v[d] = xo[d];
    }
}
/*
 * if a shell in (radius1, radius2) intersects a AABB box at xmin to xmax.
 *
 * Apply glmatrix and tileshift to the AABB first to transform it to the intended location.
 * */
int
fastpm_shell_intersects_bbox(
    double xmin[3],
    double xmax[3],
    double glmatrix[4][4],
    double tileshift[3],
    double radius1,
    double radius2
)
{
    box bbox = make_box(xmin, xmax);

    /* rotate all planes, and apply tileshift */
    for(int i = 0; i < 6; i ++) {
        rotate_vec3(bbox.planes[i].position, glmatrix, 1);
        rotate_vec3(bbox.planes[i].direction, glmatrix, 0);
        for(int d = 0; d < 3; d ++) {
            bbox.planes[i].position.v[d] += tileshift[d];
        }
    }
    for(int i = 0; i < 8; i ++) {
        rotate_vec3(bbox.corners[i], glmatrix, 1);
        for(int d = 0; d < 3; d ++) {
            bbox.corners[i].v[d] += tileshift[d];
        }
    }
    sphere s1 = make_sphere0(radius1);
    if (box_inside_sphere(bbox, s1)) return 0;
    sphere s2 = make_sphere0(radius2);
    if (sphere_outside_box(s2, bbox)) return 0;

    return 1;
}
/* FIXME:
 * the function shall take ai, af as input,
 *
 * the function shall be able to interpolate potential as
 * well as position and velocity.
 * We need a more general representation of '*drift' and '*kick'.
 *
 * */
static int
fastpm_usmesh_intersect_tile(FastPMUSMesh * mesh, double * tileshift,
        double a1,
        double a2,
        FastPMDriftFactor * drift,
        FastPMKickFactor * kick,
        FastPMStore * p,
        FastPMStore * pout
)
{
    FastPMLightCone * lc = mesh->lc;
    struct funct_params params = {
        .lc = lc,
        .drift = drift,
        .p = p,
        .i = 0,
        .a1 = a1,
        .a2 = a2,
    };

    int a1_is_outside = (params.a1 > mesh->amax) || (params.a1 < mesh->amin);
    int a2_is_outside = (params.a2 > mesh->amax) || (params.a2 < mesh->amin);

    if(a1_is_outside && a2_is_outside) {
        return 0;
    }

    {
        int d;
        for(d = 0; d < 3; d ++) {
            params.tileshift[d] = tileshift[d];
        }
    }

    params.tileshift[3] = 0;

    ptrdiff_t i;

    #pragma omp parallel firstprivate(params)
    {
        params.context = fastpm_horizon_solve_start();
        #pragma omp for
        for(i = 0; i < p->np; i ++) {
            double a_emit = 0;
            if(0 == _fastpm_usmesh_intersect_one(mesh, &params, i, &a_emit)) continue;

            /* the event is outside the region we care, skip */
            if(a_emit > mesh->amax || a_emit < mesh->amin) continue;

            double xi[4];
            double xo[4];
            int d;

            xi[3] = 1;
            if(p->v) {
                /* can we drift? if we are using a fixed grid there is no v. */
                fastpm_drift_one(drift, p, i, xi, a_emit);
            } else {
                for(d = 0; d < 3; d ++) {
                    xi[d] = p->x[i][d];
                }
            }
            for(d = 0; d < 4; d ++) {
                xi[d] += params.tileshift[d];
            }
            /* transform the coordinate */
            fastpm_gldot(lc->glmatrix, xi, xo);

            /* does it fall into the field of view? */
            if(!fastpm_lc_inside(lc, xo)) continue;

            /* A solution is found */
            /* move the particle and store it. */
            ptrdiff_t next;
            #pragma omp atomic capture
                next = pout->np++;

            if(next >= pout->np_upper) {
                continue;
            }

            /* copy the position if desired */
            if(pout->x) {
                for(d = 0; d < 3; d ++) {
                    pout->x[next][d] = xo[d];
                }
            }

            float vo[4];
            float vi[4];
            if(p->v) {
                /* can we kick? if we are using a fixed grid there is no v */
                fastpm_kick_one(kick, p, i, vi, a_emit);
                vi[3] = 0;
                /* transform the coordinate */
                fastpm_gldotf(lc->glmatrix, vi, vo);

                if(pout->v) {
                    for(d = 0; d < 3; d ++) {
                        /* convert to peculiar velocity a dx / dt in kms */
                        pout->v[next][d] = vo[d] * HubbleConstant / a_emit;
                    }
                }
            }
            if(pout->id)
                pout->id[next] = p->id[i];
            if(pout->aemit)
                pout->aemit[next] = a_emit;
            if(pout->mask)
                pout->mask[next] = p->mask[i];

            double potfactor = 1.5 * lc->cosmology->Omega_cdm / (HubbleDistance * HubbleDistance);
            /* convert to dimensionless potential */
            if(pout->potential)
                pout->potential[next] = p->potential[i] / a_emit * potfactor;

            if(pout->tidal) {
                for(d = 0; d < 6; d++) {
                    pout->tidal[next][d] = p->tidal[i][d] / a_emit * potfactor;
                }
            }
        }
        fastpm_horizon_solve_end(params.context);
    }

    if(pout->np >= pout->np_upper) {
        fastpm_raise(-1, "Too many particles in the light cone; limit = %td, wanted = %td\n", pout->np_upper, pout->np);
    }
    return 0;
}

void
fastpm_usmesh_emit(FastPMUSMesh * mesh, int whence)
{
    /* a portion of light cone is ready between a0 and a1 */
    FastPMLCEvent lcevent[1];

    lcevent->p = mesh->p;
    lcevent->ai = mesh->ai;
    lcevent->af = mesh->af;
    lcevent->whence = whence;

    fastpm_emit_event(mesh->event_handlers,
            FASTPM_EVENT_LC_READY, FASTPM_EVENT_STAGE_AFTER,
            (FastPMEvent*) lcevent, mesh);

}

int
fastpm_usmesh_intersect(FastPMUSMesh * mesh, FastPMDriftFactor * drift, FastPMKickFactor * kick,
    double a1, double a2, int whence, MPI_Comm comm)
{
    {
        double t1 = a1, t2 = a2;
        a1 = fmin(t1, t2);
        a2 = fmax(t1, t2);
    }
    CLOCK(intersect);
    if (whence == TIMESTEP_START) {
        mesh->ai = a1;
        mesh->af = a1;
        fastpm_info("usmesh start event from %0.4f to %0.4f.\n", mesh->ai, mesh->af);
        mesh->np_before = 0;
        fastpm_usmesh_emit(mesh, whence);
    } else
    if (whence == TIMESTEP_CUR) {
        double xmin[3] = {0};
        double xmax[3] = {0};
        double padding = 0.5; /* rainwoodman: add a 500 Kpc/h padding; need a better estimate. */
        fastpm_compute_bbox(mesh->source, drift, a1, a2, padding, xmin, xmax, comm);

        fastpm_info("usmesh: bounding box computed for a = (%g, %g), AABB = [ %g %g %g ] - [ %g %g %g]",
                a1, a2,
                xmin[0], xmin[1], xmin[2],
                xmax[0], xmax[1], xmax[2]
            );

        /* for each tile */
        int t;
        ENTER(intersect);
        fastpm_info("usmesh intersection from %0.4f to %0.4f with %d tiles.\n", a1, a2, mesh->ntiles);

        double r1 = mesh->lc->speedfactor * HorizonDistance(a1, mesh->lc->horizon);
        double r2 = mesh->lc->speedfactor * HorizonDistance(a2, mesh->lc->horizon);

        double volume = 4 * M_PI / 3 * (pow(r1, 3) - pow(r2, 3));

        int steps = (int) ((volume / mesh->target_volume) + 0.5);
        if(steps == 0) steps = 1;
        double da = (a1 - a2) / steps;

        for(int i = 0; i < steps; i ++) {
            double ai = a1 + da * i;
            double af = (i + 1 == steps)?a2 : a1 + da * (i + 1);
            fastpm_info("usmesh: intersection step %d / %d a = %g %g .\n", i, steps, ai, af);

            int ntiles = 0;
            for(t = 0; t < mesh->ntiles; t ++) {
                /* for spherical geometry, skip if the tile does not intersects the lightcone. */
                if(mesh->lc->fov > 0 && !fastpm_shell_intersects_bbox(
                    xmin, xmax, mesh->lc->glmatrix, &mesh->tileshifts[t][0], r2, r1)) {
                    continue;
                }
                fastpm_usmesh_intersect_tile(mesh, &mesh->tileshifts[t][0],
                        ai, af,
                        drift, kick,
                        mesh->source,
                        mesh->p); /*Store particle to get density*/
                ntiles ++;
            }
            double ntiles_max, ntiles_min, ntiles_mean;
            MPIU_stats(comm, ntiles, "<->", &ntiles_min, &ntiles_mean, &ntiles_max);
            fastpm_info("usmesh: number of bounding box intersects shell r = (%g %g), min = %g max = %g, mean=%g",
                  r2, r1, ntiles_min, ntiles_max, ntiles_mean);
            LEAVE(intersect);
            mesh->af = af;
            if(MPIU_Any(comm, mesh->p->np > 0.5 * mesh->p->np_upper)) {
                fastpm_info("usmesh cur event from %0.4f to %0.4f.\n", mesh->ai, mesh->af);
                fastpm_usmesh_emit(mesh, whence);
                mesh->np_before += mesh->p->np;
                /* now purge the store. */
                mesh->p->np = 0;
                mesh->ai = mesh->af;
            }
        }
    } else
    if (whence == TIMESTEP_END) {
        mesh->af = a2;
        fastpm_info("usmesh end event from %0.4f to %0.4f.\n", mesh->ai, mesh->af);
        fastpm_usmesh_emit(mesh, whence);
        mesh->np_before += mesh->p->np;
        /* now purge the store. */
        mesh->p->np = 0;
        mesh->ai = mesh->af;
    }
    return 0;
}

ext: c

max_stars_repo_name: sbird/FastPMRunner

max_issues_repo_name: sbird/FastPMRunner

max_forks_repo_name: sbird/FastPMRunner


text:
#ifndef LIC_CALL_STACK_FRAME
#define LIC_CALL_STACK_FRAME

#include <stack>
#include <gsl/gsl>

#include "format/Opcodes.h"
#include "loader/MethodDefinition.h"
#include "EvaluationStack.h"

namespace lic
{

class CallStackFrame
{
public:
    CallStackFrame(MethodDefinition& method, gsl::span<TypedValue> args, Opcode* returnAddress);
    CallStackFrame(MethodDefinition& method, TypedValue& thisRef, gsl::span<TypedValue> args, Opcode* returnAddress);

    MethodDefinition& Method();
    TypedValue& Arg(size_t index);
    TypedValue& Local(size_t index);
    EvaluationStack& Stack();
    Opcode* ReturnAddress() const;

private:
    MethodDefinition& method;
    std::vector<TypedValue> args;
    std::vector<TypedValue> locals;
    EvaluationStack stack;
    Opcode* returnAddress;
};

}

#endif // !LIC_CALL_STACK_FRAME

max_stars_repo_name: roberthusak/lic

max_issues_repo_name: roberthusak/lic

max_forks_repo_name: roberthusak/lic


text:
#pragma once

#include <gsl/gsl>

max_stars_repo_name: DakkyDaWolf/OpenGL

max_issues_repo_name: DakkyDaWolf/OpenGL

max_forks_repo_name: DakkyDaWolf/OpenGL


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

#include <gsl/gsl_matrix.h>

/**
* @brief Class describing the Standard Model
*
* This class holds all the parameters specifying the Standard Model, and
* methods to access them in convenient ways.
*
*/
class SM {

 public:

  /// Default value for EM coupling \f$ \alpha_{\rm{EM}}(M_Z) \f$
  const static double alpha    = 1./127.934;

  /// Default value for EM coupling at low energy
  const static double alpha0   = 1./137.0359997;

  /// Default value for the Fermi constant \f$ G_F \f$ (\f$ \rm{GeV}^{-2} \f$)
  const static double GF      = 1.16637E-5;

  /// Default value for \f$ M_Z \f$ (GeV)
  const static double MZ      = 91.15349;

  /// Default value for \f$ M_W \f$ (GeV)
  const static double MW      = 80.36951;

  /// Default value for the strong coupling \f$ \alpha_s(M_Z) \f$
  const static double alpha_s = 0.119;

  /// Default value for total width \f$ \Gamma_Z \f$ (GeV)
  const static double GammaZ  = 2.49581;

  /// Default value for total width \f$ \Gamma_W \f$ (GeV)
  const static double GammaW  = 2.08856;

  /// Default value for quark pole mass \f$ m_d \f$ (GeV)
  const static double md_p    = 0.0;

  /// Default value for quark pole mass \f$ m_u \f$ (GeV)
  const static double mu_p    = 0.0;

  /// Default value for quark pole mass \f$ m_s \f$ (GeV)
  const static double ms_p    = 0.1;

  // Scale for strange quark mass (GeV)
  const static double Q_ms    = 2.0;

  /// Default value for quark pole mass \f$ m_c \f$ (GeV)
  const static double mc_p    = 1.42;

  /// Default value for quark pole mass \f$ m_b \f$ (GeV)
  const static double mb_p    = 4.75;

  /// Default value for quark pole mass \f$ m_t \f$ (GeV)
  const static double mt_p    = 172.5;

  /// Default value for electron pole mass \f$ m_e \f$ (GeV)
  const static double me_p    = 0.510998918E-3;

  /// Default value for muon pole mass \f$ m_\mu \f$ (GeV)
  const static double mmu_p   = 0.105658367;

  /// Default value for tau pole mass \f$ m_\tau \f$ (GeV)
  const static double mtau_p  = 1.77684;

  // PDG booklet 2012
  /// Default value for CKM matrix element \f$ V_{ud} \f$
  const static double Vud     = 0.97427;
  /// Default value for CKM matrix element \f$ V_{us} \f$
  const static double Vus     = 0.22534;
  /// Default value for CKM matrix element \f$ V_{ub} \f$
  const static double Vub     = 0.00351;
  /// Default value for CKM matrix element \f$ V_{cd} \f$
  const static double Vcd     = 0.22520;
  /// Default value for CKM matrix element \f$ V_{cs} \f$
  const static double Vcs     = 0.97334;
  /// Default value for CKM matrix element \f$ V_{cb} \f$
  const static double Vcb     = 0.0412;
  /// Default value for CKM matrix element \f$ V_{td} \f$
  const static double Vtd     = 0.00867;
  /// Default value for CKM matrix element \f$ V_{ts} \f$
  const static double Vts     = 0.0404;
  /// Default value for CKM matrix element \f$ V_{tb} \f$
  const static double Vtb     = 0.999146;

  /**
  * @brief Default constructor initializing the Standard Model with default parameter
  *        values
  *
  * This constructor initializes a new SM object with default parameters. The
  * values used as defaults are specified by public constants in this class.
  */
  SM();

  /**
  * @brief Sets the EM coupling \f$ \alpha \f$
  *
  * @param alpha_in New value for \f$ \alpha(M_Z) \f$
  */
  void set_alpha(double alpha_in);

  /**
  * @brief Sets the EM coupling \f$ \alpha_0 \f at low energy$
  *
  * @param alpha_in New value for \f$ \alpha(Q^2=0) \f$
  */
  void set_alpha0(double alpha_in);

  /**
  * @brief Sets the Fermi constant \f$ G_F \f$
  *
  * @param GF_in New value for \f$ G_F \f$ (\f$ GeV^{-2} \f$)
  */
  void set_GF(double GF_in);

  /**
  * @brief Sets \f$ M_Z \f$
  *
  * @param MZ_in New value for \f$ M_Z \f$
  */
  void set_MZ(double MZ_in);

  /**
  * @brief Sets \f$ M_W \f$
  *
  * @param MW_in New value for \f$ M_W \f$
  */
  void set_MW(double MW_in);

  /**
  * @brief Sets \f$ \gamma_Z \f$
  *
  * @param GammaZ_in New value for \f$ \Gamma_Z \f$
  */
  void set_gamma_Z(double GammaZ_in);

  /**
  * @brief Sets \f$ \gamma_W \f$
  *
  * @param GammaW_in New value for \f$ \Gamma_W \f$
  */
  void set_gamma_W(double GammaW_in);

  /**
  * @brief Sets the strong coupling \f$ \alpha_s \f$
  *
  * @param alpha_s_in New value for \f$ \alpha_s(M_Z) \f$
  */
  void set_alpha_s(double alpha_s_in);

  /**
  * @brief Current value of \f$ \alpha \f$
  *
  * @returns Value of \f$ \alpha \f$
  */
  double get_alpha();
  double get_alpha0();

  /**
  * @brief Current value of \f$ G_F \f$
  *
  * @returns Value of \f$ G_F \f$ (\f$ \rm{GeV}^{-2} \f$)
  */
  double get_GF();

  /**
  * @brief Current value of \f$ \sin\theta_W \f$
  *
  * @returns Value of \f$ \sin\theta_W \f$
  */
  double get_sintw();

  /**
  * @brief Current value of \f$ \cos\theta_W \f$
  *
  * @returns Value of \f$ \cos\theta_W \f$
  */
  double get_costw();

  /**
  * @brief Current value of weak SU(2) coupling \f$ g \f$
  *
  * @returns Value of \f$ g \f$
  */
  double get_g();

  /**
  * @brief Current value of hypercharge U(1) coupling \f$ g' \f$
  *
  * @returns Value of \f$ g' \f$
  */
  double get_gprime();

  /**
  * @brief Current value of EM coupling \f$ e \f$
  *
  * @returns Value of \f$ g' \f$
  */
  double get_e();

  /**
  * @brief Current value of \f$ \alpha_s \f$
  *
  * @returns Value of \f$ \alpha_s \f$
  */
  double get_alpha_s();

  /**
  * @brief Current value of \f$ M_W \f$
  *
  * @returns Value of \f$ M_W \f$ (GeV)
  */
  double get_MW();

  /**
  * @brief Current value of \f$ \Gamma_W \f$
  *
  * @returns Value of \f$ \Gamma_W \f$ (GeV)
  */
  double get_gamma_W();

  /**
  * @brief Current value of \f$ M_Z \f$
  *
  * @returns Value of \f$ M_Z \f$ (GeV)
  */
  double get_MZ();

  /**
  * @brief Current value of \f$ \Gamma_Z \f$
  *
  * @returns Value of \f$ \Gamma_Z \f$ (GeV)
  */
  double get_gamma_Z();

  /**
  * @brief Current value of the vev \f$ v=1/\sqrt{\sqrt{2}G_F} \f$
  *
  * @returns Value of \f$ v \f$ (GeV)
  */
  double get_v();

  /**
  * @brief Current value of the vev squared \f$ v^2=1/(\sqrt{2}G_F) \f$
  *
  * @returns Value of \f$ v^2 \f$ (\f$ \rm{GeV}^2 \f$)
  */
  double get_v2();

  /**
  * @brief Vector boson mass
  *
  * This method returns the mass of a given vector boson. The numbering
  * convention used throughout 2HDMC is such that (1 = photon, 2 = Z0, 3 = W+).
  *
  * @param v Index (1-3) for the vector boson
  *
  * @returns Value of \f$ m_V \f$ (GeV)
  */
  double get_vmass(int v);

  /**
  * @brief Vector boson width
  *
  * This method returns the decay width of a given vector boson. The numbering
  * convention is such that (1 = photon, 2 = Z0, 3 = W+).
  *
  * @param v Index (1-3) for the vector boson
  *
  * @returns Value of \f$ \Gamma_V \f$ (GeV)
  */
  double get_gamma_V(int v);

 /**
  * @brief Sets the CKM matrix to a 3x3 unit matrix
  *
  * This method sets the CKM matrix to a 3x3 unit matrix, i.e. eliminating
  * all charged current interactions between different generations.
  */
  void set_diagonal_CKM();

  /**
  * @brief Full CKM matrix
  *
  * This method can be used to retrieve the full CKM matrix as a gsl_matrix.
  *
  * @returns CKM matrix
  *
  * @see get_MD, get_MU, get_ML, get_CKM_element
  */
  gsl_matrix* get_CKM_matrix();

  /**
  * @brief Element of the CKM matrix
  *
  * This method can be used to retrieve a certain element of the CKM matrix.
  *
  * @param i Row index (1,2,3 = u,c,t) of CKM matrix element
  * @param j Column index (1,2,3 = d,s,b) of CKM matrix element
  * @returns The requested CKM element
  *
  * @see get_CKM
  */
  double get_CKM_element(int i, int j);

  /**
  * @brief Sets the pole mass of a lepton
  *
  * This method sets the pole mass of lepton \a l to the supplied value.
  *
  * @param l Index of lepton (1,2,3 = \f$ e,\mu,\tau \f$)
  * @param lmass_in New mass of lepton
  */
  void set_lmass_pole(int l, double lmass_in);

  /**
  * @brief Sets the pole mass of a quark
  *
  * This method sets the pole mass of quark \a q to the supplied value.
  *
  * @param q Index of quark (1,2,3,4,5,6 = \f$ d,u,s,c,b,t \f$)
  * @param qmass_in New mass of quark
  */
  void set_qmass_pole(int q, double qmass_in);

  void set_qmass_msbar(int flav, double qmass_in);

  /**
  * @brief Sets the pole mass of a down-type quark
  *
  * This method sets the pole mass of quark \a d to the supplied value.
  *
  * @param d Index of quark (1,2,3 = \f$ d,s,b \f$)
  * @param dmass_in New mass of quark
  */
  void set_dmass_pole(int d, double dmass_in);

  /**
  * @brief Sets the pole mass of an up-type quark
  *
  * This method sets the pole mass of quark \a u to the supplied value.
  *
  * @param u Index of quark (1,2,3 = \f$ u,c,t \f$)
  * @param umass_in New mass of quark
  */
  void set_umass_pole(int u, double umass_in);

  /**
  * @brief Quark pole mass
  *
  * @param q Index (1,2,3,4,5,6 = \f$ d,u,s,c,b,t \f$) of quark
  *
  * @returns The pole mass \f$ m_q \f$
  */
  double get_qmass_pole(int q);

  /**
  * @brief Down-type quark pole mass
  *
  * @param d Index (1,2,3 = \f$ d,s,b \f$) of down-type quark
  *
  * @returns Quark pole mass \f$ m_d \f$
  */
  double get_dmass_pole(int d);

  /**
  * @brief Up-type quark pole mass
  *
  * @param u Index (1,2,3 = \f$ u,c,t \f$) of up-type quark
  *
  * @returns Quark pole mass \f$ m_u \f$
  */
  double get_umass_pole(int u);

  /**
  * @brief Lepton pole mass
  *
  * @param l Index (1,2,3 = \f$ e,\mu,\tau \f$) of lepton
  *
  * @returns The pole mass \f$ m_l \f$
  */
  double get_lmass_pole(int l);

  /**
  * @brief Mass matrix for down-type quarks
  *
  * This method can be used to retrieve, in matrix form, the current pole masses
  * used for the down-type quarks. The masses are returned as the diagonal
  * elements of a 3x3 gsl_matrix that must be allocated by the user.
  *
  * @returns Diagonal mass matrix for the down-type quarks
  *
  * @see get_MU, get_ML, get_CKM
  */
  gsl_matrix* get_MD();

  /**
  * @brief Mass matrix for up-type quarks
  *
  * This method can be used to retrieve, in matrix form, the current pole masses
  * used for the up-type quarks. The masses are returned as the diagonal
  * elements of a 3x3 gsl_matrix that must be allocated by the user.
  *
  * @returns Diagonal mass matrix for the up-type quarks
  *
  * @see get_MD, get_ML, get_CKM
  */
  gsl_matrix* get_MU();

  /**
  * @brief Mass matrix for leptons
  *
  * This method can be used to retrieve, in matrix form, the current pole masses
  * used for the leptons. The masses are returned as the diagonal
  * elements of a 3x3 gsl_matrix that must be allocated by the user.
  *
  * @returns Diagonal mass matrix for the leptons
  *
  * @see get_MD, get_MU, get_CKM
  */
  gsl_matrix* get_ML();

  /**
  * @brief Quark \f$ \overline{\rm{MS}}\f$ mass
  *
  * This method determines the \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_q(\overline{m}_q) \f$
  *
  * @param q Index (1,2,3,4,5,6 = \f$ d,u,s,c,b,t \f$) of quark
  *
  * @returns Quark \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_q(\overline{m}_q) \f$
  *
  * @see run_qmass_MSbar
  */
  double get_qmass_MSbar(int q);

  /**
  * @brief Down-type quark \f$ \overline{\rm{MS}}\f$ mass
  *
  * This method determines the \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_d(\overline{m}_d) \f$
  *
  * @param d Index (1,2,3 = \f$ d,s,b \f$) of down-type quark
  *
  * @returns Quark \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_d(\overline{m}_d) \f$
  *
  * @see get_umass_MSbar, run_qmass_MSbar
  */
  double get_dmass_MSbar(int d);

  /**
  * @brief Up-type quark \f$ \overline{\rm{MS}}\f$ mass
  *
  * This method determines the \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_u(\overline{m}_u) \f$
  *
  * @param u Index (1,2,3 = \f$ u,c,t \f$) of up-type quark
  *
  * @returns Quark \f$ \overline{\rm{MS}}\f$ mass \f$ \overline{m}_u(\overline{m}_u) \f$
  *
  * @see get_dmass_MSbar, run_qmass_MSbar
  */
  double get_umass_MSbar(int u);

  /**
  * @brief Evaluates running quark masses
  *
  * This method evaluates the running quark_mass \a quark_mass, initially specified
  * at the input scale \a Qinit, at the new scale \a Qfin. The calculation is performed
  * in the \f$ \overline{\rm{MS}} \f$ scheme with variable number of active
  * flavours and threshold matching.
  *
  * @param quark_mass The input value for the running quark mass
  * @param Qinit      Starting scale
  * @param Qfin       Final scale at which the mass is to be evaluated
  * @param mtop       \f$ \overline{\rm{MS}} \f$ top mass (for thresholds)
  * @param mbot       \f$ \overline{\rm{MS}} \f$ bottom mass (for thresholds)
  *
  * @returns The running mass at the scale \a Qfin
  */
  double run_qmass_MSbar(double quark_mass, double Qinit, double Qfin, double mtop, double mbot);

  /**
  * @brief Evaluates the running strong coupling
  *
  * This method evaluates the running strong coupling \f$ \alpha_s \f$ in the
  * \f$ \overline{\rm{MS}} \f$ scheme at the scale specified by \a Q.
  *
  * @param Q          Scale at which to evaluate the running coupling
  * @param mtop       \f$ \overline{\rm{MS}} \f$ top mass (for thresholds)
  * @param mbot       \f$ \overline{\rm{MS}} \f$ bottom mass (for thresholds)
  *
  * @returns The strong coupling \f$ \alpha_s \f$ at the scale \a Q
  */
  double run_alphas_MSbar(double Q, double mtop, double mbot);

  /**
  * @brief Number of active flavours
  *
  * This method returns the number of active quark flavours to be used at
  * a certain mass scale when decoupling the heavier quarks. The thresholds
  * are determined from the quark \f$ \overline{\rm{MS}} \f$ masses.
  *
  * @param M The scale for which the number of flavours is desired
  *
  * @returns The number of quark flavours active at this scale
  */
  int get_Nactivef(double M);

  /**
  * @brief %SM width of the top quark
  *
  * Calculates the width of the top quark in the %SM, i.e. no contributions from
  * \f$ t\to H^+ b \f$ decays etc. are included.
  *
  * @returns The top width \f$ \Gamma_t \f$
  */
  double get_gamma_top();
  double get_gamma_tWd(int d);

  // --- These masses are for compatibility test with HDECAY ---
  // Scale for HD quark masses (HD)
  const static double Q_HD    = 5.0;

  // switch for HD quark masses (HD)
  const static bool b_HD    = false;

  /// Default value for running quark mass at 5 GeV \f$ m_s(5) \f$ (GeV)
  const static double ms_5    = 0.08160951656;

  /// Default value for running quark mass at 5 GeV \f$ m_c(5) \f$ (GeV)
  const static double mc_5    = 0.8960809463;

  /// Default value for running quark mass at 5 GeV \f$ m_b(5) \f$ (GeV)
  const static double mb_5    = 4.310679848;

  /// Default value for running quark mass at 5 GeV \f$ m_t(5) \f$ (GeV)
  const static double mt_5    = 246.5552058;
  // -------------------------------------------------------------

 private:

  // Internal variables which store the current values of the SM parameters
  double m_alpha;
  double m_alpha0;
  double m_GF;
  double m_s2tW;
  double m_alpha_s;
  double m_md_p[3];
  double m_mu_p[3];
  double m_ml_p[3];
  double m_MW;
  double m_MZ;
  double m_GammaW;
  double m_GammaZ;
  double m_CKM[3][3];

  double m_qmass_ms[7];

	// Support functions for the running mass calculation
  double m_threshold(double as);
  double RQ(double as, int nf, int loops);
  double QCD_beta(int c, int nf);

  void clear_lookup();

};

#endif


text:
/* rng/rand48.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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.
 */

#include <config.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_sys.h>
#include <gsl/gsl_rng.h>

/* This is the Unix rand48() generator. The generator returns the
   upper 32 bits from each term of the sequence,

   x_{n+1} = (a x_n + c) mod m

   using 48-bit unsigned arithmetic, with a = 0x5DEECE66D , c = 0xB
   and m = 2^48. The seed specifies the upper 32 bits of the initial
   value, x_1, with the lower 16 bits set to 0x330E.

   The theoretical value of x_{10001} is 244131582646046.

   The period of this generator is ? FIXME (probably around 2^48). */

static inline void rand48_advance (void *vstate);
static unsigned long int rand48_get (void *vstate);
static double rand48_get_double (void *vstate);
static void rand48_set (void *state, unsigned long int s);

static const unsigned short int a0 = 0xE66D ;
static const unsigned short int a1 = 0xDEEC ;
static const unsigned short int a2 = 0x0005 ;

static const unsigned short int c0 = 0x000B ;

typedef struct
  {
    unsigned short int x0, x1, x2;
  }
rand48_state_t;

static inline void
rand48_advance (void *vstate)
{
  rand48_state_t *state = (rand48_state_t *) vstate;

  /* work with unsigned long ints throughout to get correct integer
     promotions of any unsigned short ints */

  const unsigned long int x0 = (unsigned long int) state->x0 ;
  const unsigned long int x1 = (unsigned long int) state->x1 ;
  const unsigned long int x2 = (unsigned long int) state->x2 ;

  unsigned long int a ;

  a = a0 * x0 + c0 ;
  state->x0 = (a & 0xFFFF) ;

  a >>= 16 ;

  /* although the next line may overflow we only need the top 16 bits
     in the following stage, so it does not matter */

  a += a0 * x1 + a1 * x0 ;
  state->x1 = (a & 0xFFFF) ;

  a >>= 16 ;
  a += a0 * x2 + a1 * x1 + a2 * x0 ;
  state->x2 = (a & 0xFFFF) ;
}

static unsigned long int
rand48_get (void *vstate)
{
  unsigned long int x1, x2;

  rand48_state_t *state = (rand48_state_t *) vstate;
  rand48_advance (state) ;

  x2 = (unsigned long int) state->x2;
  x1 = (unsigned long int) state->x1;

  return (x2 << 16) + x1;
}

static double
rand48_get_double (void * vstate)
{
  rand48_state_t *state = (rand48_state_t *) vstate;

  rand48_advance (state) ;

  return (ldexp((double) state->x2, -16)
          + ldexp((double) state->x1, -32)
          + ldexp((double) state->x0, -48)) ;
}

static void
rand48_set (void *vstate, unsigned long int s)
{
  rand48_state_t *state = (rand48_state_t *) vstate;

  if (s == 0)  /* default seed */
    {
      state->x0 = 0x330E ;
      state->x1 = 0xABCD ;
      state->x2 = 0x1234 ;
    }
  else
    {
      state->x0 = 0x330E ;
      state->x1 = s & 0xFFFF ;
      state->x2 = (s >> 16) & 0xFFFF ;
    }

  return;
}

static const gsl_rng_type rand48_type =
{"rand48",                      /* name */
 0xffffffffUL,                  /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (rand48_state_t),
 &rand48_set,
 &rand48_get,
 &rand48_get_double
};

const gsl_rng_type *gsl_rng_rand48 = &rand48_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include "CesiumUtility/Library.h"
#include <cmath>
#include <cstdint>
#include <gsl/narrow>
#include <initializer_list>
#include <map>
#include <optional>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
#include <variant>
#include <vector>

namespace CesiumUtility {

struct JsonValueMissingKey : public std::runtime_error {
  JsonValueMissingKey(const std::string& key)
      : std::runtime_error(key + " is not present in Object") {}
};

struct JsonValueNotRealValue : public std::runtime_error {
  JsonValueNotRealValue()
      : std::runtime_error("this->value was not double, uint64_t or int64_t") {}
};

template <typename T, typename U>
constexpr T losslessNarrowOrDefault(U u, T defaultValue) noexcept {
  constexpr const bool is_different_signedness =
      (std::is_signed<T>::value != std::is_signed<U>::value);

  const T t = gsl::narrow_cast<T>(u);

  if (static_cast<U>(t) != u ||
      (is_different_signedness && ((t < T{}) != (u < U{})))) {
    return defaultValue;
  }

  return t;
}

/**
 * @brief A generic implementation of a value in a JSON structure.
 *
 * Instances of this class are used to represent the common `extras` field
 * of glTF elements that extend the the {@link ExtensibleObject} class.
 */
class CESIUMUTILITY_API JsonValue final {
public:
  /**
   * @brief The type to represent a `null` JSON value.
   */
  using Null = std::nullptr_t;

  /**
   * @brief The type to represent a `Bool` JSON value.
   */
  using Bool = bool;

  /**
   * @brief The type to represent a `String` JSON value.
   */
  using String = std::string;

  /**
   * @brief The type to represent an `Object` JSON value.
   */
  using Object = std::map<std::string, JsonValue>;

  /**
   * @brief The type to represent an `Array` JSON value.
   */
  using Array = std::vector<JsonValue>;

  /**
   * @brief Default constructor.
   */
  JsonValue() : value() {}

  /**
   * @brief Creates a `null` JSON value.
   */
  JsonValue(std::nullptr_t) : value(nullptr) {}

  /**
   * @brief Creates a `Number` JSON value.
   *
   * NaN and ±Infinity are represented as {@link JsonValue::Null}.
   */
  JsonValue(double v) {
    if (std::isnan(v) || std::isinf(v)) {
      value = nullptr;
    } else {
      value = v;
    }
  }

  /**
   * @brief Creates a `std::int64_t` JSON value (Widening conversion from
   * std::int8_t).
   */
  JsonValue(std::int8_t v) : value(static_cast<std::int64_t>(v)) {}

  /**
   * @brief Creates a `std::uint64_t` JSON value (Widening conversion from
   * std::uint8_t).
   */
  JsonValue(std::uint8_t v) : value(static_cast<std::uint64_t>(v)) {}

  /**
   * @brief Creates a `std::int64_t` JSON value (Widening conversion from
   * std::int16_t).
   */
  JsonValue(std::int16_t v) : value(static_cast<std::int64_t>(v)) {}

  /**
   * @brief Creates a `std::uint64_t` JSON value (Widening conversion from
   * std::uint16_t).
   */
  JsonValue(std::uint16_t v) : value(static_cast<std::uint64_t>(v)) {}

  /**
   * @brief Creates a `std::int64_t` JSON value (Widening conversion from
   * std::int32_t).
   */
  JsonValue(std::int32_t v) : value(static_cast<std::int64_t>(v)) {}

  /**
   * @brief Creates a `std::uint64_t` JSON value (Widening conversion from
   * std::uint32_t).
   */
  JsonValue(std::uint32_t v) : value(static_cast<std::uint64_t>(v)) {}

  /**
   * @brief Creates a `std::int64_t` JSON value.
   */
  JsonValue(std::int64_t v) : value(v) {}

  /**
   * @brief Creates a `std::uint64_t` JSON value.
   */
  JsonValue(std::uint64_t v) : value(v) {}

  /**
   * @brief Creates a `Bool` JSON value.
   */
  JsonValue(bool v) : value(v) {}

  /**
   * @brief Creates a `String` JSON value.
   */
  JsonValue(const std::string& v) : value(v) {}

  /**
   * @brief Creates a `String` JSON value.
   */
  JsonValue(std::string&& v) : value(std::move(v)) {}

  /**
   * @brief Creates a `String` JSON value.
   */
  JsonValue(const char* v) : value(std::string(v)) {}

  /**
   * @brief Creates an `Object` JSON value with the given properties.
   */
  JsonValue(const std::map<std::string, JsonValue>& v) : value(v) {}

  /**
   * @brief Creates an `Object` JSON value with the given properties.
   */
  JsonValue(std::map<std::string, JsonValue>&& v) : value(std::move(v)) {}

  /**
   * @brief Creates an `Array` JSON value with the given elements.
   */
  JsonValue(const std::vector<JsonValue>& v) : value(v) {}

  /**
   * @brief Creates an `Array` JSON value with the given elements.
   */
  JsonValue(std::vector<JsonValue>&& v) : value(std::move(v)) {}

  /**
   * @brief Creates an JSON value from the given initializer list.
   */
  JsonValue(std::initializer_list<JsonValue> v)
      : value(std::vector<JsonValue>(v)) {}

  /**
   * @brief Creates an JSON value from the given initializer list.
   */
  JsonValue(std::initializer_list<std::pair<const std::string, JsonValue>> v)
      : value(std::map<std::string, JsonValue>(v)) {}

  [[nodiscard]] const JsonValue*
  getValuePtrForKey(const std::string& key) const;
  [[nodiscard]] JsonValue* getValuePtrForKey(const std::string& key);

  /**
   * @brief Gets a typed value corresponding to the given key in the
   * object represented by this instance.
   *
   * If this instance is not a {@link JsonValue::Object}, returns
   * `nullptr`. If the key does not exist in this object, returns
   * `nullptr`. If the named value does not have the type T, returns
   * nullptr.
   *
   * @tparam T The expected type of the value.
   * @param key The key for which to retrieve the value from this object.
   * @return A pointer to the requested value, or nullptr if the value
   * cannot be obtained as requested.
   */
  template <typename T>
  const T* getValuePtrForKey(const std::string& key) const {
    const JsonValue* pValue = this->getValuePtrForKey(key);
    if (!pValue) {
      return nullptr;
    }

    return std::get_if<T>(&pValue->value);
  }

  /**
   * @brief Gets a typed value corresponding to the given key in the
   * object represented by this instance.
   *
   * If this instance is not a {@link JsonValue::Object}, returns
   * `nullptr`. If the key does not exist in this object, returns
   * `nullptr`. If the named value does not have the type T, returns
   * nullptr.
   *
   * @tparam T The expected type of the value.
   * @param key The key for which to retrieve the value from this object.
   * @return A pointer to the requested value, or nullptr if the value
   * cannot be obtained as requested.
   */
  template <typename T> T* getValuePtrForKey(const std::string& key) {
    JsonValue* pValue = this->getValuePtrForKey(key);
    return std::get_if<T>(&pValue->value);
  }

  /**
   * @brief Converts the numerical value corresponding to the given key
   * to the provided numerical template type.

   * If this instance is not a {@link JsonValue::Object}, throws
   * `std::bad_variant_access`. If the key does not exist in this object, throws
   * `JsonValueMissingKey`. If the named value does not have a numerical type T,
   *  throws `JsonValueNotRealValue`, if the named value cannot be converted
   from
   * `double` / `std::uint64_t` / `std::int64_t` without precision loss, throws
   * `gsl::narrowing_error`
   * @tparam To The expected type of the value.
   * @param key The key for which to retrieve the value from this object.
   * @return The converted value.
   * @throws If unable to convert the converted value for one of the
   aforementioned reasons.
   * @remarks Compilation will fail if type 'To' is not an integral / float /
   double type.
   */
  template <
      typename To,
      typename std::enable_if<
          std::is_integral<To>::value ||
          std::is_floating_point<To>::value>::type* = nullptr>
  [[nodiscard]] To getSafeNumericalValueForKey(const std::string& key) const {
    const Object& pObject = std::get<Object>(this->value);
    const auto it = pObject.find(key);
    if (it == pObject.end()) {
      throw JsonValueMissingKey(key);
    }
    return it->second.getSafeNumber<To>();
  }

  /**
   * @brief Converts the numerical value corresponding to the given key
   * to the provided numerical template type.
   *
   * If this instance is not a {@link JsonValue::Object}, the key does not exist
   * in this object, or the named value does not have a numerical type that can
   * be represented as T without precision loss, then the default value is
   * returned.
   *
   * @tparam To The expected type of the value.
   * @param key The key for which to retrieve the value from this object.
   * @return The converted value.
   * @throws If unable to convert the converted value for one of the
   * aforementioned reasons.
   * @remarks Compilation will fail if type 'To' is not an integral / float /
   * double type.
   */
  template <
      typename To,
      typename std::enable_if<
          std::is_integral<To>::value ||
          std::is_floating_point<To>::value>::type* = nullptr>
  [[nodiscard]] To getSafeNumericalValueOrDefaultForKey(
      const std::string& key,
      To defaultValue) const {
    const Object& pObject = std::get<Object>(this->value);
    const auto it = pObject.find(key);
    if (it == pObject.end()) {
      return defaultValue;
    }
    return it->second.getSafeNumberOrDefault<To>(defaultValue);
  }

  /**
   * @brief Determines if this value is an Object and has the given key.
   *
   * @param key The key.
   * @return true if this value contains the key. false if it is not an object
   * or does not contain the given key.
   */
  [[nodiscard]] inline bool hasKey(const std::string& key) const {
    const Object* pObject = std::get_if<Object>(&this->value);
    if (!pObject) {
      return false;
    }

    return pObject->find(key) != pObject->end();
  }

  /**
   * @brief Gets the numerical quantity from the value casted to the `To`
   * type. This function should be used over `getDouble()` / `getUint64()` /
   * `getInt64()` if you plan on casting that type into another smaller type or
   * different type.
   * @returns The converted type if it can be cast without precision loss.
   * @throws If the underlying value is not a numerical type or it cannot be
   *         converted without precision loss.
   */
  template <
      typename To,
      typename std::enable_if<
          std::is_integral<To>::value ||
          std::is_floating_point<To>::value>::type* = nullptr>
  [[nodiscard]] To getSafeNumber() const {
    const std::uint64_t* uInt = std::get_if<std::uint64_t>(&this->value);
    if (uInt) {
      return gsl::narrow<To>(*uInt);
    }

    const std::int64_t* sInt = std::get_if<std::int64_t>(&this->value);
    if (sInt) {
      return gsl::narrow<To>(*sInt);
    }

    const double* real = std::get_if<double>(&this->value);
    if (real) {
      return gsl::narrow<To>(*real);
    }

    throw JsonValueNotRealValue();
  }

  /**
   * @brief Gets the numerical quantity from the value casted to the `To`
   * type or returns defaultValue if unable to do so.

   * @returns The converted type if it can be cast without precision loss
   * or `defaultValue` if it cannot be converted safely.
   */
  template <
      typename To,
      typename std::enable_if<
          std::is_integral<To>::value ||
          std::is_floating_point<To>::value>::type* = nullptr>
  [[nodiscard]] To getSafeNumberOrDefault(To defaultValue) const noexcept {
    const std::uint64_t* uInt = std::get_if<std::uint64_t>(&this->value);
    if (uInt) {
      return losslessNarrowOrDefault<To>(*uInt, defaultValue);
    }

    const std::int64_t* sInt = std::get_if<std::int64_t>(&this->value);
    if (sInt) {
      return losslessNarrowOrDefault<To>(*sInt, defaultValue);
    }

    const double* real = std::get_if<double>(&this->value);
    if (real) {
      return losslessNarrowOrDefault<To>(*real, defaultValue);
    }

    return defaultValue;
  }

  /**
   * @brief Gets the object from the value.
   * @return The object.
   * @throws std::bad_variant_access if the underlying type is not a
   * JsonValue::Object
   */

  [[nodiscard]] inline const JsonValue::Object& getObject() const {
    return std::get<JsonValue::Object>(this->value);
  }

  /**
   * @brief Gets the string from the value.
   * @return The string.
   * @throws std::bad_variant_access if the underlying type is not a
   * JsonValue::String
   */
  [[nodiscard]] inline const JsonValue::String& getString() const {
    return std::get<String>(this->value);
  }

  /**
   * @brief Gets the array from the value.
   * @return The arrayj.
   * @throws std::bad_variant_access if the underlying type is not a
   * JsonValue::Array
   */
  [[nodiscard]] inline const JsonValue::Array& getArray() const {
    return std::get<JsonValue::Array>(this->value);
  }

  /**
   * @brief Gets the bool from the value.
   * @return The bool.
   * @throws std::bad_variant_access if the underlying type is not a
   * JsonValue::Bool
   */
  [[nodiscard]] inline bool getBool() const {
    return std::get<bool>(this->value);
  }

  /**
   * @brief Gets the double from the value.
   * @return The double.
   * @throws std::bad_variant_access if the underlying type is not a double
   */
  [[nodiscard]] inline double getDouble() const {
    return std::get<double>(this->value);
  }

  /**
   * @brief Gets the std::uint64_t from the value.
   * @return The std::uint64_t.
   * @throws std::bad_variant_access if the underlying type is not a
   * std::uint64_t
   */
  [[nodiscard]] std::uint64_t getUint64() const {
    return std::get<std::uint64_t>(this->value);
  }

  /**
   * @brief Gets the std::int64_t from the value.
   * @return The std::int64_t.
   * @throws std::bad_variant_access if the underlying type is not a
   * std::int64_t
   */
  [[nodiscard]] std::int64_t getInt64() const {
    return std::get<std::int64_t>(this->value);
  }

  /**
   * @brief Gets the bool from the value or returns defaultValue
   * @return The bool or defaultValue if this->value is not a bool.
   */
  [[nodiscard]] inline bool getBoolOrDefault(bool defaultValue) const {
    const auto* v = std::get_if<bool>(&this->value);
    if (v) {
      return *v;
    }

    return defaultValue;
  }

  /**
   * @brief Gets the string from the value or returns defaultValue
   * @return The string or defaultValue if this->value is not a string.
   */
  [[nodiscard]] inline const JsonValue::String
  getStringOrDefault(String defaultValue) const {
    const auto* v = std::get_if<JsonValue::String>(&this->value);
    if (v) {
      return *v;
    }

    return defaultValue;
  }

  /**
   * @brief Gets the double from the value or returns defaultValue
   * @return The double or defaultValue if this->value is not a double.
   */
  [[nodiscard]] inline double getDoubleOrDefault(double defaultValue) const {
    const auto* v = std::get_if<double>(&this->value);
    if (v) {
      return *v;
    }

    return defaultValue;
  }

  /**
   * @brief Gets the uint64_t from the value or returns defaultValue
   * @return The uint64_t or defaultValue if this->value is not a uint64_t.
   */
  [[nodiscard]] inline std::uint64_t
  getUint64OrDefault(std::uint64_t defaultValue) const {
    const auto* v = std::get_if<std::uint64_t>(&this->value);
    if (v) {
      return *v;
    }

    return defaultValue;
  }

  /**
   * @brief Gets the int64_t from the value or returns defaultValue
   * @return The int64_t or defaultValue if this->value is not a int64_t.
   */
  [[nodiscard]] inline std::int64_t
  getInt64OrDefault(std::int64_t defaultValue) const {
    const auto* v = std::get_if<std::int64_t>(&this->value);
    if (v) {
      return *v;
    }

    return defaultValue;
  }

  /**
   * @brief Returns whether this value is a `null` value.
   */
  [[nodiscard]] inline bool isNull() const noexcept {
    return std::holds_alternative<Null>(this->value);
  }

  /**
   * @brief Returns whether this value is a `double`, `std::uint64_t` or
   * `std::int64_t`. Use this function in conjunction with `getNumber` for
   *  safely casting to arbitrary types
   */
  [[nodiscard]] inline bool isNumber() const noexcept {
    return isDouble() || isUint64() || isInt64();
  }

  /**
   * @brief Returns whether this value is a `Bool` value.
   */
  [[nodiscard]] inline bool isBool() const noexcept {
    return std::holds_alternative<Bool>(this->value);
  }

  /**
   * @brief Returns whether this value is a `String` value.
   */
  [[nodiscard]] inline bool isString() const noexcept {
    return std::holds_alternative<String>(this->value);
  }

  /**
   * @brief Returns whether this value is an `Object` value.
   */
  [[nodiscard]] inline bool isObject() const noexcept {
    return std::holds_alternative<Object>(this->value);
  }

  /**
   * @brief Returns whether this value is an `Array` value.
   */
  [[nodiscard]] inline bool isArray() const noexcept {
    return std::holds_alternative<Array>(this->value);
  }

  /**
   * @brief Returns whether this value is a `double` value.
   */
  [[nodiscard]] inline bool isDouble() const noexcept {
    return std::holds_alternative<double>(this->value);
  }

  /**
   * @brief Returns whether this value is a `std::uint64_t` value.
   */
  [[nodiscard]] inline bool isUint64() const noexcept {
    return std::holds_alternative<std::uint64_t>(this->value);
  }

  /**
   * @brief Returns whether this value is a `std::int64_t` value.
   */
  [[nodiscard]] inline bool isInt64() const noexcept {
    return std::holds_alternative<std::int64_t>(this->value);
  }

  /**
   * @brief The actual value.
   *
   * The type of the value may be queried with the `isNull`, `isDouble`,
   * `isBool`, `isString`, `isObject`, `isUint64`, `isInt64`, `isNumber`, and
   * `isArray` functions.
   *
   * The actual value may be obtained with the `getNumber`, `getBool`,
   * and `getString` functions for the respective types. For
   * `Object` values, the properties may be accessed with the
   * `getValueForKey` functions.
   */
  std::variant<
      Null,
      double,
      std::uint64_t,
      std::int64_t,
      Bool,
      String,
      Object,
      Array>
      value;
};
} // namespace CesiumUtility

max_stars_repo_name: zrkcode/cesium-native

max_issues_repo_name: zrkcode/cesium-native

max_forks_repo_name: zrkcode/cesium-native


text:
#pragma once

#include <errno.h>  // errno
#include <fcntl.h>  // open, write
#include <unistd.h> // getpid
#include <sstream>
#include <boost/filesystem.hpp>
#include <gsl/gsl_util>
#include "do_persistence.h"
#include "error_macros.h"

// RAII wrapper for creating, writing and deleting the file that will
// contain the port number for the REST interface.
class RestPortAdvertiser
{
public:
    RestPortAdvertiser(uint16_t port)
    {
        // We cannot remove the port file on shutdown because we have already dropped permissions.
        // Clean it up now when we are still running as root.
        _DeleteOlderPortFiles();

        std::stringstream ss;
        ss << docli::GetRuntimeDirectory() << '/' << _restPortFileNamePrefix << '.' << getpid();
        _outFilePath = ss.str();

        // If file already exists, it is truncated. Probably from an old instance that terminated abnormally.
        // Allow file to only be read by others.
        int fd = open(_outFilePath.data(), O_WRONLY | O_CREAT | O_TRUNC, S_IRUSR | S_IRGRP | S_IROTH);
        if (fd == -1)
        {
            THROW_HR_MSG(E_FAIL, "Failed to open file %s, errno: %d", _outFilePath.data(), errno);
        }

        try
        {
            const auto writeStr = std::to_string(port) + '\n';
            const auto cbWrite = gsl::narrow_cast<ssize_t>(writeStr.size() * sizeof(char));

            const ssize_t written = write(fd, writeStr.data(), cbWrite);
            if (written != cbWrite)
            {
                THROW_HR_MSG(E_FAIL, "Failed to write port, written: %zd, errno: %d", written, errno);
            }
        }
        catch (...)
        {
            (void)close(fd);
            throw;
        }
        (void)close(fd);
    }

    const std::string& OutFilePath() const { return _outFilePath; }

private:
    void _DeleteOlderPortFiles() try
    {
        auto& runtimeDirectory = docli::GetRuntimeDirectory();
        for (boost::filesystem::directory_iterator itr(runtimeDirectory); itr != boost::filesystem::directory_iterator(); ++itr)
        {
            auto& dirEntry = itr->path();
            if (dirEntry.filename().string().find(_restPortFileNamePrefix) != std::string::npos)
            {
                boost::system::error_code ec;
                boost::filesystem::remove(dirEntry, ec);
                if (ec)
                {
                    DoLogWarning("Failed to delete old port file (%d, %s) %s", ec.value(), ec.message().data(), dirEntry.string().data());
                }
            }
        }
    } CATCH_LOG()

    std::string _outFilePath;
    static constexpr const char* const _restPortFileNamePrefix = "restport";
};

max_stars_repo_name: QPC-database/do-client

max_issues_repo_name: QPC-database/do-client

max_forks_repo_name: QPC-database/do-client


text:
/*
  10 Sept 08: Process_Error_Code is now static and thus private to
              this library.
*/

#include <cblas.h>
#include "lbl.h"
#include "ma57.h"

#ifdef __cplusplus
extern "C" {   /* To prevent C++ compilers from mangling symbols */
#endif

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Process_Error_Code"
static int Process_Error_Code( Ma57_Data *ma57, int nerror ) {

    LOGMSG( " [%3d]", nerror );
    /* Take appropriate action according to error code */
    /* Inflation factors are from the MA57 spec sheet  */
    switch( nerror ) {

    case 0:
      break;

      /* Warnings */
    case 1:
      LOGMSG( " Found and ignored %d indices out of range.\n", ma57->info[2] );
      break;

    case 2:
      LOGMSG( " Found and summed %d duplicate entries.\n", ma57->info[3] );
      break;

    case 3:
      LOGMSG( " Found duplicate and out-of-range indices.\n" );
      break;

    case 4:
      LOGMSG( " Matrix has rank %d, deficient.\n", ma57->info[24] );
      break;

    case 5:
      LOGMSG( " Found %d pivot sign changes in definite matrix.\n",
              ma57->info[25] );
      break;

    case 8:
      LOGMSG( " Infinity norm of solution found to be zero.\n" );
      break;

    case 10:
      LOGMSG( " Insufficient real space. Increased LFACT to %d.\n",
              ma57->lfact );
      break;

    case 11:
      LOGMSG( " Insufficient integer space. Increased LIFACT to %d.\n",
              ma57->lifact );
      break;

      /* Errors */
    case -1:
      LOGMSG( " Value of n is out of range: %d.\n", ma57->info[1] );
      break;

    case -2:
      LOGMSG( " Value of nz is out of range: %d.\n", ma57->info[1] );
      break;

    case -3:
      LOGMSG( " Adjusted size of array FACT to %d.\n", ma57->lfact );
      break;

    case -4:
      LOGMSG( " Adjusting size of array IFACT to %d.\n", ma57->lifact );
      break;

    case -5:
      LOGMSG( " Small pivot encountered at pivot step %d. Threshold = %lf.\n",
              ma57->info[1], ma57->cntl[1] );

    case -6:
      LOGMSG( " Change in pivot sign detected at pivot step %d.\n",
              ma57->info[1] );
      break;

    case -7:
      LOGMSG( " Erroneous sizing of array FACT or IFACT.\n" );
      break;

    case -8:
      LOGMSG( " Iterative refinement failed to converge.\n" );
      break;

    case -9:
      LOGMSG( " Error in user-supplied permutation array in component %d.\n",
              ma57->info[1] );
      break;

    case -10:
      LOGMSG( " Unknown pivoting strategy: %d\n.", ma57->info[1] );
      break;

    case -11:
      LOGMSG( " Size of RHS must be %d. Received %d.\n",
              ma57->n, ma57->info[1] );
      break;

    case -12:
      LOGMSG( " Invalid value of JOB (%d).\n", ma57->info[1] );
      break;

    case -13:
      LOGMSG( " Invalid number of iterative refinement steps (%d).\n",
              ma57->info[1] );
      break;

    case -14:
      LOGMSG( " Failed to estimate condition number.\n" );
      break;

    case -15:
      LOGMSG( " LKEEP has value %d, less than minimum allowed.\n",
              ma57->info[1] );
      break;

    case -16:
      LOGMSG( " Invalid number of RHS (%d).\n", ma57->info[1] );
      break;

    case -17:
      LOGMSG( " Increasing size of LWORK to %d.\n", ma57->lwork );
      break;

    case -18:
      LOGMSG( " MeTiS library not available or not found.\n" );
      break;

    default:
      LOGMSG( " Unrecognized flag from Factorize()." );
      nerror = -30;
    }
    return nerror;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "MA57_Initialize"
  Ma57_Data *Ma57_Initialize( int nz, int n, FILE *logfile ) {

    /* Call initialize subroutine MA57ID and set defaults */
    Ma57_Data *ma57 = (Ma57_Data *)LBL_Calloc( 1, sizeof(Ma57_Data) );
    ma57->logfile   = logfile;

    LOGMSG( " MA57 :: Initializing..." );

    ma57->n         = n;
    ma57->nz        = nz;
    ma57->fetched   = 0;
    ma57->irn       = (int *)LBL_Calloc( nz, sizeof(int) );
    ma57->jcn       = (int *)LBL_Calloc( nz, sizeof(int) );
    ma57->lkeep     = 5*n + nz + imax(n,nz) + 42 + n; // Add n to suggested val.
    ma57->keep      = (int *)LBL_Calloc( ma57->lkeep, sizeof(int) );
    ma57->iwork     = (int *)LBL_Calloc( 5*n, sizeof(int) );
    ma57->work      = NULL; // Will be initialized in Ma57_Solve()

    LOGMSG( " Calling ma57id..." );
    MA57ID( ma57->cntl, ma57->icntl ); // Initialize all parameters

    // Ensure some default parameters are appropriate
    ma57->icntl[0] = -1;  // Stream for error messages.
    ma57->icntl[1] = -1;  // Stream for warning messages.
    ma57->icntl[2] = -1;  // Stream for monitoring printing.
    ma57->icntl[3] = -1;  // Stream for printing of statistics.
    ma57->icntl[4] =  0;  // Verbosity: 0=none, 1=errors, 2=1+warnings,
                          //            3=2+monitor, 4=3+input,output
    ma57->icntl[5] =  5;  // Pivot selection strategy:
                          //  0: AMD using MC47
                          //  1: User-supplied pivot sequence
                          //  2: AMD with dense row strategy
                          //  3: MD as in MA27
                          //  4: MeTiS
                          //  5: Automatic (4 with fallback on 2)
    ma57->icntl[6] =  1;  // Numerical pivoting strategy
                          //  1: Do pivoting using value in cntl(1)
                          //  2: No pivoting; exit on sign change or zero
                          //  3: No pivoting; exit if pivot modulus < cntl(2)
                          //  4: No pivoting; alter pivots to have same sign
    ma57->icntl[7] =  1;  // Memory will be re-allocated if necessary
    ma57->icntl[15] = 0;  // No need to scale system before factorizing
    ma57->fetched = 0;

    LOGMSG( " done\n");
    return ma57;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_Analyze"
  int Ma57_Analyze( Ma57_Data *ma57 ) {

    int finished = 0, error;
    LOGMSG( " MA57 :: Analyzing..." );

    /* Unpack data structure and call MA57AD */
    while( !finished ) {
      LOGMSG( "\n         Calling ma57ad... " );

      MA57AD( &(ma57->n), &(ma57->nz), ma57->irn, ma57->jcn,
              &(ma57->lkeep), ma57->keep, ma57->iwork, ma57->icntl,
              ma57->info, ma57->rinfo );

      error = ma57->info[0];

      if( !error )
        finished = 1;
      else {
        error = Process_Error_Code( ma57, error );
        if( error != -3 && error != -4 && error != -8 && error != -14 )
          return error;
      }
    }

    // Allocate data for Factorize()
    ma57->lfact = ceil( LFACT_GROW * ma57->info[8] );
    ma57->fact = (double *)LBL_Calloc( ma57->lfact, sizeof(double) );
    ma57->lifact = ceil( LIFACT_GROW * ma57->info[9] );
    ma57->ifact = (int *)LBL_Calloc( ma57->lifact, sizeof(int) );
    LBL_Free( ma57->iwork );
    ma57->iwork = (int *)LBL_Calloc( ma57->n, sizeof(int) );
    ma57->work = NULL;

    LOGMSG( " done\n");
    return 0;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_Factorize"
  int Ma57_Factorize( Ma57_Data *ma57, double A[] ) {

    // To ensure consistency with the MA27 interface, the user must pass
    // the array A, *including* its zero-th element. Internally, A+1 is
    // passed to MA57.

    double *newFact;
    int    *newIfact;
    int     newSize, one = 1, zero = 0;
    int     finished = 0, error;
    LOGMSG( " MA57 :: Factorizing...\n" );

    /* Unpack data structure and call MA57BD */
    while( !finished ) {
      LOGMSG( "         Calling ma57bd... " );

      MA57BD( &(ma57->n), &(ma57->nz), A, ma57->fact, &(ma57->lfact),
              ma57->ifact, &(ma57->lifact), &(ma57->lkeep), ma57->keep,
              ma57->iwork, ma57->icntl, ma57->cntl, ma57->info, ma57->rinfo );

      error = ma57->info[0];

      if( !error ) {

        finished = 1;

      } else {

        error = Process_Error_Code( ma57, error );

        if( error == -3 || error == 10 ) {

          /* Resize real workspace */
          newSize = (error == 10) ? ma57->lfact : ma57->info[16];
          newSize = ceil( LFACT_GROW * newSize );
          LOGMSG("         Resizing real workspace for factors to %d", newSize);
          newFact = (double *)LBL_Calloc( newSize, sizeof(double) );
          MA57ED( &(ma57->n), &zero, ma57->keep, ma57->fact, &(ma57->lfact),
                  newFact, &newSize, ma57->ifact, &(ma57->lifact), NULL,
                  &zero, ma57->info );
          LBL_Free( ma57->fact );
          ma57->fact = newFact; newFact = NULL;
          ma57->lfact = newSize;

        } else if( error == -4 || error == 11 ) {

          /* Resize integer workspace */
          newSize = (error == 11) ? ma57->lifact : ma57->info[17];
          newSize = ceil( LIFACT_GROW * newSize );
          LOGMSG("         Resizing int workspace for factors to %d", newSize);
          newIfact = (int *)LBL_Calloc( newSize, sizeof(int) );
          MA57ED( &(ma57->n), &one, ma57->keep, ma57->fact, &(ma57->lfact),
                  NULL, &zero, ma57->ifact, &(ma57->lifact), newIfact,
                  &newSize, ma57->info );
          LBL_Free( ma57->ifact );
          ma57->ifact = newIfact; newIfact = NULL;
          ma57->lifact = newSize;

        } else {
          finished = 1;
          if( error < 0 ) return error;
        }
      }
    }

    LOGMSG("\n");
    LOGMSG("         %-23s: %6i   %-28s: %6i\n",
           "No. of 2x2 pivots"            , ma57->info[21],
           "Pivot step for modification"  , ma57->info[26]);
    LOGMSG("         %-28s: %6i   %-28s: %6i\n",
           "No. of negative e-vals"       , ma57->info[23],
           "No. of entries in factor"     , ma57->info[13]);
    LOGMSG("         %-28s: %6i   %-28s: %6i\n",
           "Rank of factorization"        , ma57->info[24],
           "No. of pivot sign changes"    , ma57->info[25]);
    LOGMSG("         done\n");
    return 0;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_Solve"
  int Ma57_Solve( Ma57_Data *ma57, double x[] ) {

    int one = 1, finished = 0, error;

    LOGMSG( " MA57 :: Solving..." );

    ma57->job = 1;
    ma57->lrhs = ma57->n;
    ma57->nrhs = 1;
    ma57->lwork = ma57->n * ma57->nrhs;
    if( !ma57->work ) {
      LOGMSG("\n         Sizing work array to size %d ", ma57->lwork);
      ma57->work = (double *)LBL_Calloc( ma57->lwork, sizeof(double) );
    }

    while( !finished ) {
      LOGMSG( "\n         Calling ma57cd... " );

      /* Unpack data structure and call MA57CD */
      MA57CD( &(ma57->job), &(ma57->n), ma57->fact, &(ma57->lfact),
              ma57->ifact, &(ma57->lifact), &ma57->nrhs, x, &(ma57->lrhs),
              ma57->work, &(ma57->lwork), ma57->iwork, ma57->icntl,
              ma57->info );

      error = ma57->info[0];
      if( error == -17 ) {
        LBL_Free( ma57->work );
        ma57->lwork = ceil( 1.2 * ma57->lwork );
        LOGMSG("\n         Resizing work array to size %d ", ma57->lwork);
        ma57->work = (double *)LBL_Calloc( ma57->lwork, sizeof(double) );
      } else
        finished = 1;
      if( error ) error = Process_Error_Code( ma57, error );
    }

    LOGMSG( " done\n" );
    return error;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_Refine"
  int Ma57_Refine( Ma57_Data *ma57, double x[], double rhs[],
                   double A[], int maxitref, int job ) {

    int error;

    LOGMSG( " MA57 :: Performing iterative refinement..." );

    ma57->job = job;
    /* For values of 'job' that demand the presence of the residual,
     * it should be placed by the user in ma57->residual prior to calling
     * this function.
     */

    /* Allocate work space */
    ma57->icntl[8] = imax( 1, maxitref );  // Number of refinement iterations
    ma57->icntl[9] = 1; // Return estimates of condition number
    LBL_Free( ma57->iwork );
    ma57->iwork = (int *)LBL_Calloc( ma57->n, sizeof(int) );

    LBL_Free( ma57->work );
    ma57->lwork = ma57->n;
    if( ma57->icntl[8] > 1 ) {
      ma57->lwork += 2 * ma57->n;
      if( ma57->icntl[9] > 0 )
        ma57->lwork += 2 * ma57->n;
    }
    ma57->work = (double *)LBL_Calloc( ma57->lwork, sizeof(double) );
    ma57->residual = (double *)LBL_Calloc( ma57->n, sizeof(double) );

    /* Perform iterative refinement */
    MA57DD( &(ma57->job), &(ma57->n), &(ma57->nz), A, ma57->irn,
            ma57->jcn, ma57->fact, &(ma57->lfact), ma57->ifact,
            &(ma57->lifact), rhs, x, ma57->residual, ma57->work,
            ma57->iwork, ma57->icntl, ma57->cntl, ma57->info,
            ma57->rinfo );

    error = ma57->info[0];
    if( error ) error = Process_Error_Code(ma57, ma57->info[0]);
    LOGMSG( " done\n" );
    return error;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_Finalize"
  void Ma57_Finalize( Ma57_Data *ma57 ) {

    /* Free allocated memory */
    LOGMSG( " MA57 :: Deallocating data arrays..." );
    LBL_Free( ma57->irn );
    LBL_Free( ma57->jcn );
    LBL_Free( ma57->keep );
    LBL_Free( ma57->iwork );
    LBL_Free( ma57->fact );
    LBL_Free( ma57->ifact );
    LBL_Free( ma57->work );
    if( ma57->residual ) LBL_Free(ma57->residual);
    LOGMSG( " done.\n" );
    free(ma57);
    return;
  }

  /* ================================================================= */
    /* Currently, this routine doesn't do error checking on the
       allowable parameter values. */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_set_int_parm"
  void Ma57_set_int_parm( Ma57_Data *ma57, int parm, int val ) {

      switch (parm) {
      case MA57_I_PIV_SELECTION:
          ma57->icntl[MA57_I_PIV_SELECTION] = val;
          LOGMSG(" Set %-20s = %10d\n","I_PIV_SELECTION",val);
          break;
      case MA57_I_PIV_NUMERICAL:
          ma57->icntl[MA57_I_PIV_NUMERICAL] = val;
          LOGMSG(" Set %-20s = %10d\n","I_PIV_NUMERICAL",val);
          break;
      case MA57_I_SCALING:
          ma57->icntl[MA57_I_SCALING] = val;
          LOGMSG(" Set %-20s = %10d\n","I_SCALING",val);
          break;
      }
      return;
  }

  /* ================================================================= */
    /* Currently, this routine doesn't do error checking on the
       allowable parameter values. */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_set_real_parm"
  void Ma57_set_real_parm( Ma57_Data *ma57, int parm, double val ) {

      switch (parm) {
      case MA57_D_PIV_THRESH:
          ma57->cntl[MA57_D_PIV_THRESH] = val;
          LOGMSG(" Set %-20s = %10.2e\n","D_PIV_THRESH",val);
          break;
      case MA57_D_PIV_NONZERO:
          ma57->cntl[MA57_D_PIV_NONZERO] = val;
          LOGMSG(" Set %-20s = %10.2e\n","D_PIV_NONZERO",val);
          break;
      }
      return;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_get_int_parm"
  int Ma57_get_int_parm( Ma57_Data *ma57, int parm ) {

      switch (parm) {
      case MA57_I_PIV_SELECTION:
          return ma57->icntl[MA57_I_PIV_SELECTION];
      case MA57_I_PIV_NUMERICAL:
          return ma57->icntl[MA57_I_PIV_NUMERICAL];
      case MA57_I_SCALING:
          return ma57->icntl[MA57_I_SCALING];
      }
      return 0;
  }

  /* ================================================================= */

#ifdef  __FUNCT__
#undef  __FUNCT__
#endif
#define __FUNCT__ "Ma57_get_real_parm"
  double Ma57_get_real_parm( Ma57_Data *ma57, int parm ) {

      switch (parm) {
      case MA57_D_PIV_THRESH:
          return ma57->cntl[MA57_D_PIV_THRESH];
      case MA57_D_PIV_NONZERO:
          return ma57->cntl[MA57_D_PIV_NONZERO];
      }
      return 0.0;
  }

  /* ================================================================= */

#ifdef __cplusplus
}              /* Closing brace for  extern "C"  block */
#endif

max_stars_repo_name: fperignon/sandbox

max_issues_repo_name: fperignon/sandbox

max_forks_repo_name: fperignon/sandbox


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include <gsl/span>

#include "core/common/common.h"
#include "core/optimizer/graph_transformer.h"

namespace onnxruntime {

struct FreeDimensionOverride;

/**
@Class FreeDimensionOverrideTransformer

Transformer that overrides free dimensions in the graph with the specific value
that matches the denotation for that dimension.
*/
class FreeDimensionOverrideTransformer : public GraphTransformer {
 public:
  explicit FreeDimensionOverrideTransformer(gsl::span<const FreeDimensionOverride> overrides_to_apply);

 private:
  Status ApplyImpl(Graph& graph, bool& modified, int graph_level) const override;

  std::map<std::string, int64_t> dimension_override_by_denotation_;
};

}  // namespace onnxruntime