text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef rbtree_ebcaa429_edae_4ac6_b9a5_2fe774778822_h
#define rbtree_ebcaa429_edae_4ac6_b9a5_2fe774778822_h

#include <assert.h>
#include <gslib/std.h>

__gslib_begin__

struct _rbtree_trait_copy {};
struct _rbtree_trait_detach {};

enum rbtree_color
{
    rb_red,
    rb_black,
};

template<class _ty>
struct _rbtreenode_cpy_wrapper
{
    typedef _ty value;
    typedef _rbtreenode_cpy_wrapper<_ty> myref;
    typedef _rbtree_trait_copy tsf_behavior;
    typedef rbtree_color color;

    value               _value;
    myref*              _left;
    myref*              _right;
    myref*              _parent;
    color               _color;

    myref()
    {
        _left = _right = _parent = nullptr;
        _color = rb_black;
    }
    value* get_ptr() { return &_value; }
    const value* const_ptr() const { return &_value; }
    value& get_ref() { return _value; }
    const value& const_ref() const { return _value; }
    void born() {}
    void kill() {}
    template<class _ctor>
    void born() {}
    template<class _ctor>
    void kill() {}
    void copy(const myref* a) { get_ref() = a->const_ref(); }
    void attach(myref* a) { assert(0); }
    void swap_data(myref* a) { std::swap(_value, a->_value); }
    void set_color(color c) { _color = c; }
    color get_color() const { return _color; }
    bool is_red() const { return _color == rb_red; }
};

template<class _ty>
struct _rbtreenode_wrapper
{
    typedef _ty value;
    typedef _rbtreenode_wrapper<_ty> myref;
    typedef _rbtree_trait_detach tsf_behavior;
    typedef rbtree_color color;

    value*              _value;
    myref*              _left;
    myref*              _right;
    myref*              _parent;
    color               _color;

    myref()
    {
        _left = _right = _parent = nullptr;
        _value = nullptr;
        _color = rb_black;
    }
    value* get_ptr() { return _value; }
    const value* const_ptr() const { return _value; }
    value& get_ref() { return *_value; }
    const value& const_ref() const { return *_value; }
    void copy(const myref* a) { get_ref() = a->const_ref(); }
    void born() { !! }
    template<class _ctor>
    void born() { _value = new _ctor; }
    void kill() { if(_value) { delete _value; _value = nullptr; } }
    template<class _ctor>
    void kill() { if(_value) { delete _value; _value = nullptr; } }
    void attach(myref* a)
    {
        assert(a && a->_value);
        kill();
        _value = a->_value;
        a->_value = nullptr;
    }
    void swap_data(myref* a) { gs_swap(_value, a->_value); }
    void set_color(color c) { _color = c; }
    color get_color() const { return _color; }
    bool is_red() const { return _color == rb_red; }
};

template<class _wrapper>
struct _rbtree_allocator
{
    typedef _wrapper wrapper;
    static wrapper* born() { return new wrapper; }
    static void kill(wrapper* w) { delete w; }
};

template<class _val>
struct _rbtreenode_val
{
    typedef _val value;
    typedef const _val const_value;
    typedef _rbtreenode_val<_val> myref;
    typedef rbtree_color color;

    union
    {
        value*          _vptr;
        const_value*    _cvptr;
    };

    myref() { _vptr = nullptr; }
    value* get_wrapper() const { return _vptr; }
    operator bool() const { return _vptr != nullptr; }
    bool is_left() const { return (_cvptr && _cvptr->_parent) ? _cvptr->_parent->_left == _cvptr : false; }
    bool is_right() const { return (_cvptr && _cvptr->_parent) ? _cvptr->_parent->_right == _cvptr : false; }
    bool is_root() const { return _cvptr ? (!_cvptr->_parent) : false; }
    bool is_leaf() const { return _cvptr ? (!_cvptr->_left && !_cvptr->_right) : false; }
    int up_depth() const
    {
        int depth = 0;
        for(value* p = _vptr; p; p = p->_parent, depth ++);
        return depth;
    }
    int down_depth() const { return _down_depth(_vptr, 0); }
    bool operator==(const value* v) const { return _vptr == v; }
    bool operator!=(const value* v) const { return _vptr != v; }
    color get_color() const { return _vptr->_color; }
    void set_color(color c) { _vptr->_color = c; }
    void swap_data(myref& a) { _vptr->swap_data(a._vptr); }

public:
    static void connect_left_child(value* p, value* l)
    {
        assert(p);
        p->_left = l;
        if(l)
            l->_parent = p;
    }
    static void connect_right_child(value* p, value* r)
    {
        assert(p);
        p->_right = r;
        if(r)
            r->_parent = p;
    }
    static bool disconnect_parent_child(value* p, value* c)
    {
        assert(p && c && (c->_parent == p));
        if(p->_left == c) {
            p->_left = c->_parent = nullptr;
            return true;
        }
        assert(p->_right == c);
        p->_right = c->_parent = nullptr;
        return false;
    }

private:
    static int _down_depth(value* v, int ctr)
    {
        if(v == nullptr)
            return ctr;
        ctr ++;
        return gs_max(_down_depth(v->_left, ctr),
            _down_depth(v->_right, ctr)
            );
    }

protected:
    value* vleft() const { return _vptr ? _vptr->_left : nullptr; }
    value* vright() const { return _vptr ? _vptr->_right : nullptr; }
    value* vparent() const { return _vptr ? _vptr->_parent : nullptr; }
    value* vsibling() const
    {
        if(!_vptr || !_vptr->_parent)
            return nullptr;
        if(is_left())
            return _vptr->_right;
        else if(is_right())
            return _vptr->_left;
        assert(!"unexpected.");
        return nullptr;
    }
    value* vroot() const
    {
        if(!_vptr)
            return nullptr;
        value* p = _vptr;
        for( ; p->_parent; p = p->_parent);
        return p;
    }

protected:
    template<class _lambda, class _value>
    static void preorder_traversal(_lambda lam, _value* v)
    {
        assert(v);
        lam(v);
        if(v->_left)
            preorder_traversal(lam, v->_left);
        if(v->_right)
            preorder_traversal(lam, v->_right);
    }
    template<class _lambda, class _value>
    static void inorder_traversal(_lambda lam, _value* v)
    {
        assert(v);
        if(v->_left)
            inorder_traversal(lam, v->_left);
        lam(v);
        if(v->_right)
            inorder_traversal(lam, v->_right);
    }
    template<class _lambda, class _value>
    static void postorder_traversal(_lambda lam, _value* v)
    {
        assert(v);
        if(v->_left)
            postorder_traversal(lam, v->_left);
        if(v->_right)
            postorder_traversal(lam, v->-right);
        lam(v);
    }

public:
    template<class _lambda>
    void inorder_traversal(_lambda lam) { if(_vptr) inorder_traversal(lam, _vptr); }
    template<class _lambda>
    void inorder_traversal(_lambda lam) const { if(_cvptr) inorder_traversal(lam, _cvptr); }
    template<class _lambda>
    void preorder_traversal(_lambda lam) { if(_vptr) preorder_traversal(lam, _vptr); }
    template<class _lambda>
    void preorder_traversal(_lambda lam) const { if(_cvptr) preorder_traversal(lam, _cvptr); }
    template<class _lambda>
    void postorder_traversal(_lambda lam) { if(_vptr) postorder_traversal(lam, _vptr); }
    template<class _lambda>
    void postorder_traversal(_lambda lam) const { if(_cvptr) postorder_traversal(lam, _cvptr); }
};

template<class _ty,
    class _wrapper = _rbtreenode_cpy_wrapper<_ty> >
class _rbtree_const_iterator:
    public _rbtreenode_val<_wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _rbtree_const_iterator<_ty, _wrapper> iterator;

public:
    iterator(const wrapper* w = nullptr) { _cvptr = w; }
    bool is_valid() const { return _cvptr != nullptr; }
    const value* get_ptr() const { return _cvptr->const_ptr(); }
    const value* operator->() const { return _cvptr->const_ptr(); }
    const value& operator*() const { return _cvptr->const_ref(); }
    iterator left() const { return iterator(vleft()); }
    iterator right() const { return iterator(vright()); }
    iterator parent() const { return iterator(vparent()); }
    iterator sibling() const { return iterator(vsibling()); }
    iterator root() const { return iterator(vroot()); }
    bool operator==(const iterator& that) const { return _cvptr == that._cvptr; }
    bool operator!=(const iterator& that) const { return _cvptr != that._cvptr; }
};

template<class _ty,
    class _wrapper = _rbtreenode_cpy_wrapper<_ty> >
class _rbtree_iterator:
    public _rbtree_const_iterator<_ty, _wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _rbtree_const_iterator<_ty, _wrapper> const_iterator;
    typedef _rbtree_const_iterator<_ty, _wrapper> superref;
    typedef _rbtree_iterator<_ty, _wrapper> iterator;

public:
    iterator(wrapper* w): superref(w) {}
    value* get_ptr() const { return _vptr->get_ptr(); }
    value* operator->() const { return _vptr->get_ptr(); }
    value& operator*() const { return _vptr->get_ref(); }
    bool operator==(const iterator& that) const { return _vptr == that._vptr; }
    bool operator!=(const iterator& that) const { return _vptr != that._vptr; }
    bool operator==(const const_iterator& that) const { return _vptr == that._vptr; }
    bool operator!=(const const_iterator& that) const { return _vptr != that._vptr; }
    operator const_iterator() { return const_iterator(_cvptr); }
    void to_root() { _vptr = vroot(); }
    void to_left() { _vptr = vleft(); }
    void to_right() { _vptr = vright(); }
    void to_sibling() { _vptr = vsibling(); }
    void to_parent() { _vptr = vparent(); }
    iterator left() const { return iterator(vleft()); }
    iterator right() const { return iterator(vright()); }
    iterator parent() const { return iterator(vparent()); }
    iterator sibling() const { return iterator(vsibling()); }
    iterator root() const { return iterator(vroot()); }
};

