Search  

---

C:/temp/src/j2k/Beta/Math/Vector/NotDone/Gauss/Smatrix.hpp

Go to the documentation of this file.

#include <map>
#include <vector>
#include <float.h>
#include <math.h>

template <class t> inline void addEq(t &a, const t &b, bool bNeg)
{ if(bNeg) a -= b; else a += b; }

template <class t> inline void CalcEpsilon(t& maxCoef, t &epsilon, const t &d);
template <class t> inline t abs(const t& x) { return(x < 0 ? -x : x); }
template <class t> inline t max(const t& a, const t& b) { return(a > b ? a : b); }
template <class t> inline t min(const t& a, const t& b) { return(a < b ? a : b); }

typedef pair<size_t, size_t> rng;

template <class t> class smIdent;

template<class ta, class tb> ta& svAssign(ta& a, const tb &b)
{
    a.range = b.range;
    ta::iterator ia = a.lower_bound(0);
    ta::iterator iaEnd = a.lower_bound((size_t)-1);
    tb::const_iterator ib = b.lower_bound(b.range.first);
    tb::const_iterator ibEnd = b.lower_bound(b.range.second + 1);
    while(ia != iaEnd || ib != ibEnd)
    {
        if(ia == iaEnd)
        {
            ta::iterator iaNew = a.insert(ia, b.Index(ib), b.Coef(ib));
            if(iaNew != ia)
                iaEnd = ia = ++iaNew;
            ib++;
        }
        else if(ib == ibEnd)
            ia = a.erase(ia);
        else if(a.Index(ia) < b.Index(ib))
            ia = a.erase(ia);
        else if(a.Index(ia) > b.Index(ib))
        {
            ta::iterator iaNew = a.insert(ia, b.Index(ib), b.Coef(ib));
            if(iaNew != ia)
                ia = ++iaNew;
            ib++;
        }
        else
            a.Coef(ia++) = b.Coef(ib++);
    }
    return(a);
}

template<class ta, class tb> ta& svAddEq(ta& a, const tb& b, bool bNegate)
{
    ta::iterator ia = a.lower_bound(min(b.range.first, a.range.second + 1));
    ta::iterator iaEnd = a.lower_bound(a.range.second + 1);
    tb::const_iterator ib = b.lower_bound(b.range.first);
    tb::const_iterator ibEnd = b.lower_bound(b.range.second + 1);

    while(ib != ibEnd)
    {
        if(ia == iaEnd)
        {
            ta::iterator iaNew = a.insert(ia, b.Index(ib), bNegate ? -b.Coef(ib) : b.Coef(ib));
            if(iaNew != ia)
                iaEnd = ia = ++iaNew;
            ib++;
        }
        else if(a.Index(ia) > b.Index(ib))
        {
            ta::iterator iaNew = a.insert(ia, b.Index(ib), bNegate ? -b.Coef(ib) : b.Coef(ib));
            if(iaNew != ia)
                ia = ++iaNew;
            ib++;
        }
        else if(a.Index(ia) < b.Index(ib))
            ia++;
        else
        {
            addEq(a.Coef(ia), b.Coef(ib++), bNegate);
            a.CalcEpsilon(a.Coef(ia));
            ia++;
        }
    }

    a.range.first = min(a.range.first, b.range.first);
    a.range.second = max(a.range.second, b.range.second);

    return(a);
}

template <class t> class smatrix
{
public:
    class vector_map : public vector< map<size_t, size_t> *>
    {
    public:
        vector_map& operator=(const vector_map &a)
        {
            resize(a.size());
            const_iterator ap = a.begin();
            for(iterator p = begin(); p != end(); p++, ap++)
                **p = **ap;
            return(*this);
        }
        ~vector_map() { for(iterator p = begin(); p != end(); p++) delete(*p); }
        void resize(size_t NewSize)
        {
            size_t OldSize = size();
            if(NewSize < OldSize)
            {
                for(iterator p = begin(); p != &(*this)[OldSize]; p++)
                    delete(*p);
                vector< map<size_t, size_t> *>::resize(NewSize);
            }
            else if(NewSize > OldSize)
            {
                vector< map<size_t, size_t> *>::resize(NewSize, NULL);
                for(iterator p = &(*this)[OldSize]; p != end(); p++)
                    *p = new map<size_t, size_t>;
            }
        }
    };

    vector<t> e;
    t epsilon, maxCoef;
    vector_map *iMap, *jMap;
    size_t iDim, jDim, NZ, FirstBlank, FreeSlots;
    bool bOddNumPivots;
    static size_t PrintWidth, PrintPrecision;

    typedef map<size_t, size_t>::iterator iterator;
    typedef map<size_t, size_t>::const_iterator const_iterator;
    typedef pair<size_t, size_t> mapping;

    virtual ~smatrix()
    {
        delete(iMap);
        delete(jMap);
    }
    smatrix() : iMap(NULL), jMap(NULL), maxCoef(0), 
        epsilon(0) { initialize(0, 0, 0); }
    smatrix(const smatrix<t>& a) : iMap(NULL), jMap(NULL) 
    {
        *this = a;
    }
    smatrix(const t* d, size_t iDim, size_t jDim) : iMap(NULL), jMap(NULL) 
    {
        initialize(iDim, jDim, iDim * jDim);
        for(size_t r = 1; r <= iDim; r++)
        {
            for(size_t c = 1; c <= jDim; c++, d++)
            {
                if(*d != 0)
                    insert(r, c, *d);
            }
        }
    }
    smatrix(size_t iDim, size_t jDim, size_t MaxNZ) : iMap(NULL), jMap(NULL) 
    {
        initialize(iDim, jDim, MaxNZ);
    }
    smatrix<t>& operator=(const smatrix<t>& a)
    {
        e = a.e;
        epsilon = a.epsilon;
        maxCoef = a.maxCoef;

        delete(iMap);
        iMap = new vector_map;
        *iMap = *a.iMap;

        delete(jMap);
        jMap = new vector_map;
        *jMap = *a.jMap;

        NZ = a.NZ;
        iDim = a.iDim;
        jDim = a.jDim;
        FirstBlank = a.FirstBlank;
        bOddNumPivots = a.bOddNumPivots;

        return(*this);
    }
    void initialize(size_t MaxNZ)
    {
        delete(iMap);
        iMap = new vector_map;
        delete(jMap);
        jMap = new vector_map;
        e.reserve(MaxNZ);
        iMap->resize(iDim + 1);
        jMap->resize(jDim + 1);
        FirstBlank = (size_t)-1;
        NZ = 0;
        FreeSlots = 0;
        bOddNumPivots = false;
        epsilon = 0;
        maxCoef = 0;
    }
    void initialize(size_t iDim, size_t jDim, size_t MaxNZ)
    {
        this->iDim = iDim;
        this->jDim = jDim;
        initialize(MaxNZ);
    }
    void erase(size_t pos)
    {
        e[pos] = FirstBlank;
        FirstBlank = pos;
        NZ--;
        FreeSlots++;
    }
    size_t SlackSpace() const
    {
        return(FreeSlots + (e.size() - NZ));
    }
    double Sparsity() const
    {
        return((double)NZ / (iDim * jDim));
    }
    void reserve(size_t space)
    {
        size_t slack = SlackSpace();
        if(slack < space)
            e.reserve(space - FreeSlots);
    }
    bool insert(size_t i, size_t j, const t& d, iterator &iIter, iterator &jIter, bool bZeroInsert = false)
    {
        if(bZeroInsert == false && abs(d) <= epsilon)
            return(false);

        size_t pos;
        if(FirstBlank != (size_t)-1)
        {
            pos = FirstBlank;
            FirstBlank = (size_t)e[pos];
            FreeSlots--;
        }
        else
        {
            if(e.capacity() == e.size())
                e.reserve(e.size() + (e.size() / 10));
            pos = e.size();
            e.insert(e.end());
        }
        NZ++;

        e[pos] = d;
        CalcEpsilon(maxCoef, epsilon, d);

        if(i > iDim)
        {
            iDim = i;
            iMap->resize(iDim + 1);
        }
        jIter = (*iMap)[i]->insert(mapping(j, pos)).first;

        if(j > jDim)
        {
            jDim = j;
            jMap->resize(jDim + 1);
        }
        iIter = (*jMap)[j]->insert(mapping(i, pos)).first;

        return(true);
    }
    iterator insert(size_t i, size_t j, const t& d)
    {
        iterator iIter, jIter;
        if(insert(i, j, d, iIter, jIter))
            return(iIter);
        return(IterForRow(0, 0));
    }
    iterator ZeroInsert(size_t i, size_t j)
    {
        iterator iIter, jIter;
        t d(0);
        insert(i, j, d, iIter, jIter, true);
        return(iIter);
    }
    smatrix<t>& Transpose() 
    { 
        swap(iMap, jMap); 
        swap(iDim, jDim);
        return(*this);
    }
    template<class ta, class tb> static t svMultiply(const ta& a, const tb& b)
    {
        size_t first = max(a.range.first, b.range.first);
        size_t last = min(a.range.second, b.range.second);

        ta::const_iterator ia = a.lower_bound(first);
        ta::const_iterator iaEnd = a.lower_bound(last + 1);

        tb::const_iterator ib = b.lower_bound(first);
        tb::const_iterator ibEnd = b.lower_bound(last + 1);

        t d(0);
        while(ia != iaEnd && ib != ibEnd) 
        {
            if(a.Index(ia) < b.Index(ib))
                ia++;
            else if(a.Index(ia) > b.Index(ib))
                ib++;
            else
            {
                addEq(d, a.Coef(ia) * b.Coef(ib), false);
                ia++;
                ib++;
            }
        }
        return(d);
    }
    class svector : public map<size_t, t>
    {
    public:
        rng range;
        t epsilon, maxCoef;

        svector(const t &ep) : epsilon(ep), maxCoef(0) { }
        t GetEpsilon() const { return(epsilon); }
        template<class tb> t operator*(const tb &b)
        {
            epsilon *= b;
            maxCoef *= b;
            return(svMultiply(*this, b));
        }
        svector& operator=(const class const_svector_ref& a);
        t Coef(const_iterator p) const
        {
            return(p->second);
        }
        static size_t Index(iterator p) { return(p->first); }
        static size_t Index(const_iterator p) { return(p->first); }
        static t& Coef(iterator p)  { return(p->second); }
        void CalcEpsilon(const t &d)
        {
            sm::CalcEpsilon(maxCoef, epsilon, d);
        }
        iterator insert(iterator p, size_t x, const t& d)
        {
            CalcEpsilon(d);
            return(map<size_t, t>::insert(p, pair<size_t, t>(x, d)));
        }
        t& operator()(size_t n)
        {
            return(operator[](n)->second);
        }
        t operator()(size_t n) const
        {
            const_iterator p = find(n);
            return(p == end() ? 0 : p->second);
        }
    };
    class const_svector_ref
    {
    public:
        const smatrix<t> *cm;
        rng range;
        size_t n;
        bool bColVect;

        typedef map<size_t, size_t>::const_iterator const_iterator;

        const_svector_ref(const smatrix<t> *pm, const rng& i, size_t j) :
            cm(pm), range(i), n(j), bColVect(true) { }
        const_svector_ref(const smatrix<t> *pm, size_t i, const rng& j) :
            cm(pm), range(j), n(i), bColVect(false) { }
        const map<size_t, size_t> *GetMap() const
        {
            if(bColVect && n < cm->jMap->size())
                return((*cm->jMap)[n]);
            if(!bColVect && n < cm->iMap->size())
                return((*cm->iMap)[n]);
            return((*cm->iMap)[0]);
        }
        const_iterator lower_bound(size_t i) const
        {
            return(GetMap()->lower_bound(i));
        }
        svector operator*(const t& b)
        {
            svector c(GetEpsilon() * b);
            c.range = range;
            const_iterator i = lower_bound(range.first);
            const_iterator iEnd = lower_bound(range.second + 1);
            for(; i != iEnd; i++) 
                c.insert(c.end(), i->first, Coef(i) * b);
            return(c);
        }
        t GetEpsilon() const { return(cm->epsilon); }
        static size_t Index(const_iterator p) { return(p->first); }
        t Coef(const_iterator p) const { return(cm->e[p->second]); }
        t operator()(size_t n) const
        {
            const map<size_t, size_t> *Map = GetMap();
            const_iterator p = Map->find(n);
            return(p == Map->end() ? 0 : cm->e[p->second]);
        }
        bool operator!=(const const_svector_ref& b)
        {
            if(range.first != b.range.first || range.second != b.range.second) return(true);
            const_iterator ia = lower_bound(range.first);
            const_iterator iaEnd = lower_bound(range.second + 1);
            const_iterator ib = b.lower_bound(b.range.first);
            const_iterator ibEnd = b.lower_bound(b.range.second + 1);
            for(; ia != iaEnd && ib != ibEnd; ia++, ib++)
            {
                if(Index(ia) != b.Index(ib)) return(true);
                if(Coef(ia) != b.Coef(ib)) return(true);
            }
            return(ia != iaEnd || ib != ibEnd);
        }
        template<class tb> t operator*(const tb &b)
        {
            return(smatrix<t>::svMultiply(*this, b));
        }
        template<class tb> svector operator+(const tb &b) const
        {
            svector c(max(GetEpsilon(), b.GetEpsilon()));
            svAssign(c, *this);
            return(svAddEq(c, b, false));
        }
        template<class tb> svector operator-(const tb &b) const
        {
            svector c(max(GetEpsilon(), b.GetEpsilon()));
            svAssign(c, *this);
            return(svAddEq(c, b, true));
        }
    };
    class svector_ref : public const_svector_ref
    {
    public:
        smatrix<t> *m;
        typedef map<size_t, size_t>::iterator iterator;

        svector_ref(smatrix<t> *pm, const rng& i, size_t j) :
            m(pm), const_svector_ref(pm, i, j) { }
        svector_ref(smatrix<t> *pm, size_t i, const rng& j) :
            m(pm), const_svector_ref(pm, i, j) { }
        map<size_t, size_t> *GetMap()
        {
            if(bColVect && n < m->jMap->size())
                return((*m->jMap)[n]);
            if(!bColVect && n < m->iMap->size())
                return((*m->iMap)[n]);
            return((*m->iMap)[0]);
        }
        iterator lower_bound(size_t i) 
        {
            return(GetMap()->lower_bound(i));
        }
        static size_t Index(iterator p) { return(p->first); }
        t& Coef(iterator p) 
        {
            return(m->e[p->second]);
        }
        svector_ref operator=(const const_svector_ref &b)
        {
            return(svAssign(*this, b));
        }
        svector_ref& operator=(const svector &b)
        {
            return(svAssign(*this, b));
        }
        iterator insert(iterator p, size_t x, const t& d)
        {
            iterator iIter, jIter;
            if(bColVect)
            {
                iIter = p;
                m->insert(x, n, d, iIter, iIter);
                return(iIter);
            }
            jIter = p;
            m->insert(n, x, d, iIter, jIter);
            return(jIter);
        }
        iterator erase(iterator p)
        {
            m->erase(p->second);
            if(bColVect)
            {
                (*m->iMap)[p->first]->erase(n);
                return((*m->jMap)[n]->erase(p));
            }
            (*m->jMap)[p->first]->erase(n);
            return((*m->iMap)[n]->erase(p));
        }
        t& operator()(size_t n)
        {
            const map<size_t, size_t> *Map = GetMap();
            const_iterator p = Map->find(n);
            return(p == Map->end() ? insert(n, t(0)) : cm->e[p->second]);
        }
        void CalcEpsilon(const t &d)
        {
            sm::CalcEpsilon(m->maxCoef, m->epsilon, d);
        }
        svector_ref& operator+=(const const_svector_ref &b)
        {
            return(svAddEq(*this, b, false));
        }
        svector_ref& operator-=(const const_svector_ref &b)
        {
            return(svAddEq(*this, b, true));
        }
        svector_ref& operator+=(const svector &b)
        {
            return(svAddEq(*this, b, false));
        }
        svector_ref& operator-=(const svector &b)
        {
            return(svAddEq(*this, b, true));
        }
        svector_ref& operator*=(const t& d)
        {
            iterator i = lower_bound(range.first);
            iterator iEnd = lower_bound(range.second + 1);
            for(; i != iEnd; i++) 
                Coef(i) *= d;
            return(*this);
        }
    };
    t operator()(size_t i, size_t j) const
    {
        if(i > iDim || j > jDim) return(0);
        const map<size_t, size_t> *r = (*iMap)[i];
        map<size_t, size_t>::const_iterator p = r->find(j);
        if(p == r->end())
            return(0);
        return(e[p->second]);
    }
    class coefRef
    {
    public:
        smatrix<t> *m;
        size_t i, j;

        coefRef(smatrix<t> *m, size_t i, size_t j)
        {
            this->m = m;
            this->i = i;
            this->j = j;
        }
        coefRef& operator=(const t &d)
        {
            if(abs(d) > m->epsilon)
            {
                if(i > m->iDim || j > m->jDim)
                    m->insert(i, j, d);
                else
                {
                    const map<size_t, size_t> *r = (*m->iMap)[i];
                    map<size_t, size_t>::const_iterator p = r->find(j);
                    if(p == r->end())
                        m->insert(i, j, d);
                    else
                    {
                        m->e[p->second] = d;
                        CalcEpsilon(m->maxCoef, m->epsilon, d);
                    }
                }
            }
            else if(i <= m->iDim && j <= m->jDim)
                m->erase(i, j);
            return(*this);
        }
        operator t() const
        {
            if(i > m->iDim || j > m->jDim) return(0);
            const map<size_t, size_t> *r = (*m->iMap)[i];
            map<size_t, size_t>::const_iterator p = r->find(j);
            if(p == r->end())
                return(0);
            return(m->e[p->second]);
        }
    };
    coefRef operator()(size_t i, size_t j)
    {
        return(coefRef(this, i, j));
    }
    void erase(size_t i, size_t j)
    {
        map<size_t, size_t> *r = (*iMap)[i];
        map<size_t, size_t>::iterator p = r->find(j);
        if(p != r->end())
        {
            erase(p->second);
            r->erase(p);
            (*jMap)[j]->erase(i);
        }
    }
    t Coef(size_t i, size_t j) const { return((*this)(i, j)); }
    coefRef Coef(size_t i, size_t j) { return((*this)(i, j)); }
    t Coef(const_iterator p) const { return(e[p->second]); }
    t& Coef(const_iterator p) { return(e[p->second]); }
    static size_t Index(const_iterator p) { return(p->first); }
    static size_t Index(iterator p) { return(p->first); }
    static size_t CoefPtr(const_iterator p) { return(p->second); }
    static size_t CoefPtr(iterator p) { return(p->second); }

    svector_ref operator()(size_t i, const rng& range) 
        { return(svector_ref(this, i, range)); }
    svector_ref operator()(const rng& range, size_t j) 
        { return(svector_ref(this, range, j)); }
    const_svector_ref operator()(size_t i, const rng& range) const
        { return(const_svector_ref(this, i, range)); }
    const_svector_ref operator()(const rng& range, size_t j) const
        { return(const_svector_ref(this, range, j)); }

    svector_ref row(size_t i) 
        { return(svector_ref(this, i, rng(1, jDim))); }
    svector_ref row(size_t i, const rng& range) 
        { return(svector_ref(this, i, range)); }
    svector_ref col(size_t j) 
        { return(svector_ref(this, rng(1, iDim), j)); }
    svector_ref col(const rng& range, size_t j) 
        { return(svector_ref(this, range, j)); }
    const_svector_ref row(size_t i) const
        { return(const_svector_ref(this, i, rng(1, jDim))); }
    const_svector_ref row(size_t i, const rng& range) const
        { return(const_svector_ref(this, i, range)); }
    const_svector_ref col(size_t j) const
        { return(const_svector_ref(this, rng(1, iDim), j)); }
    const_svector_ref col(const rng& range, size_t j) const
        { return(const_svector_ref(this, range, j)); }

    iterator IterFor(size_t n, size_t Start, const vector_map *p)
    {
        map<size_t, size_t> *Map = (n > p->size() ? (*p)[0] : (*p)[n]);
        return(Start <= 1 ? Map->begin() : Map->lower_bound(Start));
    }
    const_iterator IterFor(size_t n, size_t Start, const vector_map *p) const
    {
        const map<size_t, size_t> *Map = (n > p->size() ? (*p)[0] : (*p)[n]);
        return(Start <= 1 ? Map->begin() : Map->lower_bound(Start));
    }
    iterator EndOf(size_t n, vector_map *p)
    {
        map<size_t, size_t> *Map = (n > p->size() ? (*p)[0] : (*p)[n]);
        return(Map->end());
    }
    const_iterator EndOf(size_t n, const vector_map *p) const
    {
        const map<size_t, size_t> *Map = (n > p->size() ? (*p)[0] : (*p)[n]);
        return(Map->end());
    }
    iterator IterForRow(size_t i, size_t jStart = 1)
        { return(IterFor(i, jStart, iMap)); }
    const_iterator IterForRow(size_t i, size_t jStart = 1) const
        { return(IterFor(i, jStart, iMap)); }
    const_iterator EndOfRow(size_t i) const
        { return(EndOf(i, iMap)); }
    iterator EndOfRow(size_t i)
        { return(EndOf(i, iMap)); }
    iterator IterForCol(size_t j, size_t iStart = 1)
        { return(IterFor(j, iStart, jMap)); }
    const_iterator IterForCol(size_t j, size_t iStart = 1) const
        { return(IterFor(j, iStart, jMap)); }
    const_iterator EndOfCol(size_t j) const
        { return(EndOf(j, jMap)); }
    iterator EndOfCol(size_t j)
        { return(EndOf(j, jMap)); }

    smatrix<t> operator+(const smatrix<t> &b) const
    {
        smatrix<t> c;
        c.initialize(max(iDim, b.iDim), max(jDim, b.jDim), NZ + b.NZ);
        for(size_t i = 1; i <= iDim; i++)
            c.row(i) = row(i) + b.row(i);
        return(c);
    }
    smatrix<t> operator-(const smatrix<t> &b) const
    {
        smatrix<t> c;
        c.initialize(max(iDim, b.iDim), max(jDim, b.jDim), NZ + b.NZ);
        for(size_t i = 1; i <= iDim; i++)
            c.row(i) = row(i) - b.row(i);
        return(c);
    }
    smatrix<t> operator*(const smatrix<t>& b) const
    {
        smatrix<t> c;

        c.iDim = iDim;
        c.jDim = b.jDim;
        size_t MaxNZ = max(NZ, b.NZ);
        if(c.iDim == 1 && c.jDim == 1) MaxNZ = 1;
        else if(c.iDim == 1 || c.jDim == 1) MaxNZ = max(c.iDim, c.jDim);
        c.initialize(MaxNZ);

        for(size_t i = 1; i <= iDim; i++)
        {
            const_svector_ref r = row(i);
            for(size_t j = 1; j <= b.jDim; j++)
            {
                t d = svMultiply(r, b.col(j));
                c.insert(i, j, d);
            }
        }
        return(c);
    }
    bool operator==(const smatrix<t> &b) const
    {
        if(iDim != b.iDim || jDim != b.jDim || NZ != b.NZ) return(false);
        for(size_t i = 1; i <= iDim; i++)
        {
            if(row(i) != b.row(i))
                return(false);
        }
        return(true);
    }
    bool operator!=(const const smatrix<t>& b)
    {
        return(operator==(b) == false);
    }
    bool operator!() { return(NZ == 0); }
    void PivotMap(size_t a, size_t b, vector_map *Source, vector_map &Dest)
    {
        bOddNumPivots = (bOddNumPivots == false);
        if(a != b)
        {
            const_iterator ia = IterFor(a, 1, Source);
            const_iterator ib = IterFor(b, 1, Source);
            while(ia != EndOf(a, Source) || ib != EndOf(b, Source))
            {
                if(ia == EndOf(a, Source))
                {
                    Dest[Index(ib)]->erase(b);
                    Dest[Index(ib)]->insert(mapping(a, CoefPtr(ib)));
                    ib++;
                }
                else if(ib == EndOf(b, Source))
                {
                    Dest[Index(ia)]->erase(a);
                    Dest[Index(ia)]->insert(mapping(b, CoefPtr(ia)));
                    ia++;
                }
                else if(Index(ia) > Index(ib))
                {
                    Dest[Index(ib)]->erase(b);
                    Dest[Index(ib)]->insert(mapping(a, CoefPtr(ib)));
                    ib++;
                }
                else if(Index(ia) < Index(ib))
                {
                    Dest[Index(ia)]->erase(a);
                    Dest[Index(ia)]->insert(mapping(b, CoefPtr(ia)));
                    ia++;
                }
                else
                {
                    map<size_t, size_t> &p = *Dest[Index(ia)];
                    swap(p[a], p[b]);
                    ia++;
                    ib++;
                }
            }
        }
        swap((*Source)[a], (*Source)[b]);
    }
    bool Pivot(size_t k, size_t &MaxIdx, t &MaxCoef, vector_map *Source, vector_map &Dest)
    {
        iterator p = Dest[k]->lower_bound(k);
        if(p == Dest[k]->end())
            return(false);
        MaxIdx = Index(p);
        MaxCoef = Coef(p);
        for(p++; p != Dest[k]->end(); p++)
        {
            if(abs(MaxCoef) < abs(Coef(p)))
            {
                MaxCoef = Coef(p);
                MaxIdx = Index(p);
            }
        }
        if(k != MaxIdx)
            PivotMap(k, MaxIdx, Source, Dest);

        return(true);
    }
    void PivotRow(size_t k, size_t MaxIdx, smatrix<t> *pMatix = NULL)
    {
        PivotMap(k, MaxIdx, iMap, *jMap);
        if(pMatix != NULL && k != MaxIdx)
            pMatix->PivotMap(k, MaxIdx, pMatix->iMap, *pMatix->jMap);
    }
    bool PivotRow(size_t k, t& MaxCoef, smatrix<t> *pMatix = NULL)
    {
        size_t MaxIdx;
        if(Pivot(k, MaxIdx, MaxCoef, iMap, *jMap))
        {
            if(pMatix != NULL && k != MaxIdx)
                pMatix->PivotMap(k, MaxIdx, pMatix->iMap, *pMatix->jMap);
            return(true);
        }
        return(false);
    }
    void PivotCol(size_t k, size_t MaxIdx, smatrix<t> *pMatix = NULL)
    {
        PivotMap(k, MaxIdx, jMap, *iMap);
        if(pMatix != NULL && k != MaxIdx)
            pMatix->PivotMap(k, MaxIdx, pMatix->jMap, *pMatix->iMap);
    }
    bool PivotCol(size_t k, t& MaxCoef, smatrix<t> *pMatix = NULL)
    {
        size_t MaxIdx;
        if(Pivot(k, MaxCoef, jMap, *iMap, L))
        {
            if(pMatix != NULL && k != MaxIdx)
                pMatix->PivotMap(k, MaxIdx, pMatix->jMap, *pMatix->iMap);
            return(true);
        }
        return(false);
    }
    bool GetDiag(size_t k, t& Diag, bool bWithPartialPivot = false, smatrix<t> *pMatix = NULL)
    {
        if(bWithPartialPivot)
        {
            if(PivotRow(k, Diag, pMatix) == false)
                return(false);
        }
        else
            Diag = Coef(k, k);
        if(abs(Diag) < epsilon)
            return(false);
        return(true);
    }
    bool FwdElim(bool bWithPartialPivot, bool bErase = true)
    {
        for(size_t k = 1; k < iDim; k++)
        {
            t Diag;
            if(!GetDiag(k, Diag, bWithPartialPivot)) return(false);
            Diag = t(-1.0) / Diag;
            iterator i = IterForCol(k, k + 1);
            while(i != EndOfCol(k))
            {
                t Scaler = Coef(i) * Diag;
                size_t r = Index(i++);
                row(r) += row(k, rng(k + 1, jDim)) * Scaler;
                if(bErase) erase(r, k);
            }
        }
        t Diag = Coef(k, k);
        if(abs(Diag) < epsilon)
            return(false);
        return(true);
    }
    bool BackElim(bool bErase = true, bool bScaleDiag = true)
    {
        for(size_t k = iDim; k >= 1; k--)
        {
            t Diag;
            if(!GetDiag(k, Diag)) return(false);
            Diag = t(-1.0) / Diag;
            iterator i = IterForCol(k);
            while(Index(i) < k)
            {
                t Scaler = Coef(i) * Diag;
                size_t r = Index(i++);
                row(r) += row(k) * Scaler;
                if(bErase) erase(r, k);
            }
            if(bScaleDiag) row(k) *= -Diag;
        }
        return(true);
    }
    bool SolveWithGaussElim(const smatrix<t> &m, bool bWithPartialPivot = true)
    {
        *this = m;
        if(FwdElim(bWithPartialPivot) == false) return(false);
        return(BackElim());
    }
    smatrix<t> SolveWithGaussElim(bool bWithPartialPivot = true) const
    {
        smatrix<t> m;
        if(m.SolveWithGaussElim(*this, bWithPartialPivot) == false)
            return(smatrix<t>());
        return(m);
    }
    t det(bool bWithPartialPivot)
    {
        if(iDim != jDim) return(t(0));
        FwdElim(bWithPartialPivot, true);
        t d(bOddNumPivots ? -1.0 : 1.0);
        for(size_t k = 1; k <= iDim; k++)
            d *= Coef(k, k);
        return(d);
    }
    int lnScaledDet(t& lnScale, bool bWithPartialPivot)
    {
        lnScale = t(0);

        if(iDim != jDim) return(t(0));
        FwdElim(bWithPartialPivot, true);
        int sign = (bOddNumPivots ? -1 : 1);

        for(size_t k = 1; k <= iDim; k++)
        {
            t d;
            if(!GetDiag(k, d, false)) return(false);
            if(d > 0) 
                lnScale += log(d);
            else
            {
                sign *= -1;
                lnScale += log(-d);
            }
        }
        return(sign);
    }
    int lnScaledDet(const smatrix<t>& m, t &lnScale, bool bWithPartialPivot)
    {
        *this = m;
        return(lnScaledDet(lnScale, bWithPartialPivot));
    }
    int lnScaledDet(t &lnScale, bool bWithPartialPivot = true) const
    {
        smatrix<t> m;
        return(m.lnScaledDet(*this, lnScale, bWithPartialPivot));
    }
    bool FwdElim(smatrix<t> &m, bool bWithPartialPivot)
    {
        for(size_t k = 1; k < iDim; k++)
        {
            t Diag;
            if(!GetDiag(k, Diag, bWithPartialPivot, &m)) return(false);
            Diag = t(-1.0) / Diag;
            iterator i = IterForCol(k, k + 1);
            while(i != EndOfCol(k))
            {
                t Scaler = Coef(i) * Diag;
                size_t r = Index(i++);
                row(r) += row(k, rng(k + 1, jDim)) * Scaler;
                m.row(r) += m.row(k) * Scaler;
            }
        }
        t Diag = Coef(k, k);
        if(abs(Diag) < epsilon)
            return(false);
        return(true);
    }
    void RecalcMaxCoef()
    {
        maxCoef = t(0);
        epsilon = 0;
        for(size_t i = 1; i < iDim; i++)
        {
            iterator j = IterForRow(i);            
            for(; j != EndOfRow(i); j++)
                CalcEpsilon(maxCoef, epsilon, e[j->second]);
        }
    }
    bool BackElim(smatrix<t> &m)
    {
        for(size_t k = iDim; k >= 1; k--)
        {
            t Diag;
            if(!GetDiag(k, Diag)) return(false);
            Diag = t(-1.0) / Diag;
            iterator i = IterForCol(k);
            while(Index(i) < k)
            {
                t Scaler = Coef(i) * Diag;
                size_t r = Index(i++);
                row(r) += row(k, rng(k, k)) * Scaler;
                m.row(r) += m.row(k) * Scaler;
            }
            m.row(k) *= -Diag;
        }
        m.RecalcMaxCoef();
        return(true);
    }
    smatrix<t> Solve(const smatrix<t> &L, const smatrix<t> &U) const
    {
        smatrix<t> m, c(*this);
        m = L;
        if(m.FwdElim(c, false))
        {
            m = U;
            if(m.BackElim(c))
                return(c);
        }
        return(smatrix<t>());
    }
    smatrix<t> Inverse(bool bWithPartialPivot) const
    {
        smIdentity<t> c(iDim);
        c.epsilon = epsilon;
        smatrix<t> m(*this);

        if(m.FwdElim(c, bWithPartialPivot))
        {
            if(m.BackElim(c))
                return(c);
        }
        return(smatrix<t>());
    }
    t FrobeniousNormal() const
    {
        if(NZ == 0) return(0);
        t d(0);
        for(vector<t>::const_iterator i = e.begin(); i != e.end(); i++)
            d += *i * *i;
        return(sqrt(d));
    }
    smatrix<t>& operator*=(const t& d)
    {
        for(vector<t>::iterator i = e.begin(); i != e.end(); i++)
            *i *= d;
        return(*this);
    }
    smatrix<t> operator*(const t& d) const
    {
        smatrix<t> c(*this);
        return(c *= d);
    }
    void PrintListing(ostream &os) const
    {
        os << "Matrix" << endl;
        os << *this << endl;
        os << "vector<t> e" << endl;
        for(size_t n = 0; n < e.size(); n++)
            os << e[n] << " ";
        os << endl << endl;

        for(size_t r = 1; r <= iDim; r++)
        {
            os << "iMap[" << r << "] ";
            for(const_iterator i = IterForRow(r); i != EndOfRow(r); i++)
                os << "(" << i->first << ", " << i->second << ") ";
            os << endl;
        }
        os << endl;
        for(size_t c = 1; c <= jDim; c++)
        {
            os << "jMap[" << c << "] ";
            for(const_iterator j = IterForCol(c); j != EndOfCol(c); j++)
                os << "(" << j->first << ", " << j->second << ") ";
            os << endl;
        }
    }
};

template <class t> smatrix<t> operator*(const t &d, const smatrix<t> &m)
{
    return(m * d);
}

template <class t> class smIdentity : public smatrix<t>
{
public:
    smIdentity(size_t dim) : smatrix<t>(dim, dim, dim)
    {
        for(size_t i = 1; i <= dim; i++)
            insert(i, i, t(1));
    }
};

template <class t> inline t det(const smatrix<t>& m, bool bWithPartialPivot = true)
{
    smatrix<t> c(m);
    return(c.det(bWithPartialPivot));
}

template <class t> inline ostream& operator<<(ostream& os, const smatrix<t>::svector& v)
{
    for(size_t n = v.range.first; n <= v.range.second; n++)
    {
        os.width(smatrix<t>::PrintWidth);
        os.precision(smatrix<t>::PrintPrecision);
        os << v(n);
    }
    return(os);
}

template <class t> inline ostream& operator<<(ostream& os, const smatrix<t>::const_svector_ref& v)
{
    for(size_t n = v.range.first; n <= v.range.second; n++)
    {
        os.width(smatrix<t>::PrintWidth);
        os.precision(smatrix<t>::PrintPrecision);
        os << v(n);
    }
    return(os);
}

template <class t> inline ostream& operator<<(ostream& os, const smatrix<t>& m)
{
    if(m.NZ > 0)
    {
        if(m.iDim < 20)
        {
            for(size_t i = 1; i <= m.iDim; i++)
                os << m(i, rng(1, m.jDim)) << endl;
        }
        else
        {
            os << m.iDim << " x " << m.jDim << ", NZ = " << m.NZ << endl;
            os << "Max coef = " << m.maxCoef << ", epsilon = " << m.epsilon << endl;
        }
    }
    else
        os << "<null>";
    return(os);
}

template <class t> inline t Residual(const smatrix<t> &a, const smatrix<t> &b)
{
    return((a - b).maxCoef);
}

inline void CalcEpsilon(double &maxCoef, double &epsilon, const double &d)
{
    double ep = abs(d);
    if(maxCoef < ep) {
        maxCoef = ep;
        epsilon = max(ep * DBL_EPSILON, epsilon); }
}

---