---

matrixSP.h

Go to the documentation of this file.

#ifndef MATRIXSP_H
#define MATRIXSP_H

// project includes
#include    "matrix.h"
// standard includes
#include    <vector>
#include    <map>
#include    <valarray>
#ifdef USE_NAMESPACE_LINAL
namespace LinAl {
#endif
/*
 * Classes
 */
class MatrixSP;

#define REF_SP(A,I,J)   ((*(((A).D)->v))[(J)][(I)])

#define FTN_OFFSET  1

class MatrixSP;
istream& operator>> (istream& s, MatrixSP& M) ;
ostream& operator<< (ostream& s, const MatrixSP& M) ;
class MatrixSP:public Matrix<double> {
    friend istream& operator>> (istream& s, MatrixSP& M) ;
    friend ostream& operator<< (ostream& s, const MatrixSP& M) ;
    typedef map<int,double>  vectorsp ;
public:
    enum status {open=0,closed=1};
protected:
    class MrepSP {
    public:
        status              s;
        valarray<double>*   m;
        valarray<int>*      ia;
        valarray<int>*      izh;
        vector<vectorsp>*   v ;
        int                 ncol ;
        int                 nref ;
        MrepSP (const int nc=0, const status S=open) : 
                s(S),m(0),ia(0),izh(0),v(0),ncol(nc) {
            nref=1;
            if (S==open) v=new vector<vectorsp>(nc);
            else {
                izh=new valarray<int>(nc);
                ia=new valarray<int>();
                m=new valarray<double>();
            }
        }
        MrepSP (const int nc, const int nel, const status _s=closed) :
                s(_s),m(0),ia(0),izh(0),v(0),ncol(nc) {
            nref=1;
            if (_s==open) v=new vector<vectorsp>(nc);
            else {
                izh=new valarray<int>(nc);
                ia=new valarray<int>(nel);
                m=new valarray<double>(nel);
            }
        }
        ~MrepSP () {
            if (s==closed) {
                delete m;
                delete ia;
                delete izh;
            }
            else
                delete v;
        }
        MrepSP* get_own_copy()
        {
            if (nref==1) return this ;
            nref-- ;
            return new MrepSP(*this) ;
        }
        void assign(const double value) {
            if (s==open) {
                for (int i=0;i<ncol;++i) {
                    vectorsp::iterator j=(*v)[i].begin(),je=(*v)[i].end();
                    while (j!=je) {
                        j->second=value;
                        j++;
                    }
                }
            }
            else {
                fill (&(*m)[0],&(*m)[m->size()],value);
            }
        }
/*
 *      
 *      void assign(const int nr, const T* value)
 *      {
 *          if (nr>(int)m.size()) m.resize(nr) ;
 *          for (int i=0;i<nr;++i) m[i]=value[i] ;
 *      }
 *      
 *      void assign(const int nr, const int nc, const T** value)
 *      {
 *          if ((nr*nc)>(int)m.size()) m.resize(nr*nc) ;
 *          int k=0;
 *          for (int j=0;j<nc;++j) 
 *              for (int i=0;i<nr;++i) m[k++]=value[i][j] ;                 
 *      }
 */
        int size() const {
            if (s==closed) return m->size();
            else {
                int is=0;
                for (int i=0;i<ncol;++i) is += (*v)[i].size() ;
                return is;
            }
        }
        void set_status(status ns=open) {
            int ftn_offset=FTN_OFFSET;
            if (ns==s) return ;
            if (ns==open)
            {
                if (!v) v=new vector<vectorsp>(ncol);
                int k=0,nel=m->size();
                for (int j=0;j<ncol;++j)
                {
                    int nr=(j<(ncol-1)?((*izh)[j+1]-(*izh)[j]):(nel-(*izh)[j]));
                    for (int i=0;i<nr;++i)
                    {
                        int key; double value;
                        key=(*ia)[k]-ftn_offset;
                        value=(*m)[k];
                        (*v)[j][key]=value;
                        k++;
                    }
                }
                delete izh; izh=0;
                delete ia; ia=0;
                delete m; m=0;
                s=open;
            }
            else if (ns==closed)
            {
                typedef map<int,double>::const_iterator ci;
                int nel=0;
                int k=0;
                int j;
                nel=size();
                if (!izh) izh=new valarray<int>(ncol);
                if (!ia) ia=new valarray<int>(nel);
                if (!m) m=new valarray<double>(nel);
                for (j=0;j<ncol;++j)
                {
                    (*izh)[j]=k+ftn_offset;
                    for (ci i=(*v)[j].begin();i!=(*v)[j].end();++i)
                    {
                        (*ia)[k] = i->first+ftn_offset ;
                        (*m)[k] = i->second ;
                         k++;
                    }
                }
                delete v; v=0;
                s=closed;
            }
            return;
        }
        status get_status() {return s;}
    private:
        MrepSP (const MrepSP& D)
        : s(D.s),m(0),ia(0),izh(0),v(0),ncol(D.ncol)
        {
            if (s==open) v=new vector<vectorsp>(*(D.v));
            else {
                izh=new valarray<int>(*(D.izh));
                ia=new valarray<int>(*(D.ia));
                m=new valarray<double>(*(D.m));
            }           
            nref=1;
        }
        MrepSP& operator= (const MrepSP&) ;
    };
    MrepSP *D ;
public:
    MatrixSP (const int nc=0,const status S=open) {D=new MrepSP(nc,S);}
    MatrixSP (const int _ncol, const int _nel, const status _s=closed) {D=new MrepSP(_ncol,_nel,_s);}
    MatrixSP (const MatrixSP& B) {B.D->nref++;D=B.D;}
    virtual ~MatrixSP () {if (--D->nref==0) delete D;} ;
    virtual MatrixSP& operator = (const MatrixSP& B)
    {
        if (&B != this) {
            B.D->nref++ ;
            if (--D->nref==0) delete D ;
            D=B.D ;
        }
        return *this ;
    }
    virtual MatrixSP& operator = (const double value)
    {
        D=D->get_own_copy();
        D->assign(value);
        return *this;
    }
    void get_own_copy() {D=D->get_own_copy();}
    void set_open() {D->set_status(open);}
    void set_closed() {D->set_status(closed);}
    bool is_open() {return D->s==open?true:false;}
    void resize (const int nc, const double value=0) {cout << "not yet implemented" <<endl;}
    //void reserve (const int nr_x_nc) { D->m.reserve(nr_x_nc) ;}
    virtual int nr (const int i=0) const {
        if (D->s==closed) {
            if (i==D->ncol-1) return D->ia->size()-(*D->izh)[i];
            else return (*D->izh)[i+1]-(*D->izh)[i];
        }
        else {
            return (*D->v)[i].size();
        }
    } ;
    virtual int nc (const int i=0) const {return D->ncol;} ;
    virtual int size() const {return D->size();}
    virtual bool empty() const {return D->v->empty();}
    //virtual int   capacity() const {return D->m.capacity();}
    double* get_m() const {assert(D->s==closed);return (&(*(D->m))[0]);}
    int* get_ia() const {assert(D->s==closed);return (&(*(D->ia))[0]);}
    int* get_izh() const {assert(D->s==closed);return (&(*(D->izh))[0]);}
protected:
    vectorsp& operator [] (const int i) const {assert(D->s==open);return (*D->v)[i];}
public:
    virtual Ref<double> operator () (const int i, const int j=0) {assert(D->s==open);return Ref<double>(*this,i,j);}
    virtual double operator () (const int i, const int j=0) const {assert(D->s==open);return REF_SP(*this,i,j);}
protected:
    double* ref(const int i, const int j) const {assert(D->s==open);return &REF_SP(*this,i,j);}
    double read(const int i, const int j) const {assert(D->s==open);return REF_SP(*this,i,j);}
    void write(const int i, const int j, const double v) {assert(D->s==open);D=D->get_own_copy();REF_SP(*this,i,j)=v;}
    void add(const int i, const int j, const double v) {assert(D->s==open);D=D->get_own_copy();REF_SP(*this,i,j)+=v;}
    void sub(const int i, const int j, const double v) {assert(D->s==open);D=D->get_own_copy();REF_SP(*this,i,j)-=v;}
    void mul(const int i, const int j, const double v) {assert(D->s==open);D=D->get_own_copy();REF_SP(*this,i,j)*=v;}
    void div(const int i, const int j, const double v) {assert(D->s==open);D=D->get_own_copy();REF_SP(*this,i,j)/=v;}
public:
/*
 *  
 *  virtual MatrixC& assign (const int nr, const T* value)
 *  {
 *      D=D->get_own_copy();
 *      D->assign(nr,value);
 *      return *this;ve
 *  }
 *  
 *  virtual MatrixC& assign (const int nr, const int nc, const T** value)
 *  {
 *      D=D->get_own_copy();
 *      D->assign(nr,nc,value);
 *      return *this;
 *  }
 */
    //virtual operator double*() const {cout << "not yet implemented" <<endl;/*return &D->m[0];*/}
    //virtual operator double*() {cout << "not yet implemented" <<endl;/*return &D->m[0];*/}
    virtual bool operator == (const MatrixSP &B)    {return (*(D->v))==(*(B.D->v));}
    virtual bool operator <  (const MatrixSP &B) {return (*(D->v))<(*(B.D->v));}
protected:
};



ostream& operator<< (ostream& s, const MatrixSP& M)
{
    typedef map<int,double>::const_iterator ci;
    s.precision(OUTPUT_PREC) ;
    s.setf(ios::scientific,ios::floatfield) ;
    //  open 
    if (M.D->s==MatrixSP::open) 
    {
        s << M.nc() << endl;
        for (int i=0;i<M.nc();++i) 
        {
            s << M[i].size() << endl;
            ci j;
            for (j=M[i].begin();j!=M[i].end();++j) s << j->first << " ";
            s << endl ;
            for (j=M[i].begin();j!=M[i].end();++j) s << j->second << " ";
            s << endl ;
        }
    }
    //  closed
    else
    {
        valarray<int>& izh=*(M.D->izh);
        valarray<int>& ia=*(M.D->ia);
        valarray<double>& m=*(M.D->m);
        int *i,*izh_begin=&izh[0],*izh_end=&izh[izh.size()],*ia_begin=&ia[0],*ia_end=&ia[ia.size()];
        double* data=&m[0];
        s << izh.size() << endl;
        for (i=izh_begin;i!=izh_end;++i)    s << *i <<" " ;
        s << endl ;
        s << ia.size() << endl;
        for (i=ia_begin;i<ia_end;++i)
        {
            s << *i <<" "<< *data <<" " ;
            data++;
        }
        s << endl ;
    }
    s.setf(0,ios::floatfield) ;
    return s ;
}
istream& operator>> (istream& s, MatrixSP& M)
{
    // open 
    if (M.D->s==MatrixSP::open) 
    {
        for (int i=0;i<M.nc();++i) 
        {
            int nr;
            s >> nr;
            for (int j=0;j<nr;++j) 
            {
                int k;double value;
                s >> k >> value;
                M(k,i)=value;
            }
        }
    }
    // closed
    else 
    {
        M.D->izh->resize(M.D->ncol);
        valarray<int>& izh=*(M.D->izh);
        valarray<int>& ia=*(M.D->ia);
        valarray<double>& m=*(M.D->m);
        int *i,*izh_begin=&izh[0],*izh_end=&izh[izh.size()],*ia_begin=&ia[0],*ia_end=&ia[ia.size()];
        double* data=&m[0];
        for (i=izh_begin;i!=izh_end;++i) {cin >> *i ;}
        for (i=ia_begin;i!=ia_end;++i) {
            cin >> *i >> *data ;
            ++data;
        }
    }
    return s;
}
#ifdef USE_NAMESPACE_LINAL
}   // namespace LinAL
#endif
#endif //MATRIXSP_H

---