template<class _ty,
    class _wrapper = _rbtreenode_cpy_wrapper<_ty>,
    class _alloc = _rbtree_allocator<_wrapper> >
class rbtree:
    public _rbtreenode_val<_wrapper>
{
public:
    typedef _ty value;
    typedef _wrapper wrapper;
    typedef _alloc alloc;
    typedef rbtree<value, wrapper, alloc> myref;
    typedef _rbtree_const_iterator<_ty, _wrapper> const_iterator;
    typedef _rbtree_iterator<_ty, _wrapper> iterator;
    typedef rbtree_color color;

public:
    rbtree() { _vptr = nullptr; }
    ~rbtree() { clear(); }
    void clear() { destroy(get_root()); }
    void destroy(iterator i)
    {
        if(!i.is_valid())
            return;
        if(iterator p = i.parent())
            disconnect_parent_child(p.get_wrapper(), i.get_wrapper());
        else {
            assert(is_root(i));
            _vptr = nullptr;
        }
        _destroy(i);
    }
    void adopt(wrapper* w)
    {
        assert(!_vptr && "use attach method.");
        _vptr = w;
    }
    iterator get_root() const { return iterator(_vptr); }
    const_iterator const_root() const { return const_iterator(_cvptr); }
    bool is_root(iterator i) const { return i ? (_vptr == i.get_wrapper()) : false; }
    bool is_valid() const { return _cvptr != nullptr; }
    bool is_mine(iterator i) const
    {
        if(!i.is_valid())
            return false;
        i.to_root();
        return i.get_wrapper() == _vptr;
    }
    int depth() const { return _cvptr->down_depth(); }
    void swap(myref& that) { gs_swap(_vptr, that._vptr); }
    iterator find(const value& v) const { return find(get_root(), v); }
    iterator find(iterator i, const value& v) const
    {
        if(!i) {
            if(!is_valid())
                return i;
            i = get_root();
        }
        assert(i);
        wrapper* f = _find(i.get_wrapper(), v);
        assert(f);
        if(!(f->const_ref() == v))
            f = nullptr;
        return iterator(f);
    }
    template<class _ctor = value>
    iterator insert(const value& v)
    {
        wrapper* f = nullptr;
        if(_vptr) {
            f = _find(_vptr, v);
            assert(f);
            if(f->const_ref() == v)
                return iterator(f);
        }
        wrapper* n = initval<_ctor, wrapper::tsf_behavior>::run(alloc::born(), v);
        n->set_color(rb_red);
        if(!f)
            _vptr = n;
        else {
            (v < f->const_ref()) ? connect_left_child(f, n) :
                connect_right_child(f, n);
        }
        _fix_insert(n);
        return iterator(n);
    }
    void erase(iterator i)
    {
        assert(i);
        wrapper* n = i.get_wrapper();
        wrapper *c, *p;
        rbtree_color cr;
        if(!n->_left)
            c = n->_right;
        else if(!n->_right)
            c = n->_left;
        else {
            wrapper* b = n, *l;
            n = n->_right;
            for(; l = n->_left; n = l);
            if(!b->_parent)
                _vptr = n;
            else {
                (b->_parent->_left == b) ? b->_parent->_left = n :
                    b->_parent->_right = n;
            }
            c = n->_right;
            p = n->_parent;
            cr = n->_color;
            if(p == b)
                p = n;
            else {
                if(c)
                    c->_parent = p;
                p->_left = c;
                n->_right = b->_right;
                b->_right->_parent = n;
            }
            n->_parent = b->_parent;
            n->_color = b->_color;
            n->_left = b->_left;
            b->_left->_parent = n;
            if(cr == rb_black)
                _fix_erase(c, p);
            return;
        }
        p = n->_parent;
        cr = n->_color;
        if(c)
            c->_parent = p;
        if(!p)
            _vptr = c;
        else {
            (p->_left == n) ? p->_left = c :
                p->_right = c;
        }
        if(cr == rb_black)
            _fix_erase(c, p);
    }
    void erase(const value& v)
    {
        if(iterator i = find(v))
            erase(i);
    }

    /* The detach and attach methods, provide subtree operations */
    myref& detach(myref& subtree, iterator i)
    {
        assert(i && is_mine(i));
        if(subtree.is_valid())
            subtree.clear();
        detach<myref>(subtree, i);
        return subtree;
    }
    template<class _cont>
    void detach(_cont& cont)
    {
        cont.adopt(_vptr);
        _vptr = nullptr;
    }
    template<class _cont>
    void detach(_cont& cont, iterator i)
    {
        assert(i && is_mine(i));
        if(i == get_root())
            return detach(cont);
        iterator p = i.parent();
        assert(p);
        disconnect_parent_child(p.get_wrapper(), i.get_wrapper());
        cont.adopt(i.get_wrapper());
    }
    iterator attach(myref& subtree, iterator i)
    {
        assert(i && is_mine(i) && i.is_leaf());
        if(i.is_root()) {
            swap(subtree);
            return get_root();
        }
        iterator p = i.parent();
        assert(p);
        bool leftp = disconnect_parent_child(p.get_wrapper(), i.get_wrapper());
        gs_swap(subtree._vptr, i._vptr);
        leftp ? connect_left_child(p.get_wrapper(), i.get_wrapper()) :
            connect_right_child(p.get_wrapper(), i.get_wrapper());
        subtree.clear();
        return i;
    }

public:
    template<class _lambda>
    void preorder_for_each(_lambda lam) { preorder_traversal([](wrapper* w) { lam(w->get_ptr()); }); }
    template<class _lambda>
    void preorder_const_for_each(_lambda lam) const { preorder_traversal([](const wrapper* w) { lam(w->const_ptr()); }); }
    template<class _lambda>
    void inorder_for_each(_lambda lam) { inorder_traversal([](wrapper* w) { lam(w->get_ptr()); }); }
    template<class _lambda>
    void inorder_const_for_each(_lambda lam) const { inorder_traversal([](const wrapper* w) { lam(w->const_ptr()); }); }
    template<class _lambda>
    void postorder_for_each(_lambda lam) { postorder_traversal([](wrapper* w) { lam(w->get_ptr()); }); }
    template<class _lambda>
    void postorder_const_for_each(_lambda lam) const { postorder_traversal([](const wrapper* w) { lam(w->const_ptr()); }); }

protected:
    void _destroy(iterator i)
    {
        if(!i.is_valid())
            return;
        _destroy(i.left());
        _destroy(i.right());
        wrapper* w = i.get_wrapper();
        w->kill();
        alloc::kill(w);
    }
    wrapper* _find(wrapper* w, const value& v) const
    {
        assert(w);
        wrapper* n = w;
        for(;;) {
            wrapper* m = nullptr;
            const value& r = n->const_ref();
            if(v < r)
                m = n->_left;
            else if(r < v)
                m = n->_right;
            if(!m)
                return n;
            n = m;
        }
        return nullptr;
    }
    void _left_rotate(wrapper* n)
    {
        assert(n);
        wrapper* r = n->_right;
        if(n->_right = r->_left)
            r->_left->_parent = n;
        r->_left = n;
        if(!(r->_parent = n->_parent))
            _vptr = r;
        else {
            (n == n->_parent->_right) ? n->_parent->_right = r :
                n->_parent->_left = r;
        }
        n->_parent = r;
    }
    void _right_rotate(wrapper* n)
    {
        assert(n);
        wrapper* l = n->_left;
        if(n->_left = l->_right)
            l->_right->_parent = n;
        l->_right = n;
        if(!(l->_parent = n->_parent))
            _vptr = l;
        else {
            (n == n->_parent->_right) ? n->_parent->_right = l :
                n->_parent->_left = l;
        }
        n->_parent = l;
    }
    void _fix_insert(wrapper* n)
    {
        assert(n);
        wrapper *p, *g, *u;
        while((p = n->_parent) && p->is_red()) {
            g = p->_parent;
            if(p == g->_left) {
                u = g->_right;
                if(u && u->is_red()) {
                    u->set_color(rb_black);
                    p->set_color(rb_black);
                    g->set_color(rb_red);
                    n = g;
                }
                else {
                    if(p->_right == n) {
                        _left_rotate(p);
                        gs_swap(p, n);
                    }
                    p->set_color(rb_black);
                    g->set_color(rb_red);
                    _right_rotate(g);
                }
            }
            else {
                u = g->_left;
                if(u && u->is_red()) {
                    u->set_color(rb_black);
                    p->set_color(rb_black);
                    g->set_color(rb_red);
                    n = g;
                }
                else {
                    if(p->_left == n) {
                        _right_rotate(p);
                        gs_swap(p, n);
                    }
                    p->set_color(rb_black);
                    g->set_color(rb_red);
                    _left_rotate(g);
                }
            }
        }
        _vptr->set_color(rb_black);
    }
    void _fix_erase(wrapper* n, wrapper* p)
    {
        wrapper *s, *sl, *sr;
        while((!n || !n->is_red()) && n != _vptr) {
            if(p->_left == n) {
                s = p->_right;
                if(s->is_red()) {
                    s->set_color(rb_black);
                    p->set_color(rb_red);
                    _left_rotate(p);
                    s = p->_right;
                }
                if((!s->_left || !s->_left->is_red()) &&
                    (!s->_right || !s->_right->is_red())
                    ) {
                    s->set_color(rb_red);
                    n = p;
                    p = n->_parent;
                }
                else {
                    if(!s->_right || !s->_right->is_red()) {
                        if(sl = s->_left)
                            sl->set_color(rb_black);
                        s->set_color(rb_red);
                        _right_rotate(s);
                        s = p->_right;
                    }
                    s->set_color(p->get_color());
                    p->set_color(rb_black);
                    if(sr = s->_right)
                        sr->set_color(rb_black);
                    _left_rotate(p);
                    n = _vptr;
                    break;
                }
            }
            else {
                s = p->_left;
                if(s->is_red()) {
                    s->set_color(rb_black);
                    p->set_color(rb_red);
                    _right_rotate(p);
                    s = p->_left;
                }
                if((!s->_left || !s->_left->is_red()) &&
                    (!s->_right || !s->_right->is_red())
                    ) {
                    s->set_color(rb_red);
                    n = p;
                    p = n->_parent;
                }
                else {
                    if(!s->_left || !s->_left->is_red()) {
                        if(sr = s->_right)
                            sr->set_color(rb_black);
                        s->set_color(rb_red);
                        _left_rotate(s);
                        s = p->_left;
                    }
                    s->set_color(p->get_color());
                    p->set_color(rb_black);
                    if(sl = s->_left)
                        sl->set_color(rb_black);
                    _right_rotate(p);
                    n = _vptr;
                    break;
                }
            }
        }
        if(n)
            n->set_color(rb_black);
    }

protected:
    friend struct initval;
    template<class _ctor, class _tsftrait>
    struct initval;
    template<class _ctor>
    struct initval<_ctor, _rbtree_trait_copy>
    {
        static wrapper* run(wrapper* w, const value& v)
        {
            assert(w);
            w->get_ref() = v;
            return w;
        }
    };
    template<class _ctor>
    struct initval<_ctor, _rbtree_trait_detach>
    {
        static wrapper* run(wrapper* w, const value& v)
        {
            assert(w);
            w->_value = &v;     /* or duplicate? */
            return w;
        }
    };
};

__gslib_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef _PYGSL_SOLVER_H_
#define _PYGSL_SOLVER_H_ 1
#include <pygsl/intern.h>
#include <pygsl/block_helpers.h>
/* Not directly needed here, but provides a lot of convienience functions */
#include <pygsl/error_helpers.h>
#include <gsl/gsl_math.h>
#include <setjmp.h>
/*
 * Many functions are "just" accessor methods. These different methods are
 * listed here.
 *
 * Convention: they all end with "m_t" short for method type. This emphasises
 * that this type is pointer to a method.
 */
typedef int (*int_m_t)(void *);
typedef size_t (*size_t_m_t)(void *);
typedef void (*void_m_t)(void *);
typedef void *(*void_a_t)(const void *);
typedef void *(*void_an_t)(const void *, size_t n);
typedef void *(*void_anp_t)(const void *, size_t n, size_t p);
typedef const char * (*name_m_t)(void *);
typedef double  (*double_m_t)(void *);
typedef gsl_vector *(*ret_vec)(void *);
typedef int (*set_m_t)(void *, void *, const gsl_vector *);
typedef int (*set_m_d_t) (void *, gsl_function *, double);
typedef int (*set_m_ddd_t) (void *, gsl_function *, double, double, double);
typedef int (*int_f_vd_t)(const gsl_vector *, double);
typedef int (*int_f_vvdd_t)(const gsl_vector *, const gsl_vector *, double, double);

struct _PyGSLSolverObject;

/*
 * GSL Methods which are implemented as C functions. The specifiy pointer to
 * the solver struct are passed as (void *) pointers.
 */
struct _GSLMethods{
     /* Method to be called to free the GSL solver */
     void_m_t free;
     /*
      * Some solvers provide a restart method. If not available set it to NULL.
      */
     void_m_t restart;
     /* returns a string with the name of the solver */
     name_m_t name;
     /* takes one more step towards the solution */
     int_m_t iterate;
};

struct _SolverStatic{
     struct _GSLMethods cmethods;
     /* How many callbacks will be used? */
     int n_cbs;
     /* Additional methods not provided by the basis solver. */
     PyMethodDef *pymethods;
     /* Describes the type of the solver. e. g. F-Minimizer */
     const char * type_name;
};

struct pygsl_array_cache{
     double * data;
     PyArrayObject * ref;
};
#define PyGSL_SOLVER_NCBS_MAX 4
#define PyGSL_SOLVER_PB_ND_MAX 2
#define PyGSL_SOLVER_N_ARRAYS 10

struct _PyGSLSolverObject{
     PyObject_HEAD
     /*
      *	Some solvers do not propagate errors. Here I use longjmp if an error
      * is raised by the evaluator.
      */
     jmp_buf buffer;
     /*
      * Array generation is a vital part of the callback process. To increase
      * the calculation speed, I want to store them here, thus I can reuse them
      * if not stored by the user for later evaluation.
      */
     struct pygsl_array_cache *cache;
     /*
      * The Python callback methods.
      */
     PyObject* cbs[PyGSL_SOLVER_NCBS_MAX];
     /*
      * Additional arguments passed to the callbacks
      */
     PyObject* args;
     /*
      * The solver itself.
      */
     void * solver;
     /*
      * The space needed to store the variables for the C functions ....
      */
     void * c_sys;
     /*
      * The dimensionality of the problem. Typically one or two numbers ...
      */
     int problem_dimensions[PyGSL_SOLVER_PB_ND_MAX];
     /*
      * The methods of the solver.
      */
     const struct _SolverStatic* mstatic;
     /* before one can iterate the solver the method set() must be called!*/
     int set_called;
     /* Used as a flag if the jmp_buf is set */
     int isset;
};
typedef struct _PyGSLSolverObject PyGSL_solver;

typedef struct{
     const void * type;
     void* alloc;
     const struct _SolverStatic* mstatic;
} solver_alloc_struct;

#ifndef _PyGSL_SOLVER_API_MODULE
#define PyGSL_SOLVER_API_EXTERN extern
#else
#define PyGSL_SOLVER_API_EXTERN static
#endif

/*
 *  Initalises a solver
 *
 *  The internal structure is set up and the
 *  nd: number of dimensions
 *     0 ... zero dimensional e.g. minimize, root
 *     1 ... one  dimensional e.g. multimin
 *     2 ... two  dimensional e.g. multifit_nlin
 *
 *     3 ... only initalise the structure
 *
 */
PyGSL_SOLVER_API_EXTERN  PyObject *
PyGSL_solver_dn_init(PyObject *self, PyObject *args, const solver_alloc_struct * alloc, int nd);

/*
 * Accessor methods.
 *
 * These methods allow to access parameters of the C structure using the access
 * methods _m_t func
 */
PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_ret_double(PyGSL_solver *self, PyObject *args, double_m_t func);

PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_ret_int(PyGSL_solver *self, PyObject *args, int_m_t func);

PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_ret_size_t(PyGSL_solver *self, PyObject *args, size_t_m_t func);

PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_ret_vec(PyGSL_solver *self, PyObject *args,  ret_vec func);

/*
 * evaluates a C function taking an vector and a double as input and returning a status.
 */

#define PyGSL_CALLABLE_CHECK(ob, name) \
(PyCallable_Check(ob) ? GSL_SUCCESS : PyGSL_Callable_Check(ob, name))

PyGSL_SOLVER_API_EXTERN int
PyGSL_Callable_Check(PyObject *f, const char * myname);

PyGSL_SOLVER_API_EXTERN int
PyGSL_function_wrap_On_O(const gsl_vector * x, PyObject *callback,
			 PyObject *arguments, double *result1,
			 gsl_vector *result2, int n, const char * c_func_name);

PyGSL_SOLVER_API_EXTERN int
PyGSL_function_wrap_OnOn_On(const gsl_vector *x, const gsl_vector *v, gsl_vector *hv, PyObject *callback,
			    PyObject *arguments, int n, const char *c_func_name);

PyGSL_SOLVER_API_EXTERN int
PyGSL_function_wrap_Op_On_Opn(const gsl_vector *x, gsl_vector *f1, gsl_matrix *f2, PyObject *callback,
			    PyObject *arguments, int n, int p, const char *c_func_name);

PyGSL_SOLVER_API_EXTERN int
PyGSL_function_wrap_Op_On(const gsl_vector * x, gsl_vector *f, PyObject *callback,
			  PyObject * arguments, int n, int p, const char *c_func_name);

PyGSL_SOLVER_API_EXTERN int
PyGSL_function_wrap_Op_Opn(const gsl_vector * x, gsl_matrix *f, PyObject *callback,
			   PyObject *arguments, int n, int p, const char * c_func_name);

/*
 * evaluates a C function taking an vector and a double as input and returning a status.
 */
PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_vd_i(PyObject * self, PyObject *args, int_f_vd_t func);

PyGSL_SOLVER_API_EXTERN PyObject *
PyGSL_solver_vvdd_i(PyObject * self, PyObject * args, int_f_vvdd_t func);

struct pygsl_solver_n_set{
     int is_fdf;
     void *c_sys;
     set_m_t set;
};

PyGSL_SOLVER_API_EXTERN PyObject *
PyGSL_solver_n_set(PyGSL_solver *self, PyObject *pyargs, PyObject *kw,
		   const struct pygsl_solver_n_set * info);

/*
 *
 */
PyGSL_SOLVER_API_EXTERN int
PyGSL_solver_func_set(PyGSL_solver *self, PyObject *args, PyObject *f,
		       PyObject *df, PyObject *fdf);

PyGSL_SOLVER_API_EXTERN PyObject*
PyGSL_solver_set_f(PyGSL_solver *self, PyObject *pyargs, PyObject *kw,
		   void *fptr, int isfdf);

#define _GET(name, cast, func) \
PyObject * PyGSL_ ## name(PyGSL_solver *self, PyObject *args) \
{ \
    return PyGSL_solver_ ## func(self, args, (cast) gsl_ ## name); \
}

#define GETDOUBLE(name) _GET(name, double_m_t, ret_double)
#define GETSIZET(name)  _GET(name, size_t_m_t, ret_size_t)
#define GETINT(name)    _GET(name, int_m_t,    ret_int)
#define GETVEC(name)    _GET(name, ret_vec,    ret_vec)

/*
 * Get or set a double ....
 */
enum PyGSL_GETSET_typemode {
     PyGSL_MODE_DOUBLE = 0,
     PyGSL_MODE_INT,
     PyGSL_MODE_SIZE_T
};

PyGSL_SOLVER_API_EXTERN PyObject *
PyGSL_solver_GetSet(PyObject *self, PyObject *args, void * address, enum PyGSL_GETSET_typemode mode);

#ifndef _PyGSL_SOLVER_API_MODULE

#define PyGSL_function_wrap_Op_On \
(* (int (*)(const gsl_vector *, gsl_vector *, PyObject *, PyObject *, int, int, const char *))\
 PyGSL_API[PyGSL_function_wrap_Op_On_NUM])

#define PyGSL_function_wrap_On_O \
(* (int (*)(const gsl_vector *, PyObject *, PyObject *, double *, gsl_vector *, int, const char *))\
 PyGSL_API[PyGSL_function_wrap_On_O_NUM])

#define PyGSL_function_wrap_OnOn_On \
(* (int (*)(const gsl_vector *, const gsl_vector *, gsl_vector *, PyObject *,  PyObject *,int, const char *))\
 PyGSL_API[PyGSL_function_wrap_OnOn_On_NUM])

#define PyGSL_function_wrap_Op_On_Opn \
(* (int (*)(const gsl_vector *, gsl_vector *, gsl_matrix *, PyObject *, PyObject *, int, int, const char *))\
 PyGSL_API[PyGSL_function_wrap_Op_On_Opn_NUM])

#define PyGSL_function_wrap_Op_Opn \
(* (int (*)(const gsl_vector *, gsl_matrix *, PyObject *, PyObject *, int, int, const char *))\
 PyGSL_API[PyGSL_function_wrap_Op_Opn_NUM])

#define PyGSL_solver_ret_int \
(*(PyObject *  (*) (PyGSL_solver *, PyObject *, int_m_t func))    PyGSL_API[PyGSL_solver_ret_int_NUM])

#define PyGSL_solver_ret_double \
(*(PyObject *  (*) (PyGSL_solver *, PyObject *, double_m_t func)) PyGSL_API[PyGSL_solver_ret_double_NUM])

#define PyGSL_solver_ret_size_t \
(*(PyObject *  (*) (PyGSL_solver *, PyObject *, size_t_m_t func)) PyGSL_API[PyGSL_solver_ret_size_t_NUM])

#define PyGSL_solver_ret_vec \
(*(PyObject *  (*) (PyGSL_solver *, PyObject *, ret_vec func))    PyGSL_API[PyGSL_solver_ret_vec_NUM])

#define  PyGSL_solver_dn_init \
(*(PyObject * (*)(PyObject *, PyObject *, const solver_alloc_struct *, int))PyGSL_API[PyGSL_solver_dn_init_NUM])

#define PyGSL_Callable_Check \
(*(int (*)(PyObject *, const char *)) PyGSL_API[PyGSL_Callable_Check_NUM])

#define PyGSL_solver_vd_i \
(*(PyObject * (*) (PyObject *, PyObject *, int_f_vd_t)) PyGSL_API[PyGSL_solver_vd_i_NUM])

#define PyGSL_solver_vvdd_i \
(*(PyObject * (*) (PyObject *, PyObject *, int_f_vvdd_t)) PyGSL_API[PyGSL_solver_vvdd_i_NUM])

#define PyGSL_solver_func_set \
(*(int (*)(PyGSL_solver *, PyObject *, PyObject *, PyObject *, PyObject *)) PyGSL_API[PyGSL_solver_func_set_NUM])

#define PyGSL_solver_n_set \
(*(PyObject * (*)(PyGSL_solver *, PyObject *, PyObject *,  const struct pygsl_solver_n_set *)) PyGSL_API[PyGSL_solver_n_set_NUM])

#define PyGSL_solver_set_f \
(* (PyObject * (*)(PyGSL_solver *, PyObject *, PyObject *, void *, int )) PyGSL_API[PyGSL_solver_set_f_NUM])

#define PyGSL_solver_GetSet \
(* (PyObject * (*) (PyObject *, PyObject *, void *, enum PyGSL_GETSET_typemode mode)) PyGSL_API[PyGSL_solver_getset_NUM])

#define PyGSL_solver_check(ob) (((ob)->ob_type) == (PyGSL_API[PyGSL_solver_type_NUM]))

#define import_pygsl_solver() \
{ \
   init_pygsl(); \
   if (PyImport_ImportModule("pygsl.testing.solver") != NULL) { \
          ;\
   } else { \
        fprintf(stderr, "failed to import pygsl solver!!\n"); \
   } \
}

#else  /* _PyGSL_API_MODULE */
#define PyGSL_solver_check(ob) ((ob)->ob_type == &PyGSL_solver_pytype)
#endif /* _PyGSL_API_MODULE */

#endif /* _PYGSL_SOLVER_H_ */

ext: h

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/**
 *
 * @file core_cpltmg.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Ichitaro Yamazaki
 * @author Julien Herrmann
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:47 2014
 *
 **/

#include <math.h>
#include <lapacke.h>
#include "common.h"

#define pi (3.1415926535897932384626433832795028841971693992)

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_cpltmg initialize a tile of a random matrix from the MatLab
 *  gallery configured with the default parameters, and a few other
 *  specific matrices.
 *
 *******************************************************************************
 *
 * @param[in] mtxtype
 *          Possible types are: PlasmaMatrixRandom, PlasmaMatrixHadamard,
 *          PlasmaMatrixParter, PlasmaMatrixRis, PlasmaMatrixKms,
 *          PlasmaMatrixMoler, PlasmaMatrixCompan, PlasmaMatrixRiemann,
 *          PlasmaMatrixLehmer, PlasmaMatrixMinij, PlasmaMatrixDorr,
 *          PlasmaMatrixDemmel, PlasmaMatrixInvhess, PlasmaMatrixCauchy,
 *          PlasmaMatrixHilb, PlasmaMatrixLotkin, PlasmaMatrixOrthog,
 *          PlasmaMatrixWilkinson, PlasmaMatrixFoster, PlasmaMatrixWright,
 *          PlasmaMatrixLangou
 *          (See further in the code for more details)
 *
 * @param[in] M
 *         The number of rows of the tile A. M >= 0.
 *
 * @param[in] N
 *         The number of columns of the tile A. N >= 0.
 *
 * @param[in,out] A
 *         On entry, the M-by-N tile to be initialized.
 *         On exit, the tile initialized in the mtxtype format.
 *
 * @param[in] LDA
 *         The leading dimension of the tile A. LDA >= max(1,M).
 *
 * @param[in] gM
 *         The global number of rows of the full matrix, A is belonging to. gM >= (m0+M).
 *
 * @param[in] gN
 *         The global number of columns of the full matrix, A is belonging to. gN >= (n0+gN).
 *
 * @param[in] m0
 *         The index of the first row of tile A in the full matrix. m0 >= 0.
 *
 * @param[in] n0
 *         The index of the first column of tile A in the full matrix. n0 >= 0.
 *
 * @param[in] seed
 *         The seed used for random generation. Must be the same for
 *         all tiles initialized with this routine.
 *
 *******************************************************************************
 *
 * @return
 *         \retval PLASMA_SUCCESS successful exit
 *         \retval <0 if INFO = -k, the k-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cpltmg = PCORE_cpltmg
#define CORE_cpltmg PCORE_cpltmg
#define CORE_cplrnt  PCORE_cplrnt
void
CORE_cplrnt( int M, int N, PLASMA_Complex32_t *A, int LDA,
             int gM, int m0, int n0,
             unsigned long long int seed );
#endif
int CORE_cpltmg( PLASMA_enum mtxtype,
                  int M, int N, PLASMA_Complex32_t *A, int LDA,
                  int gM, int gN, int m0, int n0,
                  unsigned long long int seed )
{
    int i, j;

    /* Check input arguments */
    if (M < 0) {
        coreblas_error(2, "Illegal value of M");
        return -3;
    }
    if (N < 0) {
        coreblas_error(3, "Illegal value of N");
        return -3;
    }
    if ((LDA < max(1,M)) && (M > 0)) {
        coreblas_error(5, "Illegal value of LDA");
        return -5;
    }
    if (m0 < 0) {
        coreblas_error(8, "Illegal value of m0");
        return -8;
    }
    if (n0 < 0) {
        coreblas_error(9, "Illegal value of n0");
        return -9;
    }
    if (gM < m0+M) {
        coreblas_error(6, "Illegal value of gM");
        return -6;
    }
    if (gN < n0+N) {
        coreblas_error(7, "Illegal value of gN");
        return -7;
    }

    /* Quick return */
    if ((M == 0) || (N == 0))
        return PLASMA_SUCCESS;

    switch( mtxtype ) {
    case PlasmaMatrixRandom:
    {
        CORE_cplrnt( M, N, A, LDA, gM, m0, n0, seed );
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/hadamard.html
     *
     * Initialize the tile A to create the Hadamard matrix of order gN.
     *
     * Hadamard matrices are matrices of 1's and -1's whose columns are orthogonal,
     *
     *   H'*H = gN*I
     *
     *   where [gN gN]=size(H) and I = eye(gN,gN) ,.
     *
     *  They have applications in several different areas, including
     *  combinatorics, signal processing, and numerical analysis.
     *
     *  An n-by-n Hadamard matrix with n > 2 exists only if rem(n,4) =
     *  0. This function handles only the cases where n is a power of
     *  2.
     */
    case PlasmaMatrixHadamard:
    {
        int tmp, nbone;

        /* Extra parameters check */
        if (gM != gN) {
            coreblas_error(6, "Illegal value of gM (Matrix must be square)");
            return -6;
        }

        tmp = gM;
        while ( tmp > 1 ) {
            if( tmp % 2 != 0 ) {
                coreblas_error(6, "Illegal value of gM (Matrix dimension must be a power of 2)");
                return -6;
            }
            tmp /= 2;
        }

        for (j=0; j<N; j++) {
            for (i=0; i<M; i++) {
                tmp = ((m0 + i) & (n0 + j));
                nbone = 0;
                while ( tmp != 0 )
                {
                    nbone += ( tmp & 1 );
                    tmp >>= 1;
                }
                A[j*LDA+i] = (PLASMA_Complex32_t)(1. - 2. * ( nbone % 2 ));
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000116
     *
     * Toeplitz matrix with singular values near pi.
     * Returns the tile A, such that the elment of the matrix are 1/(i-j+0.5).
     *
     * C is a Cauchy matrix and a Toeplitz matrix. Most of the
     * singular values of C are very close to pi.
     *
     */
    case PlasmaMatrixParter:
    {
        PLASMA_Complex32_t tmp;

        if (gM != gN) {
            coreblas_error(6, "Illegal value of gM (Matrix must be square for Parter)");
            return -6;
        }

        tmp = (PLASMA_Complex32_t)( .5 + m0 - n0 );
        for (j=0; j<N; j++) {
            for (i=0; i<M; i++) {
                A[j*LDA+i] = (PLASMA_Complex32_t)1. / (PLASMA_Complex32_t)( tmp + i - j );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000243
     *
     * Symmetric Hankel matrix
     * Returns a symmetric gN-by-gN Hankel matrix with elements.
     *
     *     A(i,j) = 0.5/(n-i-j+1.5)
     *
     * The eigenvalues of A cluster around π/2 and –π/2. This matrix
     * was invented by F.N. Ris.
     *
     */
     case PlasmaMatrixRis:
    {
        PLASMA_Complex32_t tmp;

        if (gM != gN) {
            coreblas_error(6, "Illegal value of gM (Matrix must be square for RIS)");
            return -6;
        }

        tmp = (PLASMA_Complex32_t)( gM - m0 - n0 - 0.5 );
        for (j=0; j<N; j++) {
            for (i=0; i<M; i++) {
                A[j*LDA+i] = (PLASMA_Complex32_t).5 / (PLASMA_Complex32_t)( tmp - i - j );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000026
     *
     * Kac-Murdock-Szego Toeplitz matrix
     *
     * Returns the n-by-n Kac-Murdock-Szego Toeplitz matrix such that
     * A(i,j) = rho^(abs(i-j)), for real rho.
     *
     * For complex rho, the same formula holds except that elements
     * below the diagonal are conjfugated. rho defaults to 0.5.
     *
     * The KMS matrix A has these properties:
     *     - An LDL' factorization with L inv(gallery('triw',n,-rho,1))',
     *       and D(i,i) (1-abs(rho)^2)*eye(n), except D(1,1) = 1.
     *     - Positive definite if and only if 0 < abs(rho) < 1.
     *     - The inverse inv(A) is tridiagonal.Symmetric Hankel matrix
     *
     * In this function, rho is set to 0.5 and cannot be changed.
     */
    case PlasmaMatrixKms:
    {
        PLASMA_Complex32_t rho;

        if (gM != gN) {
            coreblas_error(6, "Illegal value of gM (Matrix must be square for KMS)");
            return -6;
        }

        rho = .5;
        for (j=0; j<N; j++) {
            for (i=0; i<M; i++) {
                A[j*LDA+i] = (PLASMA_Complex32_t)( cpowf( rho, fabs( (float)( m0 + i - n0 - j ) ) ) );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000074
     *
     * Symmetric positive definite matrix
     *
     * Returns the symmetric positive definite n-by-n matrix U'*U,
     * where U = gallery('triw',n,alpha).
     *
     * For the default alpha = -1, A(i,j) = min(i,j)-2, and A(i,i) =
     * i. One of the eigenvalues of A is small.
     *
     */
    case PlasmaMatrixMoler:
    {
        int ii, jj;
        for (j=0,jj=n0; j<N; j++,jj++) {
            for (i=0,ii=m0; i<M; i++,ii++) {
                if ( ii == jj ) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( ii + 1. );
                } else {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( min( ii, jj ) - 1. );
                }
            }
        }
    }
    break;

    /*
     *  See http://www.mathworks.fr/fr/help/matlab/ref/compan.html
     *
     *  Companion matrix
     *
     *  A = compan(u) returns the corresponding companion matrix whose first row is
     *  -u(2:n)/u(1), where u is a vector of polynomial coefficients. The
     *  eigenvalues of compan(u) are the roots of the polynomial.
     *
     */
    case PlasmaMatrixCompan:
    {
        int jj;

        LAPACKE_claset_work(LAPACK_COL_MAJOR, 'A', M, N, 0., 0., A, LDA);

        /* First row */
        if ( m0 == 0 ) {
            PLASMA_Complex32_t v0;
            /* Get V0 */
            CORE_cplrnt( 1, 1, &v0, 1, 1, 1, 0, seed );
            v0 = 1. / v0 ;

            /* Initialize random vector */
            CORE_cplrnt( 1, N, A, LDA, 1, 1, n0, seed );
            cblas_cscal( N, CBLAS_SADDR(v0), A, LDA);

            /* Restore A(0,0) */
            if ( n0 == 0 )
                A[0] = 0.;
        }

        /* Sub diagonal*/
        for (j=0,jj=n0; j<N; j++,jj++)
        {
            i = jj + 1 - m0;
            if ( ( i > 0 ) && (i < M) )
            {
                A[j*LDA+i] = 1.;
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000232
     *
     * Matrix associated with the Riemann hypothesis
     *
     * Returns an n-by-n matrix for which the Riemann hypothesis is
     * true if and only if for every eps > 0.
     *
     * The Riemann matrix is defined by:
     *
     *    A = B(2:n+1,2:n+1)
     *
     *    where B(i,j) = i-1 if i divides j, and B(i,j) = -1 otherwise.
     *
     * The Riemann matrix has these properties:
     *     - Each eigenvalue e(i) satisfies abs(e(i)) <= m-1/m, where m = n+1.
     *     - i <= e(i) <= i+1 with at most m-sqrt(m) exceptions.
     *     - All integers in the interval (m/3, m/2] are eigenvalues.
     *
     */
    case PlasmaMatrixRiemann:
    {
        int ii, jj;
        for (j=0,jj=n0+2; j<N; j++,jj++) {
            for (i=0,ii=m0+2; i<M; i++,ii++) {
                if ( jj%ii == 0 ) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( ii - 1. );
                } else {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( -1. );
                }
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000049
     *
     * Symmetric positive definite matrix
     *
     * Returns the symmetric positive definite n-by-n matrix such that
     * A(i,j) = i/j for j >= i.
     *
     * The Lehmer matrix A has these properties:
     *     - A is totally nonnegative.
     *     - The inverse inv(A) is tridiagonal and explicitly known.
     *     - The order n <= cond(A) <= 4*n*n.Matrix associated with the
     *       Riemann hypothesis
     *
     */
     case PlasmaMatrixLehmer:
    {
        int ii, jj;
        for (j=0,jj=n0+1; j<N; j++,jj++) {
            for (i=0,ii=m0+1; i<M; i++,ii++) {
                if ( jj >= ii ) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( ii ) / (PLASMA_Complex32_t)( jj );
                } else {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( jj ) / (PLASMA_Complex32_t)( ii );
                }
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000066
     *
     * Symmetric positive definite matrix
     *
     * Returns the n-by-n symmetric positive definite matrix with
     * A(i,j) = min(i,j).
     *
     * The minij matrix has these properties:
     *     - The inverse inv(A) is tridiagonal and equal to -1 times the
     *       second difference matrix, except its (n,n) element is 1.
     *     - Givens' matrix, 2*A-ones(size(A)), has tridiagonal inverse
     *       and eigenvalues 0.5*sec((2*r-1)*pi/(4*n))^2, where r=1:n.
     *     - (n+1)*ones(size(A))-A has elements that are max(i,j) and a
     *       tridiagonal inverse.
     *
     */
    case PlasmaMatrixMinij:
    {
        int ii, jj;
        for (j=0,jj=n0+1; j<N; j++,jj++) {
            for (i=0,ii=m0+1; i<M; i++,ii++) {
                A[j*LDA+i] = (PLASMA_Complex32_t) min( ii, jj );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-999936
     *
     * Diagonally dominant, ill-conditioned, tridiagonal matrix
     *
     * Returns the n-by-n matrix, row diagonally dominant, tridiagonal
     * matrix that is ill-conditioned for small nonnegative values of
     * theta. The default value of theta is 0.01. The Dorr matrix
     * itself is the same as gallery('tridiag',c,d,e).
     *
     */
    case PlasmaMatrixDorr:
    {
        PLASMA_Complex32_t theta = 0.01;
        PLASMA_Complex32_t h     = 1. / ( gN + 1. );
        PLASMA_Complex32_t term  = theta / ( h * h );
        int jj;
        int half = (gN+1) / 2;

        LAPACKE_claset_work(LAPACK_COL_MAJOR, 'A', M, N, 0., 0., A, LDA);

        /* First half */
        for (j=0, jj=n0; (j<N) && (jj<half); j++, jj++) {
            i = jj - m0;
            if ( ( i < -1 ) || (i > M) )
                continue;

            /* Over the diagonal */
            if (i > 0)
                A[j*LDA + i-1] = - term - (0.5 - jj*h)/h;

            /* Diagonal */
            if ( i >= M )
                return PLASMA_SUCCESS;

            if ( i >=0 )
                A[j*LDA + i] = 2. * term + (0.5 - (jj+1) * h) / h;

            /* Below the diagonal */
            if (i+1 < M) {
                if (jj+1 == half)
                    A[j*LDA + i+1] = - term + (0.5 - (jj+2)*h)/h;
                else
                    A[j*LDA + i+1] = - term;
            }
        }

        /* Second half */
        for (; j<N; j++,jj++) {
            i = jj - m0;
            if ( ( i < -1 ) || (i > M))
                continue;

            if (i > 0) {
                if (jj == half)
                    A[j*LDA + i-1] = - term - (0.5 - jj*h)/h;
                else
                    A[j*LDA + i-1] = - term;
            }

            if ((i >=0) && (i < M))
                A[j*LDA + i] = 2. * term - (0.5 - (jj+1)*h)/h;

            if (i+1 < M)
                A[j*LDA + i+1] = - term + (0.5 - (jj+2)*h)/h;
        }
    }
    break;

    /*
     * See [1] J. Demmel, Applied Numerical Linear Algebra, SIAM,
     *         Philadelphia, 1997
     *
     * Returns a matrix defined by:
     *    A = D * ( I + 1e-7* rand(n)), where D = diag(10^{14*(0:n-1)/n})
     *
     */
    case PlasmaMatrixDemmel:
    {
        PLASMA_Complex32_t dii;
        int ii, jj;

        /* Randomize the matrix */
        CORE_cplrnt( M, N, A, LDA, gM, m0, n0, seed );

        for (j=0,jj=n0; j<N; j++,jj++) {
            for (i=0,ii=m0; i<M; i++,ii++) {
                dii = cpowf( 10. , 14. * ii / gM );
                A[j*LDA+i] *= dii * ( (jj == ii) ? 1. : 1.e-7 );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000000
     *
     * Inverse of an upper Hessenberg matrix
     *
     * A = gallery('invhess',x,y), where x is a length n vector and y
     * is a length n-1 vector, returns the matrix whose lower triangle
     * agrees with that of ones(n,1)*x' and whose strict upper
     * triangle agrees with that of [1 y]*ones(1,n).
     *
     * The matrix is nonsingular if x(1) ~= 0 and x(i+1) ~= y(i) for
     * all i, and its inverse is an upper Hessenberg matrix. Argument
     * y defaults to -x(1:n-1).
     *
     * If x is a scalar, invhess(x) is the same as invhess(1:x).
     *
     * Here: gallery('invhess', gM)
     */
    case PlasmaMatrixInvhess:
    {
        int ii, jj;
        for (j=0,jj=n0+1; j<N; j++,jj++) {
            for (i=0,ii=m0+1; i<M; i++,ii++) {
                if ( jj <= ii ) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( jj );
                } else {
                    A[j*LDA+i] = (PLASMA_Complex32_t)( -ii );
                }
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1019317
     *
     * Cauchy matrix
     *
     * Returns an n-by-n matrix C such that, C(i,j) = 1/(i + j).
     *
     */
    case PlasmaMatrixCauchy:
    {
        int ii, jj;
        for (j=0,jj=n0+1; j<N; j++,jj++) {
            for (i=0,ii=m0+1; i<M; i++,ii++) {
                A[j*LDA+i] = (PLASMA_Complex32_t)( 1. ) / (PLASMA_Complex32_t)( ii+jj );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/hilb.html
     *
     * Hilbert Matrix
     *
     * The Hilbert matrix is a notable example of a poorly conditioned
     * matrix. The elements of the Hilbert matrices are:
     *   H(i,j) = 1/(i * + j – 1).
     *
     */
    case PlasmaMatrixHilb:
    {
        PLASMA_Complex32_t tmp = (PLASMA_Complex32_t)( m0 + n0 + 1. );
        for (j=0; j<N; j++) {
            for (i=0; i<M; i++) {
                A[j*LDA+i] = (PLASMA_Complex32_t)1. / (PLASMA_Complex32_t)( tmp + i + j );
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000062
     *
     * Lotkin matrix
     *
     * Returns the Hilbert matrix with its first row altered to all
     * ones. The Lotkin matrix A is nonsymmetric, ill-conditioned, and
     * has many negative eigenvalues of small magnitude. Its inverse
     * has integer entries and is known explicitly.
     *
     */
    case PlasmaMatrixLotkin:
    {
        PLASMA_Complex32_t tmp = (PLASMA_Complex32_t)( m0 + n0 + 1. );
        if (m0 == 0) {
            for (j=0; j<N; j++) {
                A[j*LDA] = (PLASMA_Complex32_t)1.;
                for (i=1; i<M; i++) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)1. / (PLASMA_Complex32_t)( tmp + i + j );
                }
            }
        } else {
            for (j=0; j<N; j++) {
                for (i=0; i<M; i++) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)1. / (PLASMA_Complex32_t)( tmp + i + j );
                }
            }
        }
    }
    break;

    /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/gallery.html#f84-1000083
     *
     * Orthogonal and nearly orthogonal matrices
     *
     * Returns the matrix Q of order n, such that:
     *    Q(i,j) = sqrt(2/(n+1)) * sin(i*j*pi/(n+1))
     *
     * Symmetric eigenvector matrix for second difference matrix.
     *
     */
    case PlasmaMatrixOrthog: /* Default: k=1 */
    {
        PLASMA_Complex32_t sqrtn = (PLASMA_Complex32_t) sqrt( 2. / (gN+1.) );
        float scale = pi / (float)(gN+1.);

        int ii, jj;
        for (j=0,jj=n0+1; j<N; j++,jj++) {
            for (i=0,ii=m0+1; i<M; i++,ii++) {
                A[j*LDA+i] = sqrtn * (PLASMA_Complex32_t) sin( (float)ii * (float)jj * scale );
            }
        }
    }
    break;

     /*
     * See http://www.mathworks.fr/fr/help/matlab/ref/wilkinson.html
     *
     * Wilkinson's eigenvalue test matrix
     *
     * Returns one of J. H. Wilkinson's eigenvalue test matrices. It
     * is a symmetric, tridiagonal matrix with pairs of nearly, but
     * not exactly, equal eigenvalues.
     *
     */
   case PlasmaMatrixWilkinson:
    {
        if (gM != gN) {
            coreblas_error(6, "Illegal value of gM (Matrix must be square for Wilkinson)");
            return -6;
        }

        int ii, jj;
        for (j=0,jj=n0; j<N; j++,jj++) {
            for (i=0,ii=m0; i<M; i++,ii++) {
                if (ii == jj) {
                    PLASMA_Complex32_t tmp = (PLASMA_Complex32_t)(( (gN - 1 - ii) < ii ) ? gN - 1 - ii : ii );
                     A[j*LDA+i] = (PLASMA_Complex32_t)(gN - 2. * tmp - 1.) / 2.;
                }
                else if ( (ii == jj+1) || (ii == jj-1) ) {
                     A[j*LDA+i] = (PLASMA_Complex32_t)1.;
                }
                else {
                     A[j*LDA+i] = (PLASMA_Complex32_t)0.;
                }
           }
       }
    }
    break;

    /*
     * See [1] L. V. Foster, Gaussian Elimination with Partial
     *         Pivoting Can Fail in Practice, SIAM J. Matrix
     *         Anal. Appl., 15 (1994), pp. 1354-1362.
     *
     *     [2] L. V. Foster, The growth factor and efficiency of
     *         Gaussian elimination with rook pivoting,
     *         J. Comput. Appl. Math., 86 (1997), pp. 177-194
     *
     * A pathological case for LU with gaussian elimination.
     *
     */
    case PlasmaMatrixFoster: /* Default: k=h=c=1 */
    {
        float k=1., h=1., c=1.;

        int ii, jj;
        for (j=0,jj=n0; j<N; j++,jj++) {
            for (i=0,ii=m0; i<M; i++,ii++) {

                if (ii == jj) {
                    if (jj == 0)
                        A[j*LDA+i] = (PLASMA_Complex32_t)1.;
                    else if (jj == gN-1)
                        A[j*LDA+i] = (PLASMA_Complex32_t)(1. - (1. / c) - (k*h)/2. );
                    else
                        A[j*LDA+i] = (PLASMA_Complex32_t)(1. - (k*h)/2. );
                }
                else if (jj == 0) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-k*h/2.);
                }
                else if (jj == gN-1) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-1./c);
                }
                else if (ii > jj) {
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-k*h);
                }
                else {
                    A[j*LDA+i] = (PLASMA_Complex32_t)0.;
                }
            }
        }
    }
    break;

    /*
     * See [3] S. J. Wright, A collection of problems for which
     *         Gaussian elimination with partial pivoting is unstable,
     *         SIAM J. SCI. STATIST. COMPUT., 14 (1993), pp. 231-238.
     *
     * A pathological case for LU with gaussian elimination.
     *
     */
    /*
     * Default: h=0.01, M=[-10 -19, 19 30]. Then,
     *   exp(h*M)=[0.9048 0.8270, 1.2092 1.3499]
     */
    case PlasmaMatrixWright:
    {
        int ii, jj;
        for (j=0,jj=n0; j<N; j++,jj++) {
            for (i=0,ii=m0; i<M; i++,ii++) {

                if (ii == jj)
                    A[j*LDA+i] = (PLASMA_Complex32_t)1.;
                else if ((ii == jj + 2) && (jj % 2 == 0))
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-0.9048);
                else if ((ii == jj + 3) && (jj % 2 == 0))
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-1.2092);
                else if ((ii == jj + 2) && (jj % 2 == 1))
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-0.8270);
                else if ((ii == jj + 3) && (jj % 2 == 1))
                    A[j*LDA+i] = (PLASMA_Complex32_t)(-1.3499);
                else if ((jj == gM-2) && (ii == 0))
                    A[j*LDA+i] = (PLASMA_Complex32_t)1.;
                else if ((jj == gM-1) && (ii == 1))
                    A[j*LDA+i] = (PLASMA_Complex32_t)1.;
                else
                    A[j*LDA+i] = (PLASMA_Complex32_t)0.;
            }
        }
    }
    break;

    /*
     * Generates a pathological case for LU with gaussian elimination.
     *
     * Returns a random matrix on which, the columns from N/4 to N/2
     * are scaled down by eps.
     * These matrices fails on LU with partial pivoting, but Hybrid
     * LU-QR algorithms manage to recover the scaled down columns.
     *
     */
    case PlasmaMatrixLangou:
    {
        PLASMA_Complex32_t eps = (PLASMA_Complex32_t) LAPACKE_slamch_work( 'e' );
        int ii, jj, mm, minMN;
        int firstcol, lastcol;

         /* Create random tile */
         CORE_cplrnt( M, N, A, LDA, gM, m0, n0, seed );

         /* Scale down the columns gN/4 to gN/2 below the diagonal */
         minMN = min( gM, gN );
         firstcol = minMN / 4;
         lastcol  = minMN / 2;

         if ( (m0 >= n0) && ((n0+N) >= firstcol ) && (n0 < lastcol) ) {

             jj = max( n0, firstcol );
             j = jj - n0;

             for (; j<N && jj < lastcol; j++,jj++) {

                 ii = max( m0, jj );
                 i  = ii - m0;
                 mm = M - i;
                 cblas_cscal( mm, CBLAS_SADDR(eps), A + j*LDA + i, 1);
             }
         }
    }
    break;

    default:
        coreblas_error(1, "Illegal value of mtxtype");
        return -1;
    }

    return PLASMA_SUCCESS;
}

ext: c


text:
#ifndef libraries_h
#define libraries_h

#include <fstream>
#include <iostream>
#include <iomanip>
#include <ctime>
#include <math.h>
#include <vector>
#include <deque>
#include <queue>
#include <stack>
#include <array>
#include <algorithm>
#include <unordered_map>
#include <map>
#include <memory>
#include <utility>
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_statistics.h>
#include "types.h"
//#include <omp.h>
//#include <tbb/tbb.h>

#endif


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

#pragma once

#include <gsl/span>

namespace onnxruntime {
namespace contrib {

template <typename T>
class IAttentionMechanism {
 public:
  virtual ~IAttentionMechanism() = default;

  virtual void PrepareMemory(
    const gsl::span<const T>& memory,
    const gsl::span<const int>& memory_sequence_lengths) = 0;

  virtual void Compute(
    const gsl::span<const T>& query,
    const gsl::span<const T>& prev_alignment,
    const gsl::span<T>& output,
    const gsl::span<T>& alignment) const = 0;

  virtual const gsl::span<const T> Values() const = 0;

  virtual const gsl::span<const T> Keys() const = 0;

  virtual int GetMaxMemorySteps() const = 0;

  virtual bool NeedPrevAlignment() const = 0;
};

}
}

max_stars_repo_name: hqucms/onnxruntime

max_issues_repo_name: hqucms/onnxruntime

max_forks_repo_name: hqucms/onnxruntime


text:
#pragma once

#include <gsl\gsl>
#include <winrt\base.h>
#include <d3d11.h>
#include <DirectXTK\SpriteBatch.h>
#include "DrawableGameComponent.h"
#include "HeightmapTessellationMaterial.h"
#include "RenderStateHelper.h"

namespace Rendering
{
	class HeightmapTessellationDemo final : public Library::DrawableGameComponent
	{
	public:
		HeightmapTessellationDemo(Library::Game& game, const std::shared_ptr<Library::Camera>& camera);
		HeightmapTessellationDemo(const HeightmapTessellationDemo&) = delete;
		HeightmapTessellationDemo(HeightmapTessellationDemo&&) = default;
		HeightmapTessellationDemo& operator=(const HeightmapTessellationDemo&) = default;
		HeightmapTessellationDemo& operator=(HeightmapTessellationDemo&&) = default;
		~HeightmapTessellationDemo();

		bool AnimationEnabled() const;
		void SetAnimationEnabled(bool enabled);

		gsl::span<const float> EdgeFactors() const;
		gsl::span<const float> InsideFactors() const;
		void SetUniformFactors(float factor);

		float DisplacementScale() const;
		void SetDisplacementScale(float displacementScale);

		virtual void Initialize() override;
		virtual void Update(const Library::GameTime& gameTime) override;
		virtual void Draw(const Library::GameTime& gameTime) override;

	private:
		inline static const RECT HeightmapDestinationRectangle = { 0, 512, 256, 768 };

		Library::RenderStateHelper mRenderStateHelper;
		HeightmapTessellationMaterial mMaterial;
		winrt::com_ptr<ID3D11Buffer> mVertexBuffer;
		bool mUpdateMaterial{ true };
		std::unique_ptr<DirectX::SpriteBatch> mSpriteBatch;
		std::shared_ptr<Library::Texture2D> mHeightmap;
		bool mAnimationEnabled{ false };
	};
}


text:
#ifndef __COMPLEX_DEF_H__
#define __COMPLEX_DEF_H__

#include <complex>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_zeta.h>
#include <gsl/gsl_sf_dilog.h>
#include <gsl/gsl_sf_result.h>

std::complex<long double> operator*(const std::complex<long double>& a, const double& b);

std::complex<long double> operator*(const double& a,const std::complex<long double>& b);

std::complex<long double> operator/(const std::complex<long double>& a, const double& b);

std::complex<long double> operator/(const double& a,const std::complex<long double>& b);

std::complex<long double> operator-(const std::complex<long double>& a, const double& b);

std::complex<long double> operator-(const double& a,const std::complex<long double>& b);

std::complex<long double> operator+(const std::complex<long double>& a, const double& b);

std::complex<long double> operator+(const double& a,const std::complex<long double>& b);

//Useful Math functions
long double log_r(long double x);
long double dilog_r(long double x);

const std::complex<long double> II(0.0,1.0);
std::complex<long double> log_c(std::complex<long double> z);
std::complex<long double> dilog_c(std::complex<long double> z);
#endif /* __COMPLEX_DEF_H__ */

max_stars_repo_name: gvita/RHEHpt

max_issues_repo_name: gvita/RHEHpt

max_forks_repo_name: gvita/RHEHpt


text:
/*
 * edges.h -- élek effektív tárolása: lista az id-k szerint rendezve,
 * 	ebben gyorsan lehet keresni + binary heap az időpontok szerint, így
 * 	a legrégebbit gyorsan ki lehet választani
 *
 * Copyright 2013 Kondor Dániel <dani@thinkpad-r9r>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 * MA 02110-1301, USA.
 *
 *
 */

#ifndef EDGES_H
#define EDGES_H

#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif

#include "idlist.h"
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <ctype.h>
#include <sys/mman.h>
//~ #include <gsl/gsl_rng.h>
#include <limits.h>
//~ #include "mt19937.h"

//változás: hashmap használata helyett előre beolvassuk az összes lehetséges élt, sorbarendezve (feltétel:
//	a bemeneti fájlban sorbarendezve kell legyenek, minden él csak egyszer -- ezt egyszerűbb az SQL Serverrel
//	megcsináltatni), ebben tudunk már gyorsan keresni

typedef struct edgehelper_t {
	const idlist* il;
	uint64_t* index;
	uint64_t* outdeg;
} edgehelper;

//heap struktúra az élek tárolásához (összes él, idő szerint rendezve)
typedef struct edge_t { //egy él, heap ezekből, timestamp szerint rendezve
	uint32_t p1; //első pont -> itt az "igazi" ID-ket tároljuk (a hálózatbeli azonosítókat, ezeknek nem kell szekvenciálisnak lenni)
	uint32_t p2; //második pont (p1->p2 él)
	unsigned int offset; //heap-beli offset (heap[i].offset == i mindig) -- vagy: tranzakció ID-k
	unsigned int timestamp; //legutolsó aktivitás
} edge; //méret -- 16 bájt

typedef struct edges_t {
	edge* e; //élek tárolása itt, id-k szerint rendezve
	uint64_t nedges;
	uint64_t edges_size;
	uint64_t edges_grow;
	edgehelper* h;
} edges;

/*
 * egy tömb (unsigned int-ekből), amiben az edge-ekre mutató pointerek vannak, ezek binary heap-ben:
 * 	unsigned int* heap;
 * 	edges[heap[i]].timestamp < edges[heap[2*i+1 / 2*i+2]].timestamp
 * és edges[heap[i]].offset == i
 *
 */

//hash fv. (két 32-bites id-ből egy 64-bites szám, ezeket talán gyorsabb összehasonlítani
static inline uint64_t edgehash(const edge* e, const unsigned int i) {
		uint64_t r = e[i].p1;
		r = r << 32;
		uint64_t r2 = e[i].p2;
		return r | r2;
}

//két él összehasonlítása, eredmény: -1, ha eh2-nek e[i] előtt kell lenni, 1 ha utána, 0, ha megegyeznek
static inline int cmpedge(const edge* e, uint64_t i, uint64_t eh2) {
	uint64_t eh1 = edgehash(e,i);
	if(eh1 < eh2) return 1;
	if(eh1 > eh2) return -1;
	return 0;
}

static inline int cmpedge2(const edge* e, uint64_t i, uint64_t j) {
	uint64_t eh2 = edgehash(e,j);
	return cmpedge(e,i,eh2);
}

//éleket tároló struktúra lefoglalása / növesztése:
//	ha e == 0, akkor új struktúra lefoglalása
//	ha e != 0, akkor növesztés (mremap())
//eredmény: 0, ha hiba történt (nem sikerült memóriát lefoglalni)
//ha a bemeneti e nem volt 0, akkor azt nem szabadítottuk fel
//edges* pointer, ha rendben volt (ha e != 0, akkor e-t adjuk vissza)
edges* edges_grow(edges* e);

//élek felszabadítása
void edges_free(edges* e);

//élek beolvasása egy fájlból -- a fájlban már megfelelően rendezve kell lenniük
edges* edges_read(FILE* f);

//élek beolvasása egy fájlból -- a fájlban már megfelelően rendezve kell lenniük
edges* edges_read2(FILE* f); //+ timestamp-ok is meg vannak adva

//beolvasás általánosan:
//flags értékei:
#define EFLAGS_T1 1 //időpontokat is beolvasunk és idő szerint rendezzük az eredményt (ha nincs megadva, akkor csak
	//éleket olvasunk be és minden élből csak egyet tartunk meg)
#define EFLAGS_SELF 2 //megengedünk önmagába mutató éleket (p1->p1)
#define EFLAGS_TXID 4 //tranzakció ID-ket is beolvasunk (legutolsó oszlop, az offset tömbbe tároljuk el)
#define EFLAGS_ERROR_OVERFLOW 8 // overflow hibát jelent (különben csak kihagyjuk)
edges* edges_read0(FILE* f, int flags);

//élek átmásolása egy új tömbbe, ha r == 1, akkor fordítva (p2->p1, p1->p2)
edges* edges_copy(edges* e1, int r);

//részhalmaz átmásolása
edges* edges_copy2(edges* e1, uint64_t start, uint64_t end, int r);

//él megkeresése, eredmény: a tömbön (e->e) belüli sorszám, ha megtaláltuk, e->nedges, ha nem találtuk
uint64_t edges_find(edges* e, uint32_t p1, uint32_t p2);

//olyan él keresése, ahol az első ID megegyezik a megadottal
uint64_t edges_findfirst(edges* e, uint32_t p1);

//élek sorbarendezése az ID-k szerint
void edges_sort1(edges* e);

//véletlen cserék
//a->b, c->d helyett a->d, c->b
//N darab csere
/*
void edges_rand1(edges* e, uint64_t N, gsl_rng* r);

//véletlen cserék, a sikertelen próbálkozások is beleszámítanak a cserék számába
//a->b, c->d helyett a->d, c->b
void edges_rand2(edges* e, uint64_t N, gsl_rng* r);
*/

// create a "helper" for faster edge lookups: store the starting index for each addr, making it unnecessary to perform binary
//	search for the whole set, and also return result in O(1) time for addresses with outdeg == 1
int edges_createhelper(edges* e, const idlist* il);

//ID-k generálása N db él alapján
idlist* ids_gen(const edge* e, unsigned int N);

//N db élben az ID-k kicserélése a sorszámukra
int ids_replace(const idlist* il, edge* e, uint64_t N);

#endif


text:
#ifndef INCLUDED_TELEPATH_BLAS_1_H
#define INCLUDED_TELEPATH_BLAS_1_H

#include <complex>
#include <cblas.h>

namespace telepath{

    template< typename T >
    inline auto asum( std::size_t n, const T* x, std::size_t incX = 1 );

    template<> inline auto asum<double>( std::size_t n, const double* x,
        std::size_t incX ){ return cblas_dasum( n, x, incX ); }
    template<> inline auto asum<float>( std::size_t n, const float* x,
        std::size_t incX ){ return cblas_sasum( n, x, incX ); }
    template<> inline auto asum<std::complex<double> >( std::size_t n,
        const std::complex<double>* x, std::size_t incX ){
        return cblas_dzasum( n, (const double*) x, incX ); }
    template<> inline auto asum<std::complex<float> >( std::size_t n,
        const std::complex<float>* x, std::size_t incX ){
        return cblas_scasum( n, (const float*) x, incX ); }

} //namespace tblas

#endif

max_stars_repo_name: tbepler/telepath

max_issues_repo_name: tbepler/telepath

max_forks_repo_name: tbepler/telepath


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

#ifndef __GSL_MATRIX_LONG_DOUBLE_H__
#define __GSL_MATRIX_LONG_DOUBLE_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_long_double.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_long_double * block;
  int owner;
} gsl_matrix_long_double;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_view;

typedef _gsl_matrix_long_double_view gsl_matrix_long_double_view;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_const_view;

typedef const _gsl_matrix_long_double_const_view gsl_matrix_long_double_const_view;

/* Allocation */

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_block (gsl_block_long_double * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_matrix (gsl_matrix_long_double * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_long_double *
gsl_vector_long_double_alloc_row_from_matrix (gsl_matrix_long_double * m,
                                        const size_t i);

GSL_FUN gsl_vector_long_double *
gsl_vector_long_double_alloc_col_from_matrix (gsl_matrix_long_double * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_long_double_free (gsl_matrix_long_double * m);

/* Views */

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_submatrix (gsl_matrix_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_row (gsl_matrix_long_double * m, const size_t i);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_column (gsl_matrix_long_double * m, const size_t j);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_diagonal (gsl_matrix_long_double * m);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_superdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subrow (gsl_matrix_long_double * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_matrix_long_double_subcolumn (gsl_matrix_long_double * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector (gsl_vector_long_double * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector_with_tda (gsl_vector_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_submatrix (const gsl_matrix_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_row (const gsl_matrix_long_double * m,
                            const size_t i);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_column (const gsl_matrix_long_double * m,
                               const size_t j);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_diagonal (const gsl_matrix_long_double * m);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subdiagonal (const gsl_matrix_long_double * m,
                                    const size_t k);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_superdiagonal (const gsl_matrix_long_double * m,
                                      const size_t k);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subrow (const gsl_matrix_long_double * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subcolumn (const gsl_matrix_long_double * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector (const gsl_vector_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector_with_tda (const gsl_vector_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_long_double_set_zero (gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_set_identity (gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_set_all (gsl_matrix_long_double * m, long double x);

GSL_FUN int gsl_matrix_long_double_fread (FILE * stream, gsl_matrix_long_double * m) ;
GSL_FUN int gsl_matrix_long_double_fwrite (FILE * stream, const gsl_matrix_long_double * m) ;
GSL_FUN int gsl_matrix_long_double_fscanf (FILE * stream, gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_fprintf (FILE * stream, const gsl_matrix_long_double * m, const char * format);

GSL_FUN int gsl_matrix_long_double_memcpy(gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);
GSL_FUN int gsl_matrix_long_double_swap(gsl_matrix_long_double * m1, gsl_matrix_long_double * m2);

GSL_FUN int gsl_matrix_long_double_swap_rows(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_swap_columns(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_swap_rowcol(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_long_double_transpose (gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_transpose_memcpy (gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);

GSL_FUN long double gsl_matrix_long_double_max (const gsl_matrix_long_double * m);
GSL_FUN long double gsl_matrix_long_double_min (const gsl_matrix_long_double * m);
GSL_FUN void gsl_matrix_long_double_minmax (const gsl_matrix_long_double * m, long double * min_out, long double * max_out);

GSL_FUN void gsl_matrix_long_double_max_index (const gsl_matrix_long_double * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_long_double_min_index (const gsl_matrix_long_double * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_long_double_minmax_index (const gsl_matrix_long_double * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_long_double_isnull (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_ispos (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_isneg (const gsl_matrix_long_double * m);
GSL_FUN int gsl_matrix_long_double_isnonneg (const gsl_matrix_long_double * m);

GSL_FUN int gsl_matrix_long_double_add (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_sub (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_mul_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_div_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_FUN int gsl_matrix_long_double_scale (gsl_matrix_long_double * a, const double x);
GSL_FUN int gsl_matrix_long_double_add_constant (gsl_matrix_long_double * a, const double x);
GSL_FUN int gsl_matrix_long_double_add_diagonal (gsl_matrix_long_double * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_long_double_get_row(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t i);
GSL_FUN int gsl_matrix_long_double_get_col(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t j);
GSL_FUN int gsl_matrix_long_double_set_row(gsl_matrix_long_double * m, const size_t i, const gsl_vector_long_double * v);
GSL_FUN int gsl_matrix_long_double_set_col(gsl_matrix_long_double * m, const size_t j, const gsl_vector_long_double * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL long double   gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x);
GSL_FUN INLINE_DECL long double * gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const long double * gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
long double
gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
long double *
gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (long double *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const long double *
gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const long double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_LONG_DOUBLE_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs