text:
/* Copyright (c) 2011-2019, NVIDIA CORPORATION. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *  * Neither the name of NVIDIA CORPORATION nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
#pragma once

#include <util.h>

#include <amg.h>
#include <basic_types.h>
#include <types.h>
#include <norm.h>
#include <logger.h>
#include <matrix_distribution.h>

#include <iostream>
#include <iomanip>
#include <blas.h>
#include <multiply.h>

#include <amg_level.h>
#include <amgx_c.h>
#include <profile.h>

#include <misc.h>
#include <string>
#include <cassert>
#include <csr_multiply.h>
#include <memory_info.h>

#include <thrust/sort.h>
#include <thrust/remove.h>
#include <thrust/unique.h>
#include <thrust/binary_search.h>
#include <thrust/iterator/constant_iterator.h>

#define COARSE_CLA_CONSO 0 // used enable / disable coarse level consolidation (used in cycles files)

namespace amgx
{

/**********************************************************
 * Glue ( Consolidation)
 *********************************************************/
#ifdef AMGX_WITH_MPI

//------------------------
//---------Matrix-------------
//------------------------
template<class TConfig>
void compute_glue_info(Matrix<TConfig> &A)
{
    // Fill distributed manager fields for consolidation
    // Example
    // destination_part = [0 0 0 0 4 4 4 4 8 8 8 8] (input from manager->computeDestinationPartitions)
    // num_parts_to_consolidate = 4 for partitions 0,4,8 - (0 otherwise)
    // parts_to_consolidate (rank 0)[0 1 2 3] (rank 4)[4 5 6 7] (rank 8)[8 9 10 11]
    //coarse_part_to_fine_part = [0 4 8] num_coarse_partitions = 3
    //fine_part_to_coarse_part = [0 0 0 0 1 1 1 1 2 2 2 2]
    //ConsolidationArrayOffsets constains the offset of the nnz of each partitions in row pointer fashion : 0, pat1.NNZ, pat1.NNZ+part2.NNZ ... NNZ
    typedef typename TConfig::template setMemSpace<AMGX_host>::Type TConfig_h;
    typedef typename TConfig_h::template setVecPrec<AMGX_vecInt>::Type ivec_value_type_h;
    typedef typename ivec_value_type_h::VecPrec VecInt_t;

    bool is_root_partition = false;
    int num_parts_to_consolidate = 0;
    int num_parts = A.manager->getComms()->get_num_partitions();
    int my_id = A.manager->global_id();
    Vector<ivec_value_type_h> parts_to_consolidate;
    Vector<ivec_value_type_h> dest_partitions = A.manager->getDestinationPartitions();

    // compute is_root_partition and num_parts_to_consolidate
    for (int i = 0; i < num_parts; i++)
    {
        if (dest_partitions[i] == my_id)
        {
            is_root_partition = true;
            num_parts_to_consolidate++;
        }
    }

    parts_to_consolidate.resize(num_parts_to_consolidate);
    // parts_to_consolidate
    int count = 0;

    for (int i = 0; i < num_parts; i++)
    {
        if (dest_partitions[i] == my_id)
        {
            parts_to_consolidate[count] = i;
            count++;
        }
    }

    A.manager->setIsRootPartition(is_root_partition);
    A.manager->setNumPartsToConsolidate(num_parts_to_consolidate);
    A.manager->setPartsToConsolidate(parts_to_consolidate);
    // We don't really use the following in the latest version of the glue path but they are useful information
    // coarse_to_fine_part, fine_to_coarse_part
    Vector<ivec_value_type_h> coarse_to_fine_part, fine_to_coarse_part(num_parts);
    coarse_to_fine_part = dest_partitions;
    thrust::sort(coarse_to_fine_part.begin(), coarse_to_fine_part.end());
    cudaCheckError();
    coarse_to_fine_part.erase(thrust::unique(coarse_to_fine_part.begin(), coarse_to_fine_part.end()), coarse_to_fine_part.end());
    cudaCheckError();
    thrust::lower_bound(coarse_to_fine_part.begin(), coarse_to_fine_part.end(), dest_partitions.begin(), dest_partitions.end(), fine_to_coarse_part.begin());
    cudaCheckError();
    A.manager->setCoarseToFine(coarse_to_fine_part);
    A.manager->setFineToCoarse(fine_to_coarse_part);
    Vector<ivec_value_type_h> consolidationArrayOffsets;
    consolidationArrayOffsets.resize(num_parts);
}

template<class TConfig>
MPI_Comm compute_glue_matrices_communicator(Matrix<TConfig> &A)
{
    // Create temporary communicators for each consilidated matrix
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    if (A.manager->getDestinationPartitions().size() != 0)
    {
        int color = A.manager->getDestinationPartitions()[rank];
        MPI_Comm new_comm;
        // Split a communicator into multiple, non-overlapping communicators by color(each destiation partition has its color)
        // 1. Use MPI_Allgather to get the color and key from each process
        // 2. Count the number of processes with the same color; create a
        //    communicator with that many processes.  If this process has
        //    MPI_UNDEFINED as the color, create a process with a single member.
        // 3. Use key to order the ranks
        MPI_Comm_split(MPI_COMM_WORLD, color, rank, &new_comm);
        return new_comm;
    }
    else
    {
        FatalError("NO DESTINATION PARTIONS", AMGX_ERR_CORE);
    }
}

//function object (functor) for thrust calls (it is a unary operator to add a constant)
template<typename T>
class add_constant_op
{
        const T c;
    public:
        add_constant_op(T _c) : c(_c) {}
        __host__ __device__ T operator()(const T &x) const
        {
            return x + c;
        }
};

template<class TConfig>
int create_part_offsets(MPI_Comm &mpicm, Matrix<TConfig> &nv_mtx)
{
    /* WARNING: Notice that part_offsets_h & part_offsets have type int64_t.
                Therefore we need to use MPI_INT64_T (or MPI_LONG_LONG) in MPI_Allgather.
                Also, we need the send & recv buffers to be of the same type, therefore
                we will create a temporary variable n64 of the correct type below. */
    //create TConfig64, which is the same as TConfig, but with index type being int64_t
    typedef typename TConfig::template setVecPrec<AMGX_vecInt64>::Type TConfig64;
    typedef typename TConfig64::VecPrec t_VecPrec; //t_VecPrec = int64_t
    int n, offset, mpist;
    int nranks = 0; //nv_mtx.manager->get_num_partitions();

    if (nv_mtx.manager != NULL)
    {
        //some initializations
        nv_mtx.getOffsetAndSizeForView(OWNED, &offset, &n);
        MPI_Comm_size(mpicm, &nranks);
        nv_mtx.manager->part_offsets_h.resize(nranks + 1);
        //gather the number of rows per partition on the host (on all ranks)
        t_VecPrec n64 = n;
        nv_mtx.manager->part_offsets_h[0] = 0; //first element is zero (the # of rows is gathered afterwards)

        if (typeid(t_VecPrec) == typeid(int64_t))
        {
            mpist = MPI_Allgather(&n64, 1, MPI_INT64_T, nv_mtx.manager->part_offsets_h.raw() + 1, 1, MPI_INT64_T, mpicm);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        //perform a prefix sum
        thrust::inclusive_scan(nv_mtx.manager->part_offsets_h.begin(), nv_mtx.manager->part_offsets_h.end(), nv_mtx.manager->part_offsets_h.begin());
        //create the corresponding array on device (this is important)
        nv_mtx.manager->part_offsets.resize(nranks + 1);
        thrust::copy(nv_mtx.manager->part_offsets_h.begin(), nv_mtx.manager->part_offsets_h.end(), nv_mtx.manager->part_offsets.begin());
    }

    return 0;
}

template<class TConfig>
int glue_matrices(Matrix<TConfig> &nv_mtx, MPI_Comm &nv_mtx_com, MPI_Comm &temp_com)
{
    typedef typename TConfig::IndPrec t_IndPrec;
    typedef typename TConfig::MatPrec t_MatPrec;
    int n, nnz, offset, l, k = 0, i;
    int start, end, shift;
    int mpist, root = 0;
    //MPI call parameters
    t_IndPrec *rc_ptr, *di_ptr;
    t_IndPrec *hli_ptr, *hgi_ptr, *hgr_ptr, *i_ptr, *r_ptr;
    t_MatPrec *hlv_ptr, *hgv_ptr, *v_ptr;
    thrust::host_vector<t_IndPrec> rc;
    thrust::host_vector<t_IndPrec> di;
    //unpacked local matrix on the device and host
    device_vector_alloc<t_IndPrec> Bp;
    device_vector_alloc<t_IndPrec> Bi;
    device_vector_alloc<t_MatPrec> Bv;
    thrust::host_vector<t_IndPrec> hBp;
    thrust::host_vector<t_IndPrec> hBi;
    thrust::host_vector<t_MatPrec> hBv;
    //Consolidated matrices on the host
    thrust::host_vector<t_IndPrec> hAp;
    thrust::host_vector<t_IndPrec> hAi;
    thrust::host_vector<t_MatPrec> hAv;
    //Consolidated matrices on the device
    device_vector_alloc<t_IndPrec> Ap;
    device_vector_alloc<t_IndPrec> Ai;
    device_vector_alloc<t_MatPrec> Av;
    //WARNING: this routine currently supports matrix only with block size =1 (it can be generalized in the future, though)
    //initialize the defaults
    mpist = MPI_SUCCESS;

    if (nv_mtx.manager != NULL)
    {
        //int rank = nv_mtx.manager->global_id();
        nv_mtx.getOffsetAndSizeForView(OWNED, &offset, &n);
        nv_mtx.getNnzForView(OWNED, &nnz);

        if (nv_mtx.manager->part_offsets_h.size() == 0 || nv_mtx.manager->part_offsets.size() == 0)   // create part_offsets_h & part_offsets
        {
            create_part_offsets(nv_mtx_com, nv_mtx);  // (if needed for aggregation path)
        }

        Bp.resize(n + 1);
        Bi.resize(nnz);
        Bv.resize(nnz);
        hBp.resize(n + 1);
        hBi.resize(nnz);
        hBv.resize(nnz);
        //--- unpack the matrix ---
        nv_mtx.manager->unpack_partition(thrust::raw_pointer_cast(Bp.data()),
                                         thrust::raw_pointer_cast(Bi.data()),
                                         thrust::raw_pointer_cast(Bv.data()));
        cudaCheckError();
        //copy to host (should be able to optimize this out later on)
        hBp = Bp;
        hBi = Bi;
        hBv = Bv;
        cudaCheckError();

        // --- Glue matrices ---

        // construct global row pointers
        // compute recvcounts and displacements for MPI_Gatherv
        if (nv_mtx.manager->isRootPartition())
        {
            l = nv_mtx.manager->getNumPartsToConsolidate();      // number of partitions
            rc.resize(l);
            di.resize(l);

            //compute recvcounts and displacements for MPI_Gatherv
            for (i = 0; i < l; i++)
            {
                start = nv_mtx.manager->part_offsets_h[nv_mtx.manager->getPartsToConsolidate()[i]];
                end   = nv_mtx.manager->part_offsets_h[nv_mtx.manager->getPartsToConsolidate()[i] + 1];
                rc[i] = end - start;
                di[i] = k + 1;
                k += rc[i];
            }

            hAp.resize(k + 1); // extra +1 is needed because row_offsets have one extra element at the end
        }

        cudaCheckError();
        //alias raw pointers to thrust vector data (see thrust example unwrap_pointer for details)
        rc_ptr = thrust::raw_pointer_cast(rc.data());
        di_ptr = thrust::raw_pointer_cast(di.data());
        hli_ptr = thrust::raw_pointer_cast(hBp.data());
        hgr_ptr = thrust::raw_pointer_cast(hAp.data());
        cudaCheckError();

        //gather (on the host)
        if (typeid(t_IndPrec) == typeid(int))
        {
            mpist = MPI_Gatherv(hli_ptr + 1, n, MPI_INT,  hgr_ptr, rc_ptr, di_ptr, MPI_INT, root, temp_com);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        // Adjust row pointers, construct global column indices and values (recvcounts and displacements were computed above)
        if (nv_mtx.manager->isRootPartition())
        {
            //adjust global row pointers and setup the recvcounts & displacements for subsequent MPI calls
            for (i = 0; i < l; i++)
            {
                start = di[i] - 1;
                end   = di[i] + rc[i] - 1;
                shift = hAp[start];
                thrust::transform(hAp.begin() + start + 1, hAp.begin() + end + 1, hAp.begin() + start + 1, add_constant_op<t_IndPrec>(shift));
                cudaCheckError();
                di[i] = shift;
                rc[i] = hAp[end] - hAp[start];
            }

            //some allocations/resizing
            hAi.resize(hAp[k]);
            hAv.resize(hAp[k]);
        }

        //alias raw pointers to thrust vector data (see thrust example unwrap_pointer for details)
        rc_ptr = thrust::raw_pointer_cast(rc.data());
        di_ptr = thrust::raw_pointer_cast(di.data());
        hli_ptr = thrust::raw_pointer_cast(hBi.data());
        hgi_ptr = thrust::raw_pointer_cast(hAi.data());
        hlv_ptr = thrust::raw_pointer_cast(hBv.data());
        hgv_ptr = thrust::raw_pointer_cast(hAv.data());
        cudaCheckError();

        //gather (on the host)
        //columns indices
        if (typeid(t_IndPrec) == typeid(int))
        {
            mpist = MPI_Gatherv(hli_ptr, nnz, MPI_INT,  hgi_ptr, rc_ptr, di_ptr, MPI_INT,  root, temp_com);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        //values
        if      (typeid(t_MatPrec) == typeid(float))
        {
            mpist = MPI_Gatherv(hlv_ptr, nnz, MPI_FLOAT,  hgv_ptr, rc_ptr, di_ptr, MPI_FLOAT,  root, temp_com);
        }
        else if (typeid(t_MatPrec) == typeid(double))
        {
            mpist = MPI_Gatherv(hlv_ptr, nnz, MPI_DOUBLE, hgv_ptr, rc_ptr, di_ptr, MPI_DOUBLE, root, temp_com);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        // --- Upload matrices ---
        if (nv_mtx.manager->isRootPartition())
        {
            n = hAp.size() - 1;
            nnz = hAi.size();
            Ap.resize(hAp.size());
            Ai.resize(hAi.size());
            Av.resize(hAv.size());
            thrust::copy(hAp.begin(), hAp.end(), Ap.begin());
            thrust::copy(hAi.begin(), hAi.end(), Ai.begin());
            thrust::copy(hAv.begin(), hAv.end(), Av.begin());
            cudaCheckError();
        }
        else
        {
            n = 0;
            nnz  = 0;
            Ap.resize(1); // warning row_ponter size is expected to be n+1.
            Ap.push_back(0);
            Ai.resize(0);
            Av.resize(0);
            cudaCheckError();
        }

        r_ptr = thrust::raw_pointer_cast(Ap.data());
        i_ptr = thrust::raw_pointer_cast(Ai.data());
        v_ptr = thrust::raw_pointer_cast(Av.data());
        cudaCheckError();
        upload_matrix_after_glue(n, nnz, r_ptr, i_ptr, v_ptr, nv_mtx);
    }
    else
    {
        /* ASSUMPTION: when manager has not been allocated you are running on a single rank */
    }

    return 0;
}

template<class TConfig>
int upload_matrix_after_glue(int n, int nnz, int *r_ptr, int *i_ptr, void *v_ptr, Matrix<TConfig> &nv_mtx)
{
    // Using a path similar to AMGX_matrix_upload_all_global
    typedef typename TConfig::IndPrec t_IndPrec;
    typedef typename TConfig::MatPrec t_MatPrec;
    typedef typename TConfig::template setMemSpace<AMGX_host>::Type TConfig_h;
    typedef typename TConfig::template setVecPrec<AMGX_vecInt>::Type ivec_value_type;
    typedef typename TConfig_h::template setVecPrec<AMGX_vecInt>::Type ivec_value_type_h;
    typedef typename ivec_value_type_h::VecPrec VecInt_t;
    typedef Vector<ivec_value_type> IVector;
    // some parameters
    int block_dimx, block_dimy, num_ranks, n_global, start, end, val;
    t_IndPrec *part_vec_ptr;
    thrust::host_vector<t_IndPrec> pv;
    // set parameters
    nv_mtx.setView(ALL); // not sure about this
    n_global = nv_mtx.manager->part_offsets_h.back();
    block_dimx = nv_mtx.get_block_dimx();
    block_dimy = nv_mtx.get_block_dimy();
    //MPI_Comm* mpi_comm = nv_mtx.getResources()->getMpiComm();
    //MPI_Comm_size(*mpi_comm, &num_ranks);
    num_ranks = nv_mtx.manager->getComms()->get_num_partitions();
    // WARNING We create an artificial partition vectior that matches the new distribution
    // example n = 8, num_ranks(ie. num partitions) = 4 , DestinationPartitions[0,0,2,2], partvec = [0,0,0,0,2,2,2,2]
    // This might be an issue for the finest level if input_partvect !=NULL
    Vector<ivec_value_type_h> dest_partitions = nv_mtx.manager->getDestinationPartitions();

    for (int i = 0; i < num_ranks; i++)
    {
        val = dest_partitions[i];
        start = nv_mtx.manager->part_offsets_h[i];
        end   = nv_mtx.manager->part_offsets_h[i + 1];

        for (int j = 0; j < end - start; j++)
        {
            pv.push_back(val);
        }
    }

    part_vec_ptr = thrust::raw_pointer_cast(pv.data());
    cudaCheckError();
    // Save some glue info
    bool is_root_partition = nv_mtx.manager->isRootPartition();
    int dest_part = nv_mtx.manager->getMyDestinationPartition();
    int num_parts_to_consolidate = nv_mtx.manager->getNumPartsToConsolidate();
    Vector<ivec_value_type_h> parts_to_consolidate = nv_mtx.manager->getPartsToConsolidate();
    std::vector<int> cao;

    for (int i = 0; i < nv_mtx.manager->part_offsets_h.size(); i++)
    {
        cao.push_back(nv_mtx.manager->part_offsets_h[i]);
    }

    // WARNING
    // renumbering contains the inverse permutation to unreorder an amgx vector
    // inverse_renumbering contains the permutaion to reorder an amgx vector
    Vector<ivec_value_type> ir = nv_mtx.manager->inverse_renumbering;
    Vector<ivec_value_type> r = nv_mtx.manager->renumbering;
    // We need that to exchange the halo of unglued vectors (in coarse level consolidation)
    Vector<ivec_value_type_h> nei = nv_mtx.manager->neighbors;  // just neighbors before glue
    std::vector<std::vector<VecInt_t> > b2lr = nv_mtx.manager->getB2Lrings(); //list of boundary nodes to export to other partitions.
    Vector<ivec_value_type_h> ho = nv_mtx.manager->halo_offsets;
    std::vector<IVector > b2lm = nv_mtx.manager->getB2Lmaps();
    cudaCheckError();

    // Create a fresh distributed manager
    if (nv_mtx.manager != NULL )
    {
        delete nv_mtx.manager;
    }

    nv_mtx.manager = new DistributedManager<TConfig>(nv_mtx);
    nv_mtx.set_initialized(0);
    nv_mtx.delProps(DIAG);
    // Load distributed matrix
    MatrixDistribution mdist;
    mdist.setPartitionVec(part_vec_ptr);
    nv_mtx.manager->loadDistributedMatrix(n, nnz, block_dimx, block_dimy, r_ptr, i_ptr, (t_MatPrec *) v_ptr, num_ranks, n_global, NULL, mdist);
    // Create B2L_maps for comm
    nv_mtx.manager->renumberMatrixOneRing();
    // WARNING WE SHOULD GET THE NUMBER OF RINGS AND DO THE FOLLOWING ONLY IF THERE ARE 2 RINGS
    // Exchange 1 ring halo rows (for d2 interp)
    // if (num_import_rings == 2)
    nv_mtx.manager->createOneRingHaloRows();
    nv_mtx.manager->getComms()->set_neighbors(nv_mtx.manager->num_neighbors());
    nv_mtx.setView(OWNED);
    nv_mtx.set_initialized(1);
    cudaCheckError();
    // restore some glue info to consolidate the vectors in the future
    nv_mtx.manager->setDestinationPartitions(dest_partitions);
    nv_mtx.manager->setIsRootPartition(is_root_partition);
    nv_mtx.manager->setNumPartsToConsolidate(num_parts_to_consolidate);
    nv_mtx.manager->setPartsToConsolidate(parts_to_consolidate);
    nv_mtx.manager->setIsGlued(true);
    nv_mtx.manager->setMyDestinationPartition(dest_part);
    nv_mtx.manager->setConsolidationArrayOffsets(cao); // partions_offest before consolidation
    // set fine level data structures, this is used to match the former distribution when we upload / download from the API

    if (nv_mtx.amg_level_index == 0)
    {
        // just small copies inside fineLevelUpdate.
        nv_mtx.manager->fineLevelUpdate();
    }

    nv_mtx.manager->renumbering_before_glue = r;
    nv_mtx.manager->inverse_renumbering_before_glue = ir;
    nv_mtx.manager->neighbors_before_glue = nei;
    nv_mtx.manager->halo_offsets_before_glue = ho;
    nv_mtx.manager->B2L_rings_before_glue = b2lr;
    nv_mtx.manager->B2L_maps_before_glue = b2lm;
    cudaCheckError();
    return 0;
}

template<class TConfig>
int glue_vector(Matrix<TConfig> &nv_mtx, MPI_Comm &A_comm, Vector<TConfig> &nv_vec, MPI_Comm &temp_com)
{
    // glu vecots based on dest_partitions (which contains partitions should be merged together)
    typedef typename TConfig::IndPrec t_IndPrec;
    typedef typename TConfig::VecPrec t_VecPrec;
    int n, l, mpist, start, end, k = 0, root = 0, rank = 0;
    //MPI call parameters
    t_IndPrec *rc_ptr, *di_ptr;
    t_VecPrec *hv_ptr, *hg_ptr, *v_ptr;
    thrust::host_vector<t_IndPrec> rc;
    thrust::host_vector<t_IndPrec> di;
    //unreordered local vector on the host
    thrust::host_vector<t_VecPrec> hv;
    //constructed global vector on the host
    thrust::host_vector<t_VecPrec> hg;
    //constructed global vector on the device
    device_vector_alloc<t_VecPrec> v;
    //WARNING: this routine currently supports vectors only with block size =1 (it can be generalized in the future, though)
    //initialize the defaults
    mpist = MPI_SUCCESS;

    if (nv_mtx.manager != NULL)
    {
        // some initializations
        rank = nv_mtx.manager->global_id();

        if (nv_mtx.manager->getComms() != NULL)
        {
            nv_mtx.manager->getComms()->get_mpi_comm();
        }

        n = nv_mtx.manager->getConsolidationArrayOffsets()[rank + 1] - nv_mtx.manager->getConsolidationArrayOffsets()[rank];

        if (nv_mtx.manager->getConsolidationArrayOffsets().size() == 0)
        {
            std::cout << "ERROR part_offsets in glue path" << std::endl;
        }

        l = nv_mtx.manager->getNumPartsToConsolidate();      // number of partitions
        //some allocations/resizing
        hv.resize(nv_mtx.manager->renumbering_before_glue.size());                         // host copy of nv_vec

        if (nv_mtx.manager->isRootPartition())
        {
            // This works with neighbours only
            hg.resize(nv_mtx.manager->getConsolidationArrayOffsets()[rank + l] - nv_mtx.manager->getConsolidationArrayOffsets()[rank]); // host copy of cvec
            rc.resize(l);
            di.resize(l);
        }

        cudaCheckError();
        //--- unreorder the vector back (just like you did with the matrix, but only need to undo the interior-boundary reordering, because others do not apply) ---
        // unreorder and copy the vector
        // WARNING
        // renumbering contains the inverse permutation to unreorder an amgx vector
        // inverse_renumbering contains the permutaion to reorder an amgx vector
        thrust::copy(thrust::make_permutation_iterator(nv_vec.begin(), nv_mtx.manager->renumbering_before_glue.begin()  ),
                     thrust::make_permutation_iterator(nv_vec.begin(), nv_mtx.manager->renumbering_before_glue.begin() +  nv_mtx.manager->renumbering_before_glue.size()),
                     hv.begin());
        cudaCheckError();
        hv.resize(n);

        // --- construct global vector (rhs/sol) ---
        //compute recvcounts and displacements for MPI_Gatherv
        if (nv_mtx.manager->isRootPartition())
        {
            l = nv_mtx.manager->getNumPartsToConsolidate();      // number of partitions

            //compute recvcounts and displacements for MPI_Gatherv
            for (int i = 0; i < l; i++)
            {
                start = nv_mtx.manager->getConsolidationArrayOffsets()[nv_mtx.manager->getPartsToConsolidate()[i]];
                end   = nv_mtx.manager->getConsolidationArrayOffsets()[nv_mtx.manager->getPartsToConsolidate()[i] + 1];
                rc[i] = end - start;
                di[i] = k;
                k += rc[i];
            }
        }

        //alias raw pointers to thrust vector data (see thrust example unwrap_pointer for details)
        rc_ptr = thrust::raw_pointer_cast(rc.data());
        di_ptr = thrust::raw_pointer_cast(di.data());
        hv_ptr = thrust::raw_pointer_cast(hv.data());
        hg_ptr = thrust::raw_pointer_cast(hg.data());
        cudaCheckError();

        //gather (on the host)
        if      (typeid(t_VecPrec) == typeid(float))
        {
            mpist = MPI_Gatherv(hv_ptr, n, MPI_FLOAT,  hg_ptr, rc_ptr, di_ptr, MPI_FLOAT,  root, temp_com);
        }
        else if (typeid(t_VecPrec) == typeid(double))
        {
            mpist = MPI_Gatherv(hv_ptr, n, MPI_DOUBLE, hg_ptr, rc_ptr, di_ptr, MPI_DOUBLE, root, temp_com);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        // clean
        nv_vec.in_transfer = IDLE;

        //nv_vec.dirtybit = 0;
        if (nv_vec.buffer != NULL)
        {
            delete nv_vec.buffer;
            nv_vec.buffer = NULL;
            nv_vec.buffer_size = 0;
        }

        if (nv_vec.linear_buffers_size != 0)
        {
            amgx::memory::cudaFreeHost(&(nv_vec.linear_buffers[0]));
            nv_vec.linear_buffers_size = 0;
        }

        if (nv_vec.explicit_host_buffer)
        {
            amgx::memory::cudaFreeHost(nv_vec.explicit_host_buffer);
            nv_vec.explicit_host_buffer = NULL;
            nv_vec.explicit_buffer_size = 0;
            cudaEventDestroy(nv_vec.mpi_event);
        }

        // resize
        if (nv_mtx.manager->isRootPartition())
        {
            n = hg.size();
            v.resize(hg.size());
            thrust::copy(hg.begin(), hg.end(), v.begin());
            cudaCheckError();
        }
        else
        {
            n = 0;
            v.resize(0);
            cudaCheckError();
        }

        // upload
        v_ptr = thrust::raw_pointer_cast(v.data());
        cudaCheckError();
        upload_vector_after_glue(n, v_ptr, nv_vec, nv_mtx);
    }
    else
    {
        // ASSUMPTION: when manager has not been allocated you are running on a single rank
    }

    return 0;
}

template<class TConfig>
int upload_vector_after_glue(int n, void *v_ptr, Vector<TConfig> &nv_vec, Matrix<TConfig> &nv_mtx)
{
    typedef typename TConfig::VecPrec t_VecPrec;
    // vector bind
    nv_vec.unsetManager();
    cudaCheckError();

    if (nv_mtx.manager != NULL)
    {
        nv_vec.setManager(*(nv_mtx.manager));
    }

    cudaCheckError();

    if (nv_vec.is_transformed())
    {
        nv_vec.unset_transformed();
    }

    nv_vec.set_block_dimx(1);
    nv_vec.set_block_dimy(nv_mtx.get_block_dimy());
    // the dirtybit has to be set to one here in order to have correct results to ensure an exachange halo before the solve
    // this is particulary important when the number of consolidated partitions is greater than 1 on large matrices such as drivaer9M
    nv_vec.dirtybit = 1;

    if (nv_mtx.manager != NULL)
    {
        nv_vec.getManager()->transformAndUploadVector(nv_vec, (t_VecPrec *)v_ptr, n, nv_vec.get_block_dimy());
    }

    MPI_Barrier(MPI_COMM_WORLD);
    return 0;
}

template<class TConfig>
int unglue_vector(Matrix<TConfig> &nv_mtx, MPI_Comm &A_comm, Vector<TConfig> &nv_vec, MPI_Comm &temp_com, Vector<TConfig> &nv_vec_unglued)
{
    // glue vecots based on dest_partitions (which contains partitions should be merged together)
    typedef typename TConfig::IndPrec t_IndPrec;
    typedef typename TConfig::VecPrec t_VecPrec;
    int n_loc, l, mpist, start, end, k = 0, root = 0, rank = 0;
    //MPI call parameters
    t_IndPrec *sc_ptr, *di_ptr;
    t_VecPrec *hv_ptr, *hg_ptr;
    thrust::host_vector<t_IndPrec> sc;
    thrust::host_vector<t_IndPrec> di;
    //unreordered local vector on the host
    thrust::host_vector<t_VecPrec> hv;
    //constructed global vector on the host
    thrust::host_vector<t_VecPrec> hg;
    //constructed global vector on the device
    device_vector_alloc<t_VecPrec> v;
    //WARNING: this routine currently supports vectors only with block size =1 (it can be generalized in the future, though)
    //initialize the defaults
    mpist = MPI_SUCCESS;

    if (nv_mtx.manager != NULL)
    {
        // some initializations
        rank = nv_mtx.manager->global_id();

        if (nv_mtx.manager->getComms() != NULL)
        {
            nv_mtx.manager->getComms()->get_mpi_comm();
        }

        n_loc = nv_mtx.manager->getConsolidationArrayOffsets()[rank + 1] - nv_mtx.manager->getConsolidationArrayOffsets()[rank];

        if (nv_mtx.manager->getConsolidationArrayOffsets().size() == 0)
        {
            printf("ERROR part_offsets\n");
        }

        l = nv_mtx.manager->getNumPartsToConsolidate();      // number of partitions
        //some allocations/resizing
        hv.resize(n_loc);                         // host copy

        if (nv_mtx.manager->isRootPartition())
        {
            hg.resize(nv_vec.size()); // host copy of cvec
        }

        sc.resize(l);
        di.resize(l);
        cudaCheckError();
        // Exchange_halo before unreordering
        // Do we need that?
        nv_mtx.manager->exchange_halo(nv_vec, nv_vec.tag);

        // unreorder the vector
        if (nv_mtx.manager->isRootPartition())
        {
            // WARNING
            // renumbering contains the inverse permutation to unreorder an amgx vector
            // inverse_renumbering contains the permutaion to reorder an amgx vector
            thrust::copy(thrust::make_permutation_iterator(nv_vec.begin(), nv_mtx.manager->renumbering.begin()  ),
                         thrust::make_permutation_iterator(nv_vec.begin(), nv_mtx.manager->renumbering.begin() + nv_mtx.manager->renumbering.size()),
                         hg.begin());
            cudaCheckError();
            hg.resize(nv_mtx.manager->getConsolidationArrayOffsets()[rank + l] - nv_mtx.manager->getConsolidationArrayOffsets()[rank]);
        }

        // --- construct local vector (sol) ---
        //compute sendcounts and displacements for MPI_Gatherv
        for (int i = 0; i < l; i++)
        {
            start = nv_mtx.manager->getConsolidationArrayOffsets()[nv_mtx.manager->getPartsToConsolidate()[i]];
            end   = nv_mtx.manager->getConsolidationArrayOffsets()[nv_mtx.manager->getPartsToConsolidate()[i] + 1];
            sc[i] = end - start;
            di[i] = k;
            k += sc[i];
        }

        //alias raw pointers to thrust vector data (see thrust example unwrap_pointer for details)
        sc_ptr = thrust::raw_pointer_cast(sc.data());
        di_ptr = thrust::raw_pointer_cast(di.data());
        hv_ptr = thrust::raw_pointer_cast(hv.data());
        hg_ptr = thrust::raw_pointer_cast(hg.data());
        cudaCheckError();

        // Scatter (on the host)
        if      (typeid(t_VecPrec) == typeid(float))
        {
            mpist = MPI_Scatterv(hg_ptr, sc_ptr, di_ptr, MPI_FLOAT,  hv_ptr, n_loc, MPI_FLOAT,  root, temp_com);
        }
        else if (typeid(t_VecPrec) == typeid(double))
        {
            mpist = MPI_Scatterv(hg_ptr, sc_ptr, di_ptr, MPI_DOUBLE, hv_ptr, n_loc, MPI_DOUBLE, root, temp_com);
        }
        else
        {
            FatalError("MPI_Gatherv of the vector has failed - incorrect vector data type", AMGX_ERR_CORE);
        }

        if (mpist != MPI_SUCCESS)
        {
            FatalError("MPI_Gatherv of the vector has failed - detected incorrect MPI return code", AMGX_ERR_CORE);
        }

        // --- Manual upload ---
        // Cleaning
        nv_vec_unglued.in_transfer = IDLE;

        if (nv_vec_unglued.buffer != NULL)
        {
            delete nv_vec_unglued.buffer;
            nv_vec_unglued.buffer = NULL;
            nv_vec_unglued.buffer_size = 0;
        }

        if (nv_vec_unglued.linear_buffers_size != 0)
        {
            amgx::memory::cudaFreeHost(&(nv_vec_unglued.linear_buffers[0]));
            nv_vec_unglued.linear_buffers_size = 0;
        }

        if (nv_vec_unglued.explicit_host_buffer)
        {
            amgx::memory::cudaFreeHost(nv_vec_unglued.explicit_host_buffer);
            nv_vec_unglued.explicit_host_buffer = NULL;
            nv_vec_unglued.explicit_buffer_size = 0;
            cudaEventDestroy(nv_vec_unglued.mpi_event);
        }

        // We should avoid copies between nv_vec and hv here
        nv_vec_unglued.resize( nv_mtx.manager->inverse_renumbering_before_glue.size());
        thrust::fill( nv_vec_unglued.begin(), nv_vec_unglued.end(), 0.0 );
        thrust::copy(hv.begin(), hv.end(), nv_vec_unglued.begin());
        hv.resize( nv_mtx.manager->inverse_renumbering_before_glue.size());
        cudaCheckError();
        //  Manual reordering
        //  Upload_vector_after_glue is not going to work because matrix managers has been modified during glued matrices, and don't match the new, glued, topology.
        thrust::copy(thrust::make_permutation_iterator(nv_vec_unglued.begin(), nv_mtx.manager->inverse_renumbering_before_glue.begin()  ),
                     thrust::make_permutation_iterator(nv_vec_unglued.begin(), nv_mtx.manager->inverse_renumbering_before_glue.begin() + nv_mtx.manager->inverse_renumbering_before_glue.size()),
                     hv.begin());
        cudaCheckError();
        thrust::fill( nv_vec_unglued.begin(), nv_vec_unglued.end(), 0.0 );
        thrust::copy(hv.begin(), hv.end(), nv_vec_unglued.begin());
    }
    else
    {
        // ASSUMPTION: when manager has not been allocated you are running on a single rank
        printf("Glue was called on a single rank\n");
    }

    return 0;
}

#if 0

// The folowing code is to perform an exhange halo using data that doesn't matches the topology of the matrix stored in its distributed manager
// We use instead other containers stored in the distributed manager. They are suffixed by "_before_glue"
// This allows to exchange vector halo between unglued vectors from glued matrices

template <class TConfig>
void exchange_halo_after_unglue(const Matrix<TConfig> &A, Vector<TConfig>  &data, int tag, int num_ring = 1)
{
    setup_after_unglue(data, A, num_ring);  //set pointers to buffer
    gather_B2L_after_unglue(A, data, num_ring);             //write values to buffer
    exchange_halo_after_unglue(data, A, num_ring); //exchange buffers
    //scatter_L2H(data);            //NULL op
}
/*
template <class TConfig>
void CommsMPIHostBufferStream<T_Config>::setup(DVector &b, const Matrix<TConfig> &m, int num_rings) { do_setup_after_unglue((b, m, num_rings);}
template <class T_Config>
void CommsMPIHostBufferStream<T_Config>::exchange_halo(DVector &b, const Matrix<TConfig> &m, cudaEvent_t event, int tag, int num_rings) {  do_exchange_halo_after_unglue((b, m, num_rings);}
template <class T_Config>
void CommsMPIHostBufferStream<T_Config>::setup(FVector &b, const Matrix<TConfig> &m, int tag, int num_rings) { do_setup_after_unglue((b, m, num_rings);}
template <class T_Config>
void CommsMPIHostBufferStream<T_Config>::exchange_halo(FVector &b, const Matrix<TConfig> &m, cudaEvent_t event, int tag, int num_rings) {  do_exchange_halo_after_unglue((b, m, num_rings);}
*/
template <class TConfig>
void setup_after_unglue(Vector<TConfig> &b, const Matrix<TConfig> &m, int num_rings)
{
    /*
    thrust::copy( m.manager->neighbors_before_glue.begin(),  m.manager->neighbors_before_glue.end(), std::ostream_iterator<int64_t>(std::cout, " "));
    thrust::copy( m.manager->halo_offsets_before_glue.begin(),  m.manager->halo_offsets_before_glue.end(), std::ostream_iterator<int64_t>(std::cout, " "));
    for (int i = 0; i < m.manager->B2L_rings_before_glue.size(); ++i)
    {
      thrust::copy( m.manager->B2L_rings_before_glue[i].begin(),  m.manager->B2L_rings_before_glue[i].end(), std::ostream_iterator<int64_t>(std::cout, " "));
    }

    for (int i = 0; i < m.manager->B2L_maps_before_glue.size(); ++i)
    {
      thrust::copy( m.manager->B2L_maps_before_glue[i].begin(),  m.manager->B2L_maps_before_glue[i].end(), std::ostream_iterator<int64_t>(std::cout, " "));
    }
    */
    if (TConfig::memSpace == AMGX_host)
    {
        FatalError("MPI Comms module no implemented for host", AMGX_ERR_NOT_IMPLEMENTED);
    }
    else
    {
#ifdef AMGX_WITH_MPI
        int bsize = b.get_block_size();
        int num_cols = b.get_num_cols();

        if (bsize != 1 && num_cols != 1)
            FatalError("Error: vector cannot have block size and subspace size.",
                       AMGX_ERR_INTERNAL);

        // set num neighbors = size of B2L_rings_before_glue
        // need to do this because comms might have more neighbors than our matrix knows about
        int neighbors = m.manager->B2L_rings_before_glue.size();
        m.manager->getComms()->set_neighbors(m.manager->B2L_rings_before_glue.size());

        if (b.in_transfer & SENDING)
        {
            b.in_transfer = IDLE;
        }

        typedef typename TConfig::template setVecPrec<(AMGX_VecPrecision)AMGX_GET_MODE_VAL(AMGX_MatPrecision, TConfig::mode)>::Type value_type;
        b.requests.resize(2 * neighbors); //first part is sends second is receives
        b.statuses.resize(2 * neighbors);

        for (int i = 0; i < 2 * neighbors; i++)
        {
            b.requests[i] = MPI_REQUEST_NULL;
        }

        int total = 0;

        for (int i = 0; i < neighbors; i++)
        {
            total += m.manager->B2L_rings_before_glue[i][num_rings] * bsize * num_cols;
        }

        b.buffer_size = total;

        if (b.buffer == NULL)
        {
            b.buffer = new Vector<TConfig>(total);
        }
        else
        {
            if (total > b.buffer->size())
            {
                b.buffer->resize(total);
            }
        }

        if (b.linear_buffers_size < neighbors)
        {
            if (b.linear_buffers_size != 0) { amgx::memory::cudaFreeHost(b.linear_buffers); }

            amgx::memory::cudaMallocHost((void **) & (b.linear_buffers), neighbors * sizeof(value_type *));
            b.linear_buffers_size = neighbors;
        }

        cudaCheckError();
        total = 0;
        bool linear_buffers_changed = false;

        for (int i = 0; i < neighbors; i++)
        {
            if (b.linear_buffers[i] != b.buffer->raw() + total)
            {
                linear_buffers_changed = true;
            }

            b.linear_buffers[i] = b.buffer->raw() + total;
            total += m.manager->B2L_rings_before_glue[i][num_rings] * bsize * num_cols;
        }

        // Copy to device
        if (linear_buffers_changed)
        {
            b.linear_buffers_ptrs.resize(neighbors);
            //thrust::copy(b.linear_buffers.begin(),b.linear_buffers.end(),b.linear_buffers_ptrs.begin());
            cudaMemcpyAsync(thrust::raw_pointer_cast(&b.linear_buffers_ptrs[0]), &(b.linear_buffers[0]), neighbors * sizeof(value_type *), cudaMemcpyHostToDevice);
            cudaCheckError();
        }

        int size = 0;
        size  = total + (m.manager->halo_offsets_before_glue[num_rings * neighbors] - m.manager->halo_offsets_before_glue[0]) * bsize * num_cols;

        if (size > 0)
        {
            if (b.explicit_host_buffer == NULL)
            {
                b.host_buffer.resize(1);
                cudaEventCreateWithFlags(&b.mpi_event, cudaEventDisableTiming);
                cudaCheckError();
                amgx::memory::cudaMallocHost((void **)&b.explicit_host_buffer, size * sizeof(value_type));
                cudaCheckError();
            }
            else if (size > b.explicit_buffer_size)
            {
                amgx::memory::cudaFreeHost(b.explicit_host_buffer);
                cudaCheckError();
                amgx::memory::cudaMallocHost((void **)&b.explicit_host_buffer, size * sizeof(value_type));
                cudaCheckError();
            }

            cudaCheckError();
            b.explicit_buffer_size = size;
        }

#else
        FatalError("MPI Comms module requires compiling with MPI", AMGX_ERR_NOT_IMPLEMENTED);
#endif
    }
}
template <class TConfig>
void gather_B2L_after_unglue(const Matrix<TConfig> &m, Vector<TConfig> &b, int num_rings = 1)
{
    if (TConfig::memSpace == AMGX_host)
    {
        if (m.manager->neighbors_before_glue.size() > 0)
        {
            FatalError("Distributed solve only supported on devices", AMGX_ERR_NOT_IMPLEMENTED);
        }
    }
    else
    {
        for (int i = 0; i < m.manager->neighbors_before_glue.size(); i++)
        {
            int size = m.manager->B2L_rings_before_glue[i][num_rings];
            int num_blocks = min(4096, (size + 127) / 128);

            if ( size != 0)
            {
                if (b.get_num_cols() == 1)
                {
                    gatherToBuffer <<< num_blocks, 128>>>(b.raw(), m.manager->B2L_maps_before_glue[i].raw(), b.linear_buffers[i], b.get_block_size(), size);
                }
                else
                {
                    gatherToBufferMultivector <<< num_blocks, 128>>>(b.raw(), m.manager->B2L_maps_before_glue[i].raw(), b.linear_buffers[i], b.get_num_cols(), b.get_lda(), size);
                }

                cudaCheckError();
            }
        }
    }
}
template <class TConfig>
void exchange_halo_after_unglue(Vector<TConfig> &b, const Matrix<TConfig> &m, int num_rings)
{
    if (TConfig::memSpace == AMGX_host)
    {
        FatalError("Halo exchanges not implemented for host", AMGX_ERR_NOT_IMPLEMENTED);
    }
    else
    {
#ifdef AMGX_WITH_MPI
        typedef typename TConfig::VecPrec VecPrec;
        cudaCheckError();
        int bsize = b.get_block_size();
        int num_cols = b.get_num_cols();
        int offset = 0;
        int neighbors = m.manager->B2L_rings_before_glue.size();
        MPI_Comm mpi_comm = m.manager->getComms()->get_mpi_comm();

        if (b.buffer_size != 0)
        {
            cudaMemcpy(&(b.explicit_host_buffer[0]), b.buffer->raw(), b.buffer_size * sizeof(typename TConfig::VecPrec), cudaMemcpyDeviceToHost);
        }

        for (int i = 0; i < neighbors; i++)
        {
            int size = m.manager->B2L_rings_before_glue[i][num_rings] * bsize * num_cols;

            if (size != 0)
            {
                MPI_Isend(&(b.explicit_host_buffer[offset]), size * sizeof(typename TConfig::VecPrec), MPI_BYTE, m.manager->neighbors_before_glue[i], m.manager->global_id(), mpi_comm, &b.requests[i]);
            }
            else
            {
                MPI_Isend(&(b.host_buffer[0]), size * sizeof(typename TConfig::VecPrec), MPI_BYTE, m.manager->neighbors_before_glue[i], m.manager->global_id(), mpi_comm, &b.requests[i]);
            }

            offset += size;
        }

        b.in_transfer = RECEIVING | SENDING;
        offset = 0;

        for (int i = 0; i < neighbors; i++)
        {
            // Count total size to receive from one neighbor
            int size = 0;

            for (int j = 0; j < num_rings; j++)
            {
                size += m.manager->halo_offsets_before_glue[j * neighbors + i + 1] * bsize * num_cols - m.manager->halo_offsets_before_glue[j * neighbors + i] * bsize * num_cols;
            }

            if (size != 0)
            {
                MPI_Irecv(&(b.explicit_host_buffer[b.buffer_size + offset]), size * sizeof(typename TConfig::VecPrec), MPI_BYTE, m.manager->neighbors_before_glue[i], m.manager->neighbors_before_glue[i], mpi_comm, &b.requests[neighbors + i]);
            }
            else
            {
                MPI_Irecv(&(b.host_buffer[0]), size * sizeof(typename TConfig::VecPrec), MPI_BYTE, m.manager->neighbors_before_glue[i], m.manager->neighbors_before_glue[i], mpi_comm, &b.requests[neighbors + i]);
            }

            offset += size;
            int required_size = m.manager->halo_offsets_before_glue[0] * bsize * num_cols + offset;

            if (required_size > b.size())
            {
                // happen because we have 2 ring
                // required_size correspond to "n" in the FULL view of the unconsolidated matrix.
                // In exchange halo this is a fatal error since it should never happen.
                b.resize(required_size);
            }
        }

        MPI_Waitall(2 * neighbors, &b.requests[0], /*&b.statuses[0]*/ MPI_STATUSES_IGNORE); //I only wait to receive data, I can start working before all my buffers were sent
        b.dirtybit = 0;
        b.in_transfer = IDLE;

        // copy on host ring by ring
        if (num_rings == 1)
        {
            if (num_cols == 1)
            {
                if (offset != 0)
                {
                    cudaMemcpy(b.raw() + m.manager->halo_offsets_before_glue[0]*bsize, &(b.explicit_host_buffer[b.buffer_size]), offset * sizeof(typename TConfig::VecPrec), cudaMemcpyHostToDevice);
                }
            }
            else
            {
                int lda = b.get_lda();
                VecPrec *rank_start = &(b.explicit_host_buffer[b.buffer_size]);

                for (int i = 0; i < neighbors; ++i)
                {
                    int halo_size = m.manager->halo_offsets_before_glue[i + 1] - m.manager->halo_offsets_before_glue[i];

                    for (int s = 0; s < num_cols; ++s)
                    {
                        VecPrec *halo_start = b.raw() + lda * s + m.manager->halo_offsets_before_glue[i];
                        VecPrec *received_halo = rank_start + s * halo_size;
                        cudaMemcpy(halo_start, received_halo, halo_size * sizeof(VecPrec), cudaMemcpyHostToDevice);
                    }

                    rank_start += num_cols * halo_size;
                }
            }
        }
        else
        {
            if (num_cols == 1)
            {
                offset = 0;

                // Copy into b, one neighbor at a time, one ring at a time
                for (int i = 0 ; i < neighbors ; i++)
                {
                    for (int j = 0; j < num_rings; j++)
                    {
                        int size = m.manager->halo_offsets_before_glue[j * neighbors + i + 1] * bsize - m.manager->halo_offsets_before_glue[j * neighbors + i] * bsize;

                        if (size != 0)
                        {
                            cudaMemcpy(b.raw() + m.manager->halo_offsets_before_glue[j * neighbors + i]*bsize, &(b.explicit_host_buffer[b.buffer_size + offset]), size * sizeof(typename TConfig::VecPrec), cudaMemcpyHostToDevice);
                        }

                        offset += size;
                    }
                }
            }
            else
            {
                FatalError("num_rings != 1 && num_cols != 1 not supported\n", AMGX_ERR_NOT_IMPLEMENTED);
            }
        }

#else
        FatalError("MPI Comms module requires compiling with MPI", AMGX_ERR_NOT_IMPLEMENTED);
#endif
    }
}
#endif
// if 0
#endif
//MPI
} // namespace amgx

max_stars_repo_name: fizmat/AMGX

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
#include <math.h>
#include <string.h>

#include <cblas.h>

#include "cnn.h"
#include "cnn_builtin_math.h"
#include "cnn_builtin_math_inline.h"
#include "cnn_config.h"
#include "cnn_init.h"
#include "cnn_macro.h"

#ifdef CNN_WITH_CUDA
#include <cuda_runtime.h>
#include "cnn_builtin_math_cu.h"
#endif

CNN_ACTIV_DEF((*cnn_activ_list[])) = {
    cnn_softmax,        //
    cnn_relu,           //
    cnn_swish,          //
    cnn_sigmoid,        //
    cnn_tanh,           //
    cnn_gaussian,       //
    cnn_bent_identity,  //
    cnn_softplus,       //
    cnn_softsign,       //
    cnn_sinc,           //
    cnn_sinusoid,       //
    cnn_identity        //
};

CNN_ACTIV_GRAD_DEF((*cnn_activ_grad_list[])) = {
    cnn_softmax_grad,        //
    cnn_relu_grad,           //
    cnn_swish_grad,          //
    cnn_sigmoid_grad,        //
    cnn_tanh_grad,           //
    cnn_gaussian_grad,       //
    cnn_bent_identity_grad,  //
    cnn_softplus_grad,       //
    cnn_softsign_grad,       //
    cnn_sinc_grad,           //
    cnn_sinusoid_grad,       //
    cnn_identity_grad        //
};

const char* cnn_activ_name[] = {
    "Softmax",             //
    "ReLU",                //
    "Swish",               //
    "Sigmoid",             //
    "Hyperbolic Tangent",  //
    "Gaussian",            //
    "Bent Identity",       //
    "SoftPlus",            //
    "SoftSign",            //
    "Sinc",                //
    "Sinusoid",            //
    "Identity"             //
};

CNN_ACTIV_DEF(cnn_softmax)
{
    float max, sum;

#ifdef CNN_WITH_CUDA
    // Find max value
    cnn_max_gpu(&max, src, len, buf);

    // Find shifted vector
    cnn_add_gpu(buf, src, len, -max);

    // Find exponential vector
    cnn_exp_gpu(buf, buf, len);

    // Find sum
    cnn_sum_gpu(&sum, buf, len, dst);

    // Find softmax
    cnn_div_gpu(dst, buf, len, sum);
#else
    int i;

    // Find max value
    max = src[0];
    for (i = 1; i < len; i++)
    {
        if (src[i] > max)
        {
            max = src[i];
        }
    }

    // Find exponential summation
    sum = 0;
    for (i = 0; i < len; i++)
    {
        dst[i] = src[i] - max;
        sum += exp(dst[i]);
    }

    // Find softmax output
    for (i = 0; i < len; i++)
    {
        dst[i] = exp(dst[i]) / sum;
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_softmax_grad)
{
#ifdef CNN_WITH_CUDA
    float alpha = 1.0;
    float beta = 0.0;

    // Find derivative matrix
    cnn_smax_grad_gpu(buf, cache, len);

    // Find layer gradient
    cnn_assert_cu(cublasSgemm(cnnInit.blasHandle, CUBLAS_OP_N, CUBLAS_OP_N,  //
                              len, 1, len,                                   //
                              &alpha,                                        //
                              buf, len,                                      //
                              gradIn, len,                                   //
                              &beta,                                         //
                              gradOut, len));

#else
    int i, j;

    // Find softmax gradient matrix
    for (i = 0; i < len; i++)
    {
        for (j = 0; j < len; j++)
        {
            buf[i * len + j] = cache[i] * ((float)(i == j) - cache[j]);
        }
    }

    // Find layer gradient
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,  //
                1, len, len,                                //
                1.0,                                        //
                gradIn, len,                                //
                buf, len,                                   //
                0.0,                                        //
                gradOut, len);
#endif
}

CNN_ACTIV_DEF(cnn_relu)
{
#ifdef CNN_WITH_CUDA
    cnn_relu_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_relu(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_relu_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_relu_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find relu gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_relu_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_swish)
{
#ifdef CNN_WITH_CUDA
    cnn_swish_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_swish(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_swish_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_swish_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find swish gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_swish_grad(gradOut + i, gradIn + i, src + i, cache + i);
    }
#endif
}

CNN_ACTIV_DEF(cnn_sigmoid)
{
#ifdef CNN_WITH_CUDA
    cnn_sigmoid_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_sigmoid(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_sigmoid_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_sigmoid_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find sigmoid gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_sigmoid_grad(gradOut + i, gradIn + i, NULL, cache + i);
    }
#endif
}

CNN_ACTIV_DEF(cnn_tanh)
{
#ifdef CNN_WITH_CUDA
    cnn_tanh_gpu(dst, src, len);
#else
#pragma omp parallel
    for (int i = 0; i < len; i++)
    {
        __cnn_tanh(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_tanh_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_tanh_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find tanh gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_tanh_grad(gradOut + i, gradIn + i, NULL, cache + i);
    }
#endif
}

CNN_ACTIV_DEF(cnn_gaussian)
{
#ifdef CNN_WITH_CUDA
    cnn_gaussian_gpu(dst, src, len);
#else
#pragma omp parallel
    for (int i = 0; i < len; i++)
    {
        __cnn_gaussian(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_gaussian_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_gaussian_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find gaussian gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_gaussian_grad(gradOut + i, gradIn + i, src + i, cache + i);
    }
#endif
}

CNN_ACTIV_DEF(cnn_bent_identity)
{
#ifdef CNN_WITH_CUDA
    cnn_bent_identity_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_bent_identity(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_bent_identity_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_bent_identity_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find bent indentity gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_bent_identity_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_softplus)
{
#ifdef CNN_WITH_CUDA
    cnn_softplus_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_softplus(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_softplus_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_softplus_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find softplus gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_softplus_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_softsign)
{
#ifdef CNN_WITH_CUDA
    cnn_softsign_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_softsign(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_softsign_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_softsign_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find softsign gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_softsign_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_sinc)
{
#ifdef CNN_WITH_CUDA
    cnn_sinc_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_sinc(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_sinc_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_sinc_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find sinc gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_sinc_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_sinusoid)
{
#ifdef CNN_WITH_CUDA
    cnn_sin_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_sinusoid(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_sinusoid_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_sin_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find sinusoid gradient
#pragma omp parallel for shared(gradOut, gradIn, src)
    for (int i = 0; i < len; i++)
    {
        __cnn_sinusoid_grad(gradOut + i, gradIn + i, src + i, NULL);
    }
#endif
}

CNN_ACTIV_DEF(cnn_identity)
{
#ifdef CNN_WITH_CUDA
    cnn_identity_gpu(dst, src, len);
#else
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_identity(dst + i, src + i);
    }
#endif
}

CNN_ACTIV_GRAD_DEF(cnn_identity_grad)
{
#ifdef CNN_WITH_CUDA
    cnn_identity_grad_gpu(gradOut, gradIn, src, len, cache);
#else
    // Find identity gradient
#pragma omp parallel for
    for (int i = 0; i < len; i++)
    {
        __cnn_identity_grad(gradOut + i, gradIn + i, NULL, NULL);
    }
#endif
}

int cnn_get_activ_id(const char* name)
{
    int i;
    int ret = CNN_PARSE_FAILED;

    if (name != NULL)
    {
        for (i = 0; i < CNN_ACTIV_AMOUNT; i++)
        {
            ret = strcmp(name, cnn_activ_name[i]);
            if (ret == 0)
            {
                ret = i;
                goto RET;
            }
        }
    }

RET:
    return ret;
}

ext: c

max_stars_repo_name: jamesljlster/cnn

max_issues_repo_name: jamesljlster/cnn

max_forks_repo_name: jamesljlster/cnn


text:
/* Odeint solver */
#include <pygsl/utils.h>
#include <pygsl/block_helpers.h>
#include <pygsl/error_helpers.h>
#include <Python.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

char odeiv_module_doc[] = "XXX odeiv module doc missing!\n";

static char this_file[] = __FILE__;
static PyObject *module = NULL; /* set by initodeiv */

#if 0
static void			/* generic instance destruction */
generic_dealloc (PyObject *self)
{
  DEBUG_MESS(1, " *** generic_dealloc %p\n", (void *) self);
  PyMem_Free(self);
}
#endif

typedef struct {
     PyObject_HEAD
     gsl_odeiv_step * step;
     gsl_odeiv_system system;
     PyObject *py_func;
     PyObject *py_jac;
     PyObject *arguments;
     jmp_buf buffer;
}
PyGSL_odeiv_step;

typedef struct {
     PyObject_HEAD
     PyGSL_odeiv_step * step;
     gsl_odeiv_control * control;
} PyGSL_odeiv_control;

typedef struct {
     PyObject_HEAD
     PyGSL_odeiv_step * step;
     PyGSL_odeiv_control * control;
     gsl_odeiv_evolve * evolve;
} PyGSL_odeiv_evolve;

/*---------------------------------------------------------------------------
 * Declaration of the various Methods
 *---------------------------------------------------------------------------*/
/*
 * stepper
 */
static int
PyGSL_odeiv_func(double t, const double y[], double f[], void *params);
static int
PyGSL_odeiv_jac(double t, const double y[], double *dfdy, double dfdt[],
		void *params);
static PyObject *
PyGSL_odeiv_step_apply(PyGSL_odeiv_step *self, PyObject *args);
static PyObject *
PyGSL_odeiv_step_reset(PyGSL_odeiv_step *self, PyObject *args);
static PyObject *
PyGSL_odeiv_step_name(PyGSL_odeiv_step *self, PyObject *args);
static PyObject *
PyGSL_odeiv_step_order(PyGSL_odeiv_step *self, PyObject *args);
static void
PyGSL_odeiv_step_free(PyGSL_odeiv_step * self);

/*
 * control
 */
static PyObject *
PyGSL_odeiv_control_hadjust(PyGSL_odeiv_control *self, PyObject *args);
static PyObject *
PyGSL_odeiv_control_name(PyGSL_odeiv_control *self, PyObject *args);
static void
PyGSL_odeiv_control_free(PyGSL_odeiv_control * self);

/*
 * evolve
 */
static PyObject *
PyGSL_odeiv_evolve_apply(PyGSL_odeiv_evolve *self, PyObject *args);
static PyObject *
PyGSL_odeiv_evolve_apply_array(PyGSL_odeiv_evolve *self, PyObject *args);
static PyObject *
PyGSL_odeiv_evolve_reset(PyGSL_odeiv_evolve *self, PyObject *args);
static void
PyGSL_odeiv_evolve_free(PyGSL_odeiv_evolve * self);
/*---------------------------------------------------------------------------*/

static char  PyGSL_odeiv_step_doc[] = "XXX Documentation missing\n";
static char  PyGSL_odeiv_control_doc[] = "XXX Documentation missing\n";
static char  PyGSL_odeiv_evolve_doc[] = "XXX Documentation missing\n";

static char odeiv_step_apply_doc[] =  "XXX Documentation missing\n";
static char odeiv_step_name_doc[] =  "XXX Documentation missing\n";
static char odeiv_step_order_doc[] =  "XXX Documentation missing\n";
static char odeiv_step_reset_doc[] =  "XXX Documentation missing\n";

static char odeiv_control_name_doc[] =  "XXX Documentation missing\n";
static char odeiv_control_hadjust_doc[] =  "XXX Documentation missing\n";

static char odeiv_evolve_apply_doc[] =  "XXX Documentation missing\n";
static char odeiv_evolve_apply_array_doc[] =  "XXX Documentation missing\n";
static char odeiv_evolve_reset_doc[] =  "XXX Documentation missing\n";

static char * odeiv_step_init_err_msg = "odeiv_step.__init__";

static struct PyMethodDef PyGSL_odeiv_step_methods[] = {
     {"apply", (PyCFunction) PyGSL_odeiv_step_apply, METH_VARARGS, odeiv_step_apply_doc},
     {"reset", (PyCFunction) PyGSL_odeiv_step_reset, METH_VARARGS, odeiv_step_reset_doc},
     {"name",  (PyCFunction) PyGSL_odeiv_step_name,  METH_VARARGS, odeiv_step_name_doc},
     {"order", (PyCFunction) PyGSL_odeiv_step_order, METH_VARARGS, odeiv_step_order_doc},
     {NULL, NULL}
};
static struct PyMethodDef PyGSL_odeiv_control_methods[] = {
     {"hadjust", (PyCFunction) PyGSL_odeiv_control_hadjust, METH_VARARGS, odeiv_control_hadjust_doc},
     {"name",  (PyCFunction) PyGSL_odeiv_control_name,  METH_VARARGS, odeiv_control_name_doc},
     {NULL, NULL}
};

static struct PyMethodDef PyGSL_odeiv_evolve_methods[] = {
     {"apply", (PyCFunction) PyGSL_odeiv_evolve_apply, METH_VARARGS, odeiv_evolve_apply_doc},
     {"reset", (PyCFunction) PyGSL_odeiv_evolve_reset, METH_VARARGS, odeiv_evolve_reset_doc},
     {"apply_array", (PyCFunction) PyGSL_odeiv_evolve_apply_array, METH_VARARGS, odeiv_evolve_apply_array_doc},
     {NULL, NULL}
};

#define PyGSL_ODEIV_GENERIC_ALL                                                           \
static PyObject *									  \
PyGSL_ODEIV_GENERIC_GETATTR(PyGSL_ODEIV_GENERIC *self, char *name)		          \
{											  \
     PyObject *tmp = NULL;								  \
											  \
     FUNC_MESS_BEGIN();									  \
     tmp = Py_FindMethod(PyGSL_ODEIV_GENERIC_METHODS, (PyObject *) self, name);	          \
     if(NULL == tmp){	  								  \
	  PyGSL_add_traceback(module, __FILE__, "odeiv.__attr__", __LINE__ - 1);	  \
	  return NULL;									  \
     }											  \
     FUNC_MESS_END();									  \
     return tmp;                                                                          \
}                                                                                \
static PyTypeObject PyGSL_ODEIV_GENERIC_PYTYPE = {					  \
  PyObject_HEAD_INIT(NULL)	/* fix up the type slot in initcrng */			  \
  0,				/* ob_size */						  \
  PyGSL_ODEIV_GENERIC_NAME,			/* tp_name */			          \
  sizeof(PyGSL_ODEIV_GENERIC),  /* tp_basicsize */					  \
  0,				/* tp_itemsize */					  \
											  \
  /* standard methods */								  \
  (destructor)  PyGSL_ODEIV_GENERIC_DELETE,       /* tp_dealloc  ref-count==0  */	  \
  (printfunc)   0,		                      /* tp_print    "print x"     */	  \
  (getattrfunc) PyGSL_ODEIV_GENERIC_GETATTR,       /* tp_getattr  "x.attr"      */	  \
  (setattrfunc) 0,		   /* tp_setattr  "x.attr=v"    */			  \
  (cmpfunc)     0,		   /* tp_compare  "x > y"       */			  \
  (reprfunc)    0,                 /* tp_repr     `x`, print x  */			  \
											  \
  /* type categories */									  \
  0,				/* tp_as_number   +,-,*,/,%,&,>>,pow...*/		  \
  0,				/* tp_as_sequence +,[i],[i:j],len, ...*/		  \
  0,				/* tp_as_mapping  [key], len, ...*/			  \
											  \
  /* more methods */									  \
  (hashfunc)     0,		/* tp_hash    "dict[x]" */				  \
  (ternaryfunc)  0,             /* tp_call    "x()"     */				  \
  (reprfunc)     0,             /* tp_str     "str(x)"  */				  \
  (getattrofunc) 0,		/* tp_getattro */					  \
  (setattrofunc) 0,		/* tp_setattro */					  \
  0,				/* tp_as_buffer */					  \
  0L,				/* tp_flags */						  \
  PyGSL_ODEIV_GENERIC_DOC       /* doc */                                                 \
};

#define PyGSL_ODEIV_GENERIC            PyGSL_odeiv_step
#define PyGSL_ODEIV_GENERIC_NAME      "PyGSL_odeiv_step"
#define PyGSL_ODEIV_GENERIC_PYTYPE     PyGSL_odeiv_step_pytype
#define PyGSL_ODEIV_GENERIC_DOC        PyGSL_odeiv_step_doc
#define PyGSL_ODEIV_GENERIC_GETATTR    PyGSL_odeiv_step_getattr
#define PyGSL_ODEIV_GENERIC_METHODS    PyGSL_odeiv_step_methods
#define PyGSL_ODEIV_GENERIC_DELETE     PyGSL_odeiv_step_free
PyGSL_ODEIV_GENERIC_ALL
/**/;
#undef PyGSL_ODEIV_GENERIC
#undef PyGSL_ODEIV_GENERIC_NAME
#undef PyGSL_ODEIV_GENERIC_PYTYPE
#undef PyGSL_ODEIV_GENERIC_DOC
#undef PyGSL_ODEIV_GENERIC_GETATTR
#undef PyGSL_ODEIV_GENERIC_METHODS
#undef PyGSL_ODEIV_GENERIC_DELETE
#define PyGSL_ODEIV_GENERIC            PyGSL_odeiv_control
#define PyGSL_ODEIV_GENERIC_NAME      "PyGSL_odeiv_control"
#define PyGSL_ODEIV_GENERIC_PYTYPE     PyGSL_odeiv_control_pytype
#define PyGSL_ODEIV_GENERIC_DOC        PyGSL_odeiv_control_doc
#define PyGSL_ODEIV_GENERIC_GETATTR    PyGSL_odeiv_control_getattr
#define PyGSL_ODEIV_GENERIC_METHODS    PyGSL_odeiv_control_methods
#define PyGSL_ODEIV_GENERIC_DELETE     PyGSL_odeiv_control_free
PyGSL_ODEIV_GENERIC_ALL
/**/;
#undef PyGSL_ODEIV_GENERIC
#undef PyGSL_ODEIV_GENERIC_NAME
#undef PyGSL_ODEIV_GENERIC_PYTYPE
#undef PyGSL_ODEIV_GENERIC_DOC
#undef PyGSL_ODEIV_GENERIC_GETATTR
#undef PyGSL_ODEIV_GENERIC_METHODS
#undef PyGSL_ODEIV_GENERIC_DELETE
#define PyGSL_ODEIV_GENERIC            PyGSL_odeiv_evolve
#define PyGSL_ODEIV_GENERIC_NAME      "PyGSL_odeiv_evolve"
#define PyGSL_ODEIV_GENERIC_PYTYPE     PyGSL_odeiv_evolve_pytype
#define PyGSL_ODEIV_GENERIC_DOC        PyGSL_odeiv_evolve_doc
#define PyGSL_ODEIV_GENERIC_GETATTR    PyGSL_odeiv_evolve_getattr
#define PyGSL_ODEIV_GENERIC_METHODS    PyGSL_odeiv_evolve_methods
#define PyGSL_ODEIV_GENERIC_DELETE     PyGSL_odeiv_evolve_free
PyGSL_ODEIV_GENERIC_ALL
/**/;

#define PyGSL_ODEIV_STEP_Check(v)    ((v)->ob_type == &PyGSL_odeiv_step_pytype)
#define PyGSL_ODEIV_CONTROL_Check(v) ((v)->ob_type == &PyGSL_odeiv_control_pytype)
#define PyGSL_ODEIV_EVOLVE_Check(v)  ((v)->ob_type == &PyGSL_odeiv_evolve_pytype)

/*---------------------------------------------------------------------------
  Wrapper functions to push call the approbriate Python Objects
  ---------------------------------------------------------------------------*/
static int
PyGSL_odeiv_func(double t, const double y[], double f[], void *params)
{
    int dimension, flag = GSL_FAILURE;

    PyObject *arglist = NULL, *result = NULL;
    PyArrayObject *yo = NULL;
    PyGSL_odeiv_step * step;
    gsl_vector_view yv, fv;
    PyGSL_error_info  info;

    FUNC_MESS_BEGIN();

    step = (PyGSL_odeiv_step *) params;
    if(!PyGSL_ODEIV_STEP_Check(step)){
	  PyGSL_add_traceback(module, this_file, __FUNCTION__,
			      __LINE__ - 2);
	  gsl_error("Param not a step type!",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
    }
    dimension =  step->system.dimension;

    /* Do I need to copy the array ??? */
    yv = gsl_vector_view_array((double *) y, dimension);
    yo = PyGSL_copy_gslvector_to_pyarray(&yv.vector);
    if (yo == NULL) goto fail;

    FUNC_MESS("\t\tBuild args");
    arglist = Py_BuildValue("(dOO)", t, yo, step->arguments);
    FUNC_MESS("\t\tEnd Build args");

    info.callback = step->py_func;
    info.message  = "odeiv_func";
    result  = PyEval_CallObject(step->py_func, arglist);

    if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 1, &info)) != GSL_SUCCESS){
	  goto fail;
     }
    info.argnum = 1;
    fv = gsl_vector_view_array(f, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslvector(&fv.vector, result, dimension,
					       &info)) != GSL_SUCCESS){
	  goto fail;
     }

    Py_DECREF(arglist);    arglist = NULL;
    Py_DECREF(yo);         yo = NULL;
    Py_DECREF(result);     result = NULL;
    FUNC_MESS_END();
    return GSL_SUCCESS;

 fail:
    FUNC_MESS("    IN Fail BEGIN");
    Py_XDECREF(yo);
    Py_XDECREF(result);
    Py_XDECREF(arglist);
    FUNC_MESS("    IN Fail END");
    assert(flag != GSL_SUCCESS);
    longjmp(step->buffer, flag);
    return flag;
}

static int
PyGSL_odeiv_jac(double t, const double y[], double *dfdy, double dfdt[],
		void *params)
{
    int dimension, flag = GSL_FAILURE;
    PyGSL_odeiv_step *step = NULL;
    PyGSL_error_info  info;

    PyObject *arglist = NULL, *result = NULL, *tmp=NULL;
    PyArrayObject *yo = NULL;

    gsl_vector_view yv, dfdtv;
    gsl_matrix_view dfdyv;

    FUNC_MESS_BEGIN();

    step = (PyGSL_odeiv_step *) params;
    if(!PyGSL_ODEIV_STEP_Check(step)){
	  PyGSL_add_traceback(module, this_file, __FUNCTION__,
			      __LINE__ - 2);
	  gsl_error("Param not a step type!",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
    }
    dimension = step->system.dimension;

    yv = gsl_vector_view_array((double *) y, dimension);
    yo = PyGSL_copy_gslvector_to_pyarray(&yv.vector);
    if (yo == NULL) goto fail;

    arglist = Py_BuildValue("(dOO)", t, yo, step->arguments);

    assert(step->py_jac);
    result  = PyEval_CallObject(step->py_jac, arglist);

    info.callback = step->py_jac;
    info.message  = "odeiv_jac";
    if((flag = PyGSL_CHECK_PYTHON_RETURN(result, 2, &info)) != GSL_SUCCESS){
	  goto fail;
     }

    info.argnum = 1;
    tmp = PyTuple_GET_ITEM(result, 0);
    dfdyv = gsl_matrix_view_array((double *) dfdy, dimension, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslmatrix(&dfdyv.matrix, tmp, dimension, dimension, &info)) != GSL_SUCCESS){
	  goto fail;
    }

    info.argnum = 2;
    tmp = PyTuple_GET_ITEM(result, 1);
    dfdtv = gsl_vector_view_array((double *) dfdt, dimension);
    if((flag = PyGSL_copy_pyarray_to_gslvector(&dfdtv.vector, tmp, dimension, &info)) != GSL_SUCCESS){
	  goto fail;
    }

    Py_DECREF(arglist);    arglist = NULL;
    Py_DECREF(result);     result = NULL;
    Py_DECREF(yo);         yo = NULL;
    FUNC_MESS_END();
    return GSL_SUCCESS;

 fail:
    FUNC_MESS("IN Fail");
    assert(flag != GSL_SUCCESS);
    longjmp(step->buffer, flag);
    return flag;
}

/* Wrappers for the evaluation of the system */
static PyObject *
PyGSL_odeiv_step_apply(PyGSL_odeiv_step *self, PyObject *args)
{
    PyObject *result = NULL;
    PyObject *y0_o = NULL, *dydt_in_o = NULL;
    PyArrayObject *volatile y0 = NULL, * volatile yerr = NULL,
	 *volatile dydt_in = NULL, *volatile dydt_out = NULL,
	 *volatile yout = NULL;

    double t=0, h=0, *volatile dydt_in_d;
    int dimension, r, flag;

    FUNC_MESS_BEGIN();
    assert(PyGSL_ODEIV_STEP_Check(self));
    if(! PyArg_ParseTuple(args, "ddOOO", &t, &h, &y0_o,  &dydt_in_o)){
      return NULL;
    }

    dimension = self->system.dimension;
    y0 = PyGSL_PyArray_PREPARE_gsl_vector_view(y0_o, PyArray_DOUBLE, 1, dimension, 1, NULL);
    if(y0 == NULL) goto fail;

    if (Py_None == dydt_in_o){
	 dydt_in_d = NULL;
    }else{
	 dydt_in = PyGSL_PyArray_PREPARE_gsl_vector_view(dydt_in_o, PyArray_DOUBLE, 1, dimension, 2, NULL);
	 if(dydt_in == NULL) goto fail;
	 dydt_in_d = (double *) dydt_in->data;
    }

    dydt_out = (PyArrayObject *)  PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
    if (dydt_out == NULL) goto fail;

    yerr = (PyArrayObject *) PyArray_FromDims(1, &dimension, PyArray_DOUBLE);
    if(yerr == NULL) goto fail;

    yout = (PyArrayObject *) PyArray_CopyFromObject((PyObject * ) y0, PyArray_DOUBLE, 1, 1);
    if(yout == NULL) goto fail;

    if((flag=setjmp(self->buffer)) == 0){
	  FUNC_MESS("\t\t Setting Jmp Buffer");
    } else {
	  FUNC_MESS("\t\t Returning from Jmp Buffer");
	  goto fail;
    }

    r = gsl_odeiv_step_apply(self->step, t, h,
			     (double *) yout->data,
			     (double *) yerr->data,
			     dydt_in_d,
			     (double *) dydt_out->data,
			     &(self->system));
    if (GSL_SUCCESS != r){
	PyErr_SetString(PyExc_TypeError, "Error While evaluating gsl_odeiv");
      goto fail;
    }

    FUNC_MESS("    Returnlist create ");
    assert(yout != NULL);
    assert(yerr != NULL);
    assert(dydt_out != NULL);

    result = Py_BuildValue("(OOO)", yout, yerr, dydt_out);

    FUNC_MESS("    Memory free ");
    /* Deleting the arrays */
    Py_DECREF(y0);           y0 = NULL;
    Py_DECREF(yout);         yout = NULL;
    Py_DECREF(yerr);         yerr = NULL;
    Py_DECREF(dydt_out);     dydt_out = NULL;
    /* This array does not need to exist ... */
    Py_XDECREF(dydt_in);	 dydt_in=NULL;

    FUNC_MESS_END();
    return result;

    fail:
    FUNC_MESS("IN Fail");
    Py_XDECREF(y0);
    Py_XDECREF(yout);
    Py_XDECREF(yerr);
    Py_XDECREF(dydt_in);
    Py_XDECREF(dydt_out);
    FUNC_MESS("IN Fail End");
    return NULL;
}

static void
PyGSL_odeiv_step_free(PyGSL_odeiv_step * self)
{
     assert(PyGSL_ODEIV_STEP_Check(self));
     Py_DECREF(self->py_func);
     Py_XDECREF(self->py_jac);
     Py_DECREF(self->arguments);
     gsl_odeiv_step_free(self->step);
     PyMem_Free(self);
}

static PyObject *
PyGSL_odeiv_step_reset(PyGSL_odeiv_step *self, PyObject *args)
{
     assert(PyGSL_ODEIV_STEP_Check(self));
     gsl_odeiv_step_reset(self->step);
     Py_INCREF(Py_None);
     return Py_None;
}

static PyObject *
PyGSL_odeiv_step_name(PyGSL_odeiv_step *self, PyObject *args)
{
     assert(PyGSL_ODEIV_STEP_Check(self));
     return PyString_FromString(gsl_odeiv_step_name(self->step));
}

static PyObject *
PyGSL_odeiv_step_order(PyGSL_odeiv_step *self, PyObject *args)
{
     assert(PyGSL_ODEIV_STEP_Check(self));
     return PyInt_FromLong((long) gsl_odeiv_step_order(self->step));
}

/* --------------------------------------------------------------------------- */
/* control_hadjust needs a few arrays */
/*
extern int
gsl_odeiv_control_hadjust (gsl_odeiv_control * c, gsl_odeiv_step * s,
			   const double y0[],  const double yerr[],
			   const double dydt[], double * h);
*/
static PyObject *
PyGSL_odeiv_control_hadjust(PyGSL_odeiv_control *self, PyObject *args)
{

  PyObject *result = NULL;
  PyObject *y0_o = NULL, *yerr_o = NULL, *dydt_o = NULL;
  PyArrayObject *y0 = NULL, *yerr = NULL, *dydt = NULL;
  double h = 0;
  int r = 0;

  size_t dimension = 0;

  FUNC_MESS_BEGIN();
  assert(PyGSL_ODEIV_CONTROL_Check(self));
  if(!PyArg_ParseTuple(args, "OOOd",  &y0_o, &yerr_o, &dydt_o, &h)){
    return NULL;
  }

  dimension = self->step->system.dimension;

  y0 = PyGSL_PyArray_PREPARE_gsl_vector_view(y0_o, PyArray_DOUBLE, 1, dimension,  1, NULL);
  if(y0 == NULL)   goto fail;
  yerr = PyGSL_PyArray_PREPARE_gsl_vector_view(yerr_o, PyArray_DOUBLE, 1, dimension, 2, NULL);
  if(yerr == NULL) goto fail;
  dydt = PyGSL_PyArray_PREPARE_gsl_vector_view(yerr_o, PyArray_DOUBLE, 1, dimension, 3, NULL);
  if(dydt == NULL) goto fail;

  FUNC_MESS("      Array Pointers End");

  r = gsl_odeiv_control_hadjust(self->control, self->step->step,
				(double *) y0->data,
				(double *) yerr->data,
				(double *) dydt->data, &h);

  FUNC_MESS("      Function End");
  Py_DECREF(y0);       y0 = NULL;
  Py_DECREF(yerr);     yerr = NULL;
  Py_DECREF(dydt);     dydt = NULL;

  result = Py_BuildValue("di",h,r);
  FUNC_MESS_END();
  return result;

 fail:
  FUNC_MESS("IN Fail");
  Py_XDECREF(y0);
  Py_XDECREF(yerr);
  Py_XDECREF(dydt);
  FUNC_MESS("IN Fail END");
  return NULL;
}

static void
PyGSL_odeiv_control_free(PyGSL_odeiv_control * self)
{
     assert(PyGSL_ODEIV_CONTROL_Check(self));
     Py_DECREF(self->step);
     //gsl_odeiv_control_free(self->control);
     PyMem_Free(self);
}

static PyObject *
PyGSL_odeiv_control_name(PyGSL_odeiv_control *self, PyObject *args)
{
     assert(PyGSL_ODEIV_CONTROL_Check(self));
     return PyString_FromString(gsl_odeiv_control_name(self->control));
}

static void
PyGSL_odeiv_evolve_free(PyGSL_odeiv_evolve * self)
{
     assert(PyGSL_ODEIV_EVOLVE_Check(self));
     Py_DECREF(self->step);
     Py_DECREF(self->control);
     gsl_odeiv_evolve_free(self->evolve);
     PyMem_Free(self);
}

static PyObject *
PyGSL_odeiv_evolve_apply(PyGSL_odeiv_evolve *self, PyObject *args)
{
    PyObject *result = NULL;
    PyObject *y0_o = NULL;
    PyArrayObject *volatile y0 = NULL, *volatile yout = NULL;

    double t=0, h=0, t1 = 0, flag;

    int dimension, r;

    assert(PyGSL_ODEIV_EVOLVE_Check(self));
    FUNC_MESS_BEGIN();

    if(! PyArg_ParseTuple(args, "dddO",
			  &t, &t1, &h, &y0_o)){
      return NULL;
    }

    dimension = self->step->system.dimension;

    y0 = PyGSL_PyArray_PREPARE_gsl_vector_view(y0_o, PyArray_DOUBLE, 1, dimension, 1, NULL);
    if(y0 == NULL) goto fail;

    yout = (PyArrayObject *)  PyArray_CopyFromObject((PyObject * ) y0, PyArray_DOUBLE, 1, 1);
    if(yout == NULL) goto fail;

    if((flag=setjmp(self->step->buffer)) == 0){
	  FUNC_MESS("\t\t Setting Jmp Buffer");
     } else {
	  FUNC_MESS("\t\t Returning from Jmp Buffer");
	  goto fail;
     }

    r = gsl_odeiv_evolve_apply(self->evolve,
			       self->control->control,
			       self->step->step,
			       &(self->step->system), &t, t1, &h,
			       (double * )yout->data);

    if (GSL_SUCCESS != r){
	 goto fail;
    }

    assert(yout != NULL);

    result = Py_BuildValue("(ddO)", t, h, yout);

    /* Deleting the arrays */
    /* Do I need to do that??? Not transfered Py_DECREF(yout);     yout = NULL; */
    Py_DECREF(y0);       y0=NULL;
    FUNC_MESS_END();
    return result;

 fail:
    FUNC_MESS("IN Fail");
    PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__);
    Py_XDECREF(y0);
    Py_XDECREF(yout);
    FUNC_MESS("IN Fail End");
    return NULL;
}

static PyObject *
PyGSL_odeiv_evolve_apply_array(PyGSL_odeiv_evolve *self, PyObject *args)
{
    PyObject *result = NULL;
    PyObject *y0_o = NULL;
    PyObject *t_o = NULL;
    PyArrayObject *volatile y0 = NULL, *volatile yout = NULL, *volatile ta = NULL;

    double t=0, h=0, t1 = 0, flag;
    double *y0_data = NULL, *yout_data = NULL;
    int dimension, r, dims[2];
    int nlen=-1, i, j;

    assert(PyGSL_ODEIV_EVOLVE_Check(self));
    FUNC_MESS_BEGIN();

    if(! PyArg_ParseTuple(args, "OdO", &t_o, &h, &y0_o)){
      return NULL;
    }

    dimension = self->step->system.dimension;

    ta = PyGSL_PyArray_PREPARE_gsl_vector_view(t_o, PyArray_DOUBLE, 1, -1, 1, NULL);
    if(ta == NULL) goto fail;
    nlen = ta->dimensions[0];

    y0 = PyGSL_PyArray_PREPARE_gsl_vector_view(y0_o, PyArray_DOUBLE, 1, dimension, 1, NULL);
    if(y0 == NULL) goto fail;

    dims[0] = nlen;
    dims[1] = dimension;
    yout = (PyArrayObject *)  PyArray_FromDims(2, dims, PyArray_DOUBLE);
    if(yout == NULL) goto fail;

    if((flag=setjmp(self->step->buffer)) == 0){
	  FUNC_MESS("\t\t Setting Jmp Buffer");
     } else {
	  FUNC_MESS("\t\t Returning from Jmp Buffer");
	  goto fail;
     }

    DEBUG_MESS(5, "\t\t Yout Layout %p, shape=[%d,%d], strides=[%d,%d]",
	       yout->data, yout->dimensions[0], yout->dimensions[1], yout->strides[0], yout->strides[1]);
    r = GSL_CONTINUE;

    yout_data = (double *)(yout->data);
    /* Copy the start vector */
    for(j=0; j<dimension; j++){
	 yout_data[j] = *((double *)(y0->data + y0->strides[0] * j));
    }
    for(i=1; i<nlen; i++){
	 y0_data   = (double *)(yout->data + yout->strides[0]*(i-1));
	 yout_data = (double *)(yout->data + yout->strides[0]*(i)  );
	 /* Copy the start vector */
	 for(j=0; j<dimension; j++){
	      yout_data[j] = y0_data[j];
	 }
	 t =  *((double *) (ta->data + ta->strides[0]*(i-1)) );
	 t1 = *((double *) (ta->data + ta->strides[0]*(i)  ) );

	 while(t<t1){
	      r = gsl_odeiv_evolve_apply(self->evolve,
					 self->control->control,
					 self->step->step,
					 &(self->step->system), &t, t1, &h,
					 yout_data);
	      /* All GSL Errors are > 0. */
	      if (r != GSL_SUCCESS){
		   goto fail;
	      }
	 }
	 assert(r == GSL_SUCCESS);
    }

    result = (PyObject *) yout;

    /* Deleting the arrays */
    Py_DECREF(y0); y0=NULL;
    Py_DECREF(ta);
    FUNC_MESS_END();
    return result;

 fail:
    FUNC_MESS("IN Fail");
    PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__);
    Py_XDECREF(ta);  ta = NULL;
    Py_XDECREF(y0);  y0=NULL;
    Py_XDECREF(yout);
    FUNC_MESS("IN Fail End");
    return NULL;
}

static PyObject *
PyGSL_odeiv_evolve_reset(PyGSL_odeiv_evolve *self, PyObject *args)
{
     assert(PyGSL_ODEIV_EVOLVE_Check(self));
     gsl_odeiv_evolve_reset(self->evolve);
     Py_INCREF(Py_None);
     return Py_None;
}

static PyObject *
PyGSL_odeiv_step_init(PyObject *self, PyObject *args, PyObject *kwdict, const gsl_odeiv_step_type * odeiv_type)
{

     PyObject *func=NULL, *jac=NULL, *o_params=NULL;
     PyGSL_odeiv_step *odeiv_step = NULL;

     static char * kwlist[] = {"dimension", "func", "jac", "args", NULL};
     int dim, has_jacobian = 0;

     FUNC_MESS_BEGIN();
     assert(args);
     if (0 == PyArg_ParseTupleAndKeywords(args, kwdict, "iOOO:odeiv_step.__init__", kwlist,
					  &dim, &func, &jac, &o_params)){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 2);
	  return NULL;
     }
     if (dim <= 0){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	  gsl_error("The dimension of the problem must be at least 1",
		    this_file, __LINE__ -2, GSL_EDOM);
	  return NULL;
     }
     if(!PyCallable_Check(func)){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	  gsl_error("The function object is not callable!",
		    this_file, __LINE__ -2, GSL_EBADFUNC);
	  goto fail;
     }

     if(jac == Py_None){
	  if(odeiv_type == gsl_odeiv_step_bsimp){
	       PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	       gsl_error("The bsimp method needs a jacobian! You supplied None.",
			 this_file, __LINE__ -2, GSL_EBADFUNC);
	       goto fail;
	  }
     }else{
	  if(!PyCallable_Check(jac)){
	       PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg, __LINE__ - 1);
	       gsl_error("The jacobian object must be None or callable!",
			 this_file, __LINE__ -2, GSL_EBADFUNC);
	       goto fail;
	  }
	  has_jacobian = 1;

     }

     odeiv_step =  (PyGSL_odeiv_step *) PyObject_NEW(PyGSL_odeiv_step, &PyGSL_odeiv_step_pytype);
     if(odeiv_step == NULL){
	  PyErr_NoMemory();
	  goto fail;
     }
     odeiv_step->step=NULL;
     odeiv_step->py_func=NULL;
     odeiv_step->py_jac=NULL;
     odeiv_step->arguments=NULL;

     odeiv_step->system.dimension = dim;
     if(has_jacobian)
	  odeiv_step->system.jacobian = PyGSL_odeiv_jac;
     else
	  odeiv_step->system.jacobian = NULL;

     odeiv_step->system.function = PyGSL_odeiv_func;
     odeiv_step->system.params = (void *) odeiv_step;

     odeiv_step->step =  gsl_odeiv_step_alloc(odeiv_type, dim);
     if(odeiv_step->step == NULL){
	  Py_DECREF(odeiv_step);
	  PyErr_NoMemory();
	  return NULL;
     }

     odeiv_step->py_func=func;
     if(has_jacobian)
	  odeiv_step->py_jac=jac;

     odeiv_step->arguments = o_params;
     Py_INCREF(odeiv_step->py_func);
     Py_INCREF(odeiv_step->arguments);

     Py_XINCREF(odeiv_step->py_jac);

     FUNC_MESS_END();
     return (PyObject *) odeiv_step;
 fail:
     return NULL;
}
#define ADD_ODESTEPPER(mytype)                                                      \
static PyObject *                                                                   \
PyGSL_odeiv_step_init_ ## mytype (PyObject * self, PyObject * args, PyObject *kwdic)\
{                                                                                \
     return PyGSL_odeiv_step_init(self, args, kwdic, gsl_odeiv_step_ ## mytype);    \
}
ADD_ODESTEPPER(rk2)
ADD_ODESTEPPER(rk4)
ADD_ODESTEPPER(rkf45)
ADD_ODESTEPPER(rkck)
ADD_ODESTEPPER(rk8pd)
ADD_ODESTEPPER(rk2imp)
ADD_ODESTEPPER(rk4imp)
ADD_ODESTEPPER(bsimp)
ADD_ODESTEPPER(gear1)
ADD_ODESTEPPER(gear2)

static PyObject *
PyGSL_odeiv_control_init(PyObject *self, PyObject *args, void * type)
{
     int nargs = -1, tmp=0;
     double eps_abs, eps_rel, a_y, a_dydt;
     PyGSL_odeiv_step *step=NULL;
     PyGSL_odeiv_control *a_con=NULL;

     gsl_odeiv_control *(*evaluator_5)(double , double , double , double ) = NULL;
     gsl_odeiv_control *(*evaluator_3)(double , double ) = NULL;

     FUNC_MESS_BEGIN();
     /* The arguments depend on the type of control */
     if(type == (void *) gsl_odeiv_control_standard_new){
	  /* step, eps_abs, eps_rel, a_y, a_dydt */
	  nargs = 5;
     }else if(type == (void *) gsl_odeiv_control_y_new ||
	      type == (void *) gsl_odeiv_control_yp_new){
	  nargs = 3;
     }else{
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg,
			      __LINE__ - 2);
	  gsl_error("Unknown control type",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail;
     }
     assert(nargs > -1);

     switch(nargs){
     case 5:
	  tmp == PyArg_ParseTuple(args, "O!dddd:odeiv_control.__init__",
				  &PyGSL_odeiv_step_pytype, &step, &eps_abs, &eps_rel, &a_y, &a_dydt);
	  break;
     case 3:
	  tmp == PyArg_ParseTuple(args, "O!dd:odeiv_control.__init__",
				  &PyGSL_odeiv_step_pytype, &step, &eps_abs, &eps_rel);
	  break;
     default:
	  fprintf(stderr, "nargs = %d\n", nargs);
	  gsl_error("Unknown number of arguments",
		    this_file, __LINE__ -2, GSL_EFAULT);
	  goto fail; break;
     }

     if(!step){
	  PyErr_SetString(PyExc_TypeError, "The first argument must be a step object!");
	  goto fail;
     }

     if(tmp){
	  PyGSL_add_traceback(module, this_file, odeiv_step_init_err_msg,
			      __LINE__ - 2);
	  return NULL;
     }

     a_con =  PyObject_NEW(PyGSL_odeiv_control, &PyGSL_odeiv_control_pytype);
     if (NULL == a_con){
	  PyErr_NoMemory();
	  goto fail;
     }
     a_con->control=NULL;
     switch(nargs){
     case 5:
	  evaluator_5 = type;
	  a_con->control = evaluator_5(eps_abs, eps_rel, a_y, a_dydt);
	  break;
     case 3:
	  evaluator_3 = type;
	  a_con->control = evaluator_3(eps_abs, eps_rel);
	  break;
     default:
	  goto fail;
     }
     if (NULL == a_con->control){
	  PyErr_NoMemory();
	  goto fail;
     }
     assert(step);
     a_con->step =  step;
     Py_INCREF(step);
     FUNC_MESS_END();
     return (PyObject *) a_con;

 fail:
     FUNC_MESS("FAIL");
     Py_XDECREF(a_con);
     return NULL;

}

#define ADD_ODECONTROL(name)                                                  \
static PyObject *                                                             \
PyGSL_odeiv_control_init_ ## name (PyObject * self, PyObject * args)          \
{                                                                             \
     return PyGSL_odeiv_control_init(self, args, (void *) gsl_odeiv_control_ ## name);    \
}
ADD_ODECONTROL(standard_new)
ADD_ODECONTROL(y_new)
ADD_ODECONTROL(yp_new)

static PyObject *
PyGSL_odeiv_evolve_init(PyObject *self, PyObject *args)
{
     PyGSL_odeiv_step *step = NULL;
     PyGSL_odeiv_control *control = NULL;
     PyGSL_odeiv_evolve *a_ev = NULL;

     /* step, control */
     FUNC_MESS_BEGIN();
     if(0== PyArg_ParseTuple(args, "O!O!:odeiv_evolve.__init__",
			     &PyGSL_odeiv_step_pytype, &step,
			     &PyGSL_odeiv_control_pytype, &control)){
	  return NULL;
     }

     if(!step){
	  PyErr_SetString(PyExc_TypeError, "The first argument must be a step object!");
	  goto fail;
     }

     if(!control){
	  PyErr_SetString(PyExc_TypeError, "The second argument must be a control object!");
	  goto fail;
     }

     a_ev =  (PyGSL_odeiv_evolve *) PyObject_NEW(PyGSL_odeiv_evolve, &PyGSL_odeiv_evolve_pytype);
     if(NULL == a_ev){
	  PyErr_NoMemory();
	  return NULL;
     }

     a_ev->step = step;
     a_ev->control = control;

     a_ev->evolve = gsl_odeiv_evolve_alloc(step->system.dimension);
     if(NULL == a_ev){
	  PyErr_NoMemory();
	  goto fail;
     }
     Py_INCREF(step);
     Py_INCREF(control);
     FUNC_MESS_END();
     return (PyObject *) a_ev;
 fail:
     FUNC_MESS("FAIL");
     Py_DECREF(a_ev);
     return NULL;
}

static PyMethodDef PyGSL_odeiv_module_functions[] = {
     {"step_rk2",    PyGSL_odeiv_step_init_rk2,    METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rk4",    PyGSL_odeiv_step_init_rk4,    METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rkf45",  PyGSL_odeiv_step_init_rkf45,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rkck",   PyGSL_odeiv_step_init_rkck,   METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rk8pd",  PyGSL_odeiv_step_init_rk8pd,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rk2imp", PyGSL_odeiv_step_init_rk2imp, METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_rk4imp", PyGSL_odeiv_step_init_rk4imp, METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_bsimp",  PyGSL_odeiv_step_init_bsimp,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_gear1",  PyGSL_odeiv_step_init_gear1,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"step_gear2",  PyGSL_odeiv_step_init_gear2,  METH_VARARGS|METH_KEYWORDS, NULL},
     {"control_standard_new", PyGSL_odeiv_control_init_standard_new, METH_VARARGS, NULL},
     {"control_y_new",        PyGSL_odeiv_control_init_y_new,        METH_VARARGS, NULL},
     {"control_yp_new",       PyGSL_odeiv_control_init_yp_new,       METH_VARARGS, NULL},
     {"evolve", PyGSL_odeiv_evolve_init, METH_VARARGS, NULL},
     {NULL, NULL, 0}        /* Sentinel */
};

void
initodeiv(void)
{
     PyObject *m=NULL, *item=NULL, *dict=NULL;

     FUNC_MESS_BEGIN();
     fprintf(stderr, "Compiled at %s %s\n", __DATE__, __TIME__);
     m = Py_InitModule("odeiv", PyGSL_odeiv_module_functions);
     assert(m);
     module = m;
     import_array();
     init_pygsl();

     PyGSL_odeiv_step_pytype.ob_type = &PyType_Type;
     PyGSL_odeiv_control_pytype.ob_type = &PyType_Type;
     PyGSL_odeiv_evolve_pytype.ob_type = &PyType_Type;

     dict = PyModule_GetDict(m);
     /* create_odeiv_step_types(dict); */
     if(!dict)
	  goto fail;

     if (!(item = PyString_FromString(odeiv_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }
     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();
     return;
 fail:
     FUNC_MESS("Fail");
     fprintf(stderr, "Import of module odeiv failed!\n");
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


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

#pragma once

#include "arcana/functional/inplace_function.h"
#include "arcana/sentry.h"

#include "affinity.h"
#include "blocking_concurrent_queue.h"

#include <gsl/gsl>
#include <vector>

namespace mira
{
    template<size_t WorkSize>
    class dispatcher
    {
    public:
        using callback_t = stdext::inplace_function<void(), WorkSize>;
        static constexpr size_t work_size = WorkSize;

        template<typename T>
        void queue(T&& work)
        {
            m_work.push(std::forward<T>(work));
        }

        affinity get_affinity() const
        {
            return m_affinity;
        }

        dispatcher(const dispatcher&) = delete;
        dispatcher& operator=(const dispatcher&) = delete;

        virtual ~dispatcher() = default;

    protected:
        dispatcher() = default;
        dispatcher& operator=(dispatcher&&) = default;
        dispatcher(dispatcher&&) = default;

        bool tick(const cancellation& token)
        {
            return internal_tick(token, false);
        }

        bool blocking_tick(const cancellation& token)
        {
            return internal_tick(token, true);
        }

        /*
        Sets the dispatcher's tick thread affinity. Once this is set the methods
        on this instance will need to be called by that thread.
        */
        void set_affinity(const affinity& aff)
        {
            m_affinity = aff;
        }

        void cancelled()
        {
            m_work.cancelled();
        }

        void clear()
        {
            m_work.clear();
        }

    private:
        bool internal_tick(const cancellation& token, bool block)
        {
            GSL_CONTRACT_CHECK("thread affinity", m_affinity.check());

            if (block)
            {
                if (!m_work.blocking_drain(m_workload, token))
                    return false;
            }
            else
            {
                if (!m_work.try_drain(m_workload, token))
                    return false;
            }

            for (auto& work : m_workload)
            {
                work();
            }

            m_workload.clear();

            return true;
        }

        blocking_concurrent_queue<callback_t> m_work;
        affinity m_affinity;
        std::vector<callback_t> m_workload;
    };

    template<size_t WorkSize>
    class manual_dispatcher : public dispatcher<WorkSize>
    {
    public:
        using dispatcher<WorkSize>::blocking_tick;
        using dispatcher<WorkSize>::cancelled;
        using dispatcher<WorkSize>::clear;
        using dispatcher<WorkSize>::set_affinity;
        using dispatcher<WorkSize>::tick;
    };

    template<size_t WorkSize>
    class background_dispatcher : public dispatcher<WorkSize>
    {
    public:
        background_dispatcher()
            : m_registration{ m_cancellation.add_listener([this] { this->cancelled(); }) }
        {
            m_thread = std::thread{ [&]() {

                // TODO: Set the affinity when usage bugs are fixed.
                constexpr bool should_set_affinity{ false };
                if (should_set_affinity)
                {
                    this->set_affinity(std::this_thread::get_id());
                }

                while (!m_cancellation.cancelled())
                {
                    this->blocking_tick(m_cancellation);
                }
            } };
        }

        void cancel()
        {
            m_cancellation.cancel();

            if (m_thread.joinable())
            {
                m_thread.join();
            }

            this->clear();
        }

        ~background_dispatcher()
        {
            cancel();
        }

    private:
        std::thread m_thread;
        cancellation_source m_cancellation;
        cancellation_source::ticket m_registration;
    };
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
/** @file */
#ifndef __CCL_UTILS_H_INCLUDED__
#define __CCL_UTILS_H_INCLUDED__

#include <gsl/gsl_spline.h>

#define CCL_MIN(a, b)  (((a) < (b)) ? (a) : (b))
#define CCL_MAX(a, b)  (((a) > (b)) ? (a) : (b))

CCL_BEGIN_DECLS

typedef enum ccl_integration_t {
  ccl_integration_qag_quad = 500,  // GSL's quad
  ccl_integration_spline = 501,  // Spline integral
} ccl_integration_t;

/**
 * Compute bin edges of N-1 linearly spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_linear_spacing(double xmin, double xmax, int N);

/**
 * Compute bin edges of N-1 logarithmically and then linearly spaced bins on the interval [xmin,xmax]
 * @param xminlog minimum value of spacing
 * @param xmin value when logarithmical ends and linear spacing begins
 * @param xmax maximum value of spacing
 * @param Nlin number of linear bins plus one (number of bin edges)
 * @param Nlog number of logarithmic bins plus one (number of bin edges)
 * @return x, bin edges in range [xminlog, xmax]
 */
double * ccl_linlog_spacing(double xminlog, double xmin, double xmax, int Nlin, int Nlog);

/**
 * Compute bin edges of N-1 logarithmically spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_log_spacing(double xmin, double xmax, int N);
//Returns array of N logarithmically-spaced values between xmin and xmax

double ccl_j_bessel(int l,double x);
//Spherical Bessel function of order l (adapted from CAMB)

/**
 * Compute spline integral.
 * @param nx number of elements in input array.
 * @param ny number of y arrays.
 * @param x input x-values.
 * @param y input y-values (ny arrays with nx elements).
 * @param a lower end of integration range.
 * @param b upper end of integration range (use b<a if you just want to integrate over all of y).
 * @param result array of output spline integral values.
 * @param T spline type.
 * @param status status flag.
 */
void ccl_integ_spline(int ny, int nx,double *x,double **y,
                      double a, double b, double *result,
                      const gsl_interp_type *T, int *status);

CCL_END_DECLS

#endif

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL