symmat.hxx

Go to the documentation of this file.

#ifndef CH_MATRIX_CLASSES__SYMMAT_HXX
#define CH_MATRIX_CLASSES__SYMMAT_HXX

#ifndef CH_MATRIX_CLASSES__MATRIX_HXX
#include "matrix.hxx"
#endif

namespace CH_Matrix_Classes {


// **************************************************************************
//                            class definition
// **************************************************************************


class Symmatrix: protected Memarrayuser
{
    friend class Matrix; 
    friend class Sparsesym;
    friend class Sparsemat;
    
private:

    static const Mtype mtype;   
    Integer mem_dim;   
    Integer nr;        
    Real *m;     

    bool is_init;   

    inline void init_to_zero();        
    
    Integer tred2(Integer nm, Integer n, Real *a, Real *d, Real *e, Real *z) const;
    Integer imtql2(Integer nm, Integer n, Real *d, Real *e, Real *z) const;

public:

  //----------------------------------------------------
  //----  constructors, destructor, and initialization
  //----------------------------------------------------



  inline Symmatrix(); 
  inline Symmatrix(const Symmatrix& A,double d=1.);
  inline Symmatrix(Integer nr);
  inline Symmatrix(Integer nr,Real d);   
  inline Symmatrix(Integer nr,const Real* dp);
  inline ~Symmatrix();                   
  
#if (CONICBUNDLE_DEBUG>=1)
  void set_init(bool i){is_init=i;}
 int get_init() const {return is_init;}
#else 
  void set_init(bool /* i */){}
  bool get_init() const {return true;}
#endif

  inline Symmatrix& init(const Symmatrix& A,double d=1.);
  inline Symmatrix& init(const Matrix& A,double d=1.);
  inline Symmatrix& init(const Indexmatrix& A,double d=1.);
  inline Symmatrix& init(const Sparsesym& A,Real d=1.);
  inline Symmatrix& init(Integer nr,Real d);
  inline Symmatrix& init(Integer nr,const Real* dp);
    
  void newsize(Integer n);      //resize matrix without Initialization




  inline Symmatrix(const Matrix&,double d=1.);
  inline Symmatrix(const Indexmatrix&,double d=1.);
  inline Symmatrix(const Sparsesym& A,Real d=1.); 
  

  //----------------------------------------------------
  //----  size and type information
  //----------------------------------------------------


  void dim(Integer& _nr, Integer& _nc) const {_nr=nr; _nc=nr;}

  Integer dim() const {return nr;}

  Integer rowdim() const {return nr;}

  Integer coldim() const {return nr;}
  
  Mtype get_mtype() const {return mtype;}

    

  //--------------------------------
  //----  Indexing and Submatrices
  //--------------------------------


  inline Real& operator()(Integer i,Integer j);

  inline Real& operator()(Integer i);

  inline Real operator()(Integer i,Integer j) const;

  inline Real operator()(Integer i) const;
  
  Matrix col(Integer i) const;     
  Matrix row(Integer i) const;       
  Matrix cols(const Indexmatrix &vec) const;  //returns cols as indexed by vec
  Matrix rows(const Indexmatrix &vec) const;  //returns rows as indexed by vec

  Symmatrix& swapij(Integer i,Integer j);
 
  Symmatrix& pivot_permute(const Indexmatrix& piv,bool inverse=false);

  Symmatrix& principal_submatrix(const Indexmatrix& ind,Symmatrix& S) const; 

  inline Symmatrix principal_submatrix(const Indexmatrix& ind) const;


  Symmatrix& delete_principal_submatrix(const Indexmatrix& ind,bool sorted_increasingly=false);

  Symmatrix& enlarge_below(Integer addn);

  Symmatrix& enlarge_below(Integer addn,Real d);


 
   Real* get_store() {return m;}   //use cautiously, do not use delete!
  const Real* get_store() const {return m;}   //use cautiously




  friend Matrix diag(const Symmatrix& A);    //=(A(1,1),A(2,2),...)^t
  friend Symmatrix Diag(const Matrix& A);    //vec(A) on the diagonal

  friend inline void swap(Symmatrix& A,Symmatrix& B);


  //------------------------------
  //----  BLAS-like Routines
  //------------------------------


  Symmatrix& xeya(const Symmatrix& A,Real d=1.);   //*this=d*A
  Symmatrix& xpeya(const Symmatrix& A,Real d=1.);  //*this+=d*A;

    


  friend Symmatrix& rankadd(const Matrix& A,Symmatrix& C,
                            Real alpha,Real beta,int trans);
  friend Symmatrix& scaledrankadd(const Matrix& A,const Matrix& D,Symmatrix& C,
                            Real alpha,Real beta,int trans);
  friend Symmatrix& rank2add(const Matrix& A, const Matrix& B, Symmatrix& C,
                             Real alpha,Real beta,int trans);
   
  
  friend Symmatrix& xbpeya(Symmatrix& x,const Symmatrix& y,Real alpha,Real beta);
  

  friend Symmatrix& xeyapzb(Symmatrix& x,const Symmatrix& y,const Symmatrix& z,Real alpha,Real beta);
  
  friend Matrix& genmult(const Symmatrix& A,const Matrix& B,Matrix& C,
                         Real alpha, Real beta, int btrans);

  friend Matrix& genmult(const Matrix& A,const Symmatrix& B,Matrix& C,
                         Real alpha, Real beta, int atrans);


  //------------------------------
  //----  usual operators
  //------------------------------


  inline Symmatrix& operator=(const Symmatrix &A); 
  inline Symmatrix& operator+=(const Symmatrix &A);
  inline Symmatrix& operator-=(const Symmatrix &A);
  inline Symmatrix& operator%=(const Symmatrix &A);  //Hadamard product
  inline Symmatrix operator-() const; 
  
  inline Symmatrix& operator*=(Real d);
  inline Symmatrix& operator/=(Real d);
  inline Symmatrix& operator+=(Real d);
  inline Symmatrix& operator-=(Real d);

  Symmatrix& transpose(){return *this;} 



  friend inline Matrix operator*(const Symmatrix &A,const Symmatrix &B);
  friend inline Symmatrix operator%(const Symmatrix &A,const Symmatrix &B);
  friend inline Symmatrix operator+(const Symmatrix &A,const Symmatrix &B);
  friend inline Symmatrix operator-(const Symmatrix &A,const Symmatrix &B);
  friend inline Matrix operator*(const Symmatrix &A,const Matrix &B); 
  friend inline Matrix operator*(const Matrix &A,const Symmatrix &B);
  friend inline Matrix operator+(const Symmatrix &A,const Matrix &B);
  friend inline Matrix operator+(const Matrix &A,const Symmatrix &B);
  friend inline Matrix operator-(const Symmatrix &A,const Matrix &B);
  friend inline Matrix operator-(const Matrix &A,const Symmatrix &B);

  friend inline Symmatrix operator*(const Symmatrix &A,Real d);
  friend inline Symmatrix operator*(Real d,const Symmatrix &A);
  friend inline Symmatrix operator/(const Symmatrix &A,Real d);
  friend inline Symmatrix operator+(const Symmatrix &A,Real d);
  friend inline Symmatrix operator+(Real d,const Symmatrix &A);
  friend inline Symmatrix operator-(const Symmatrix &A,Real d);
  friend inline Symmatrix operator-(Real d,const Symmatrix &A);

  friend inline Symmatrix transpose(const Symmatrix& A);


  //------------------------------------------
  //----  Connections to other Matrix Classes
  //------------------------------------------
  

  Symmatrix& xeya(const Matrix& A,Real d=1.);   //*this=d*A
  Symmatrix& xpeya(const Matrix& A,Real d=1.);  //*this+=d*A;
  Symmatrix& xeya(const Indexmatrix& A,Real d=1.);   //*this=d*A
  Symmatrix& xpeya(const Indexmatrix& A,Real d=1.);  //*this+=d*A;
  Symmatrix& xeya(const Sparsesym& A,Real d=1.);   //*this=d*A
  Symmatrix& xpeya(const Sparsesym& A,Real d=1.);  //*this+=d*A;

  Symmatrix& xetriu_yza(const Matrix& A,const Matrix& B,Real d=1.);
  Symmatrix& xpetriu_yza(const Matrix& A,const Matrix& B,Real d=1.);
  Symmatrix& xetriu_yza(const Sparsemat& A,const Matrix& B,Real d=1.);
  Symmatrix& xpetriu_yza(const Sparsemat& A,const Matrix& B,Real d=1.);

  inline Symmatrix& operator=(const Sparsesym& A);
  inline Symmatrix& operator+=(const Sparsesym& A);
  inline Symmatrix& operator-=(const Sparsesym& A);



  friend Matrix& genmult(const Symmatrix& A,const Sparsemat& B,Matrix& C,
                         Real alpha,Real beta, int btrans);

  friend Matrix& genmult(const Sparsemat& A,const Symmatrix& B,Matrix& C,
                           Real alpha,Real beta, int atrans);

  friend Symmatrix& rankadd(const Sparsemat& A,Symmatrix& C,
                            Real alpha,Real beta,int trans);
    
  friend Symmatrix& scaledrankadd(const Sparsemat& A,const Matrix& D,Symmatrix& C,
                            Real alpha,Real beta,int trans);

  friend Symmatrix& rank2add(const Sparsemat& A,const Matrix& B,Symmatrix& C,
                             Real alpha,Real beta,int trans);

   
  //------------------------------
  //----  Elementwise Operations
  //------------------------------



  friend Symmatrix abs(const Symmatrix& A);                 


  //----------------------------
  //----  Numerical Methods
  //----------------------------


  Symmatrix& shift_diag(Real s);             //add s to all diagonal elements

  //----- LDL Factorization (for positive definite matrices, no pivots so far!)
  int LDLfactor(Real tol=1e-10);     
  int LDLsolve(Matrix& x) const;      //call LDLfactor before using LDLsolve!
  int LDLinverse(Symmatrix& S) const; //call LDLfactor before using LDLinverse!
    
  //----- Cholesky Factorization with pivoting
  int Chol_factor(Real tol=1e-10); //stores fact. in *this
  int Chol_solve(Matrix& x) const; //call _factor before
  int Chol_inverse(Symmatrix& S) const; // --- " ---
  int Chol_Lsolve(Matrix& rhs) const; // --- " ---, solve only Ly=x
  int Chol_Ltsolve(Matrix& rhs) const; // --- " ---, solve only L^Ty=x
  int Chol_scaleLi(Symmatrix& S) const;
  int Chol_scaleLt(Symmatrix& S) const;
  int Chol_Lmult(Matrix& rhs) const; 
  int Chol_Ltmult(Matrix& rhs) const; 
  int Chol_factor(Indexmatrix &piv, Real tol=1e-10); //stores fact. in *this
  int Chol_solve(Matrix& x,const Indexmatrix& piv) const; //call _factor before
  int Chol_inverse(Symmatrix& S,const Indexmatrix & piv) const; // --- " ---

  //----- Aasen factorization with pivoting
  int Aasen_factor(Indexmatrix &piv); 
  int Aasen_Lsolve(Matrix& x) const;
  int Aasen_Ltsolve(Matrix& x) const;
  int Aasen_tridiagsolve(Matrix& x) const;
  int Aasen_solve(Matrix& x, const Indexmatrix& piv) const; //call Aasen_factor before!

  //----- Eigenvalue decomposition
  Integer eig(Matrix& P,Matrix& d,bool sort_non_decreasingly=true) const;
  //if return value ==0 then P contains eigenvectors as columns
  //and d is a column vector containing the eigenvalues.



  friend Real trace(const Symmatrix& A);                  //=sum(diag(A))
  friend Real ip(const Symmatrix& A, const Symmatrix& B); //=trace(B^t*A)
  friend Real ip(const Matrix& A, const Symmatrix& B);    //=trace(B^t*A)
  friend Real ip(const Symmatrix& A, const Matrix& B);    //=trace(B^t*A)
  friend inline Real norm2(const Symmatrix& A); //{return sqrt(ip(A,A));}

  friend Matrix sumrows(const Symmatrix& A); //=(1 1 1 ... 1)*A
  friend Matrix sumcols(const Symmatrix& A); //=A*(1 1 ... 1)^t
  friend Real sum(const Symmatrix& A);       //=(1 1 ... 1)*A*(1 1 ... 1)^t
  
  //---- svec and skron

  friend void svec(const Symmatrix& A,Matrix& v,Real a,bool add,Integer startindex_vec,Integer startindex_A,Integer blockdim); 
  friend inline Matrix svec(const Symmatrix& A);
  friend void sveci(const Matrix& v,Symmatrix& A,Real a,bool add,Integer startindex_vec,Integer startindex_A,Integer blockdim);

  // initialize from an svec stored in a real array (or matrix)
  void init_svec(Integer nr,const Real *dp,Integer incr=1,Real d=1.);
  // store in the form of an svec in the real array (or matrix)
  void store_svec(Real *dp,Integer incr=1,Real d=1.) const;
  
  friend inline Symmatrix skron(const Symmatrix& A,const Symmatrix& B,Real alpha,bool add,Integer startindex_S);
  friend void skron(const Symmatrix& A,const Symmatrix& B,Symmatrix& S,Real a,bool add,Integer startindex_S);
  friend void symscale(const Symmatrix &A,const Matrix& B,Symmatrix& S,Real a,Real b,int btrans);


  //---------------------------------------------
  //----  Comparisons / Max / Min / sort / find
  //---------------------------------------------


  friend Matrix minrows(const Symmatrix& A); //minimum of each column (over rows)
  friend Matrix mincols(const Symmatrix& A); //minimum of each row (over colums)
  friend Real min(const Symmatrix& A);       //minimum of matrix
  friend Matrix maxrows(const Symmatrix& A); //similar
  friend Matrix maxcols(const Symmatrix& A); //similar
 friend Real max(const Symmatrix& A);       //similar



  //--------------------------------
  //----  Input / Output
  //--------------------------------



  void display(std::ostream& out, 
                 int precision=0,   
                 int width=0,       
                 int screenwidth=0  
               ) const;

    void mfile_output(std::ostream& out, 
                 int precision=16,   
                 int width=0       
                 ) const;




  friend std::ostream& operator<<(std::ostream& o,const Symmatrix &A);

  friend std::istream& operator>>(std::istream& i,Symmatrix &A);
  


};



// **************************************************************************
//                make non inline friends available outside
// **************************************************************************
  
Matrix diag(const Symmatrix& A);    //=(A(1,1),A(2,2),...)^t

Symmatrix Diag(const Matrix& A);    //vec(A) on the diagonal

Symmatrix& rankadd(const Matrix& A,Symmatrix& C,
                            Real alpha=1.,Real beta=0.,int trans=0);

Symmatrix& scaledrankadd(const Matrix& A,const Matrix& D,Symmatrix& C,
                            Real alpha=1.,Real beta=0.,int trans=0);

Symmatrix& rank2add(const Matrix& A, const Matrix& B, Symmatrix& C,
                             Real alpha=1.,Real beta=0.,int trans=0);

Matrix& genmult(const Symmatrix& A,const Matrix& B,Matrix& C,
                         Real alpha=1., Real beta=0., int btrans=0);

Matrix& genmult(const Matrix& A,const Symmatrix& B,Matrix& C,
                         Real alpha=1., Real beta=0., int atrans=0);

Matrix& genmult(const Symmatrix& A,const Sparsemat& B,Matrix& C,
                Real alpha=1.,Real beta=0., int btrans=0);

Matrix& genmult(const Sparsemat& A,const Symmatrix& B,Matrix& C,
                Real alpha=1.,Real beta=0., int atrans=0);

Symmatrix& rankadd(const Sparsemat& A,Symmatrix& C,
                   Real alpha=1.,Real beta=0.,int trans=0);

Symmatrix& scaledrankadd(const Sparsemat& A,const Matrix& D,Symmatrix& C,
                            Real alpha=1.,Real beta=0.,int trans=0);

Symmatrix& rank2add(const Sparsemat& A,const Matrix& B,Symmatrix& C,
                    Real alpha=1.,Real beta=0.,int trans=0);

Symmatrix abs(const Symmatrix& A);

Real trace(const Symmatrix& A);                  //=sum(diag(A))

Real ip(const Symmatrix& A, const Symmatrix& B); //=trace(B^t*A)

Real ip(const Matrix& A, const Symmatrix& B);    //=trace(B^t*A)

Real ip(const Symmatrix& A, const Matrix& B);    //=trace(B^t*A)

Matrix sumrows(const Symmatrix& A); //=(1 1 1 ... 1)*A

Matrix sumcols(const Symmatrix& A); //=A*(1 1 ... 1)^t

Real sum(const Symmatrix& A);       //=(1 1 ... 1)*A*(1 1 ... 1)^t


  void svec(const Symmatrix& A,Matrix& sv,Real a=1.,bool add=false,Integer startindex_vec=-1,Integer startindex_A=0,Integer blockdim=-1); //stores svec in sv

void sveci(const Matrix& sv,Symmatrix& A,Real a=1.,bool add=false,Integer startindex_vec=0,Integer startindex_A=-1,Integer blockdim=-1);

void skron(const Symmatrix& A,const Symmatrix& B,Symmatrix& S,Real alpha=1.,bool add=false,Integer startindex_S=-1);

void symscale(const Symmatrix &A,const Matrix& B,Symmatrix& S,Real alpha=1.,Real beta=0.,int btrans=0);

Matrix minrows(const Symmatrix& A); //minimum of each column (over rows)

Matrix mincols(const Symmatrix& A); //minimum of each row (over colums)

Real min(const Symmatrix& A);       //minimum of matrix

Matrix maxrows(const Symmatrix& A); //similar

Matrix maxcols(const Symmatrix& A); //similar

Real max(const Symmatrix& A);       //similar

  std::ostream& operator<<(std::ostream& o,const Symmatrix &A);

std::istream& operator>>(std::istream& i,Symmatrix &A);


// **************************************************************************
//                    implemenation of inline functions
// **************************************************************************

inline void Symmatrix::init_to_zero()
{
 nr=0;mem_dim=0;m=0;
#if (CONICBUNDLE_DEBUG>=1)
 is_init=true;
#endif
} 

inline Symmatrix& Symmatrix::init(const Symmatrix& A,double d) 
{ return xeya(A,d); }

inline Symmatrix& Symmatrix::init(const Matrix& A,double d) 
{ return xeya(A,d); }

inline Symmatrix& Symmatrix::init(const Indexmatrix& A,double d) 
{ return xeya(A,d); }

inline Symmatrix& Symmatrix::init(Integer inr,Real d)
{
 newsize(inr);
 mat_xea((nr*(nr+1))/2,m,d);
 chk_set_init(*this,1);
 return *this;
}

inline Symmatrix& Symmatrix::init(Integer inr,const Real* pd)
{
 newsize(inr);
 mat_xey((nr*(nr+1))/2,m,pd);
 chk_set_init(*this,1);
 return *this;
}

inline Symmatrix::Symmatrix() 
{ init_to_zero(); }
    
inline Symmatrix::Symmatrix(const Symmatrix& A,double d):Memarrayuser()
{init_to_zero(); xeya(A,d);}     

inline Symmatrix::Symmatrix(const Matrix &M,double d)
{ init_to_zero(); xeya(M,d); }

inline Symmatrix::Symmatrix(const Indexmatrix &M,double d)
{ init_to_zero(); xeya(M,d); }

inline Symmatrix::Symmatrix(Integer i)
{ init_to_zero(); newsize(i); }

inline Symmatrix::Symmatrix(Integer i,Real d)
{ init_to_zero(); init(i,d); }

inline Symmatrix::Symmatrix(Integer i,const Real *pd)
{ init_to_zero(); init(i,pd); }

inline Symmatrix::~Symmatrix()
{memarray->free(m);}

inline Real& Symmatrix::operator()(Integer i,Integer j) 
{
 chk_range(i,j,nr,nr);
 //if (j>=i) return m[((i*((nr<<1)-(i+1)))>>1)+j];
 //return m[((j*((nr<<1)-(j+1)))>>1)+i];
 //now unsigned is used to avoid "cc1plus: warning: assuming signed overflow does not occur when assuming that (X + c) < X is always false [-Wstrict-overflow]"
 if (unsigned(j)>=unsigned(i)) return m[((unsigned(i)*((unsigned(nr)<<1)-(unsigned(i)+1)))>>1)+unsigned(j)];
 return m[((unsigned(j)*((unsigned(nr)<<1)-(unsigned(j)+1)))>>1)+unsigned(i)];
}

inline Real Symmatrix::operator()(Integer i,Integer j) const
{
 chk_range(i,j,nr,nr);
 if (unsigned(j)>=unsigned(i)) return m[((unsigned(i)*((unsigned(nr)<<1)-(unsigned(i)+1)))>>1)+unsigned(j)];
 return m[((unsigned(j)*((unsigned(nr)<<1)-(unsigned(j)+1)))>>1)+unsigned(i)];
}

inline Real& Symmatrix::operator()(Integer i)
{
 chk_range(i,0,nr*nr,1);
 // Integer j=i/nr;
 // i=i%nr;
 // if (j>=i) return m[((i*((nr<<1)-(i+1)))>>1)+j];
 // return m[((j*((nr<<1)-(j+1)))>>1)+i];
 const unsigned long uj=unsigned(i/nr);
 const unsigned long ui=unsigned(i%nr);
 if (uj>=ui) return m[((ui*((unsigned(nr)<<1)-(ui+1)))>>1)+uj];
 return m[((uj*((unsigned(nr)<<1)-(uj+1)))>>1)+ui];
}

inline Real Symmatrix::operator()(Integer i) const
{
 chk_range(i,0,nr*nr,1);
 const unsigned long uj=unsigned(i/nr);
 const unsigned long ui=unsigned(i%nr);
 if (uj>=ui) return m[((ui*((unsigned(nr)<<1)-(ui+1)))>>1)+uj];
 return m[((uj*((unsigned(nr)<<1)-(uj+1)))>>1)+ui];
}

inline Symmatrix Symmatrix::principal_submatrix(const Indexmatrix& ind) const
{
  Symmatrix S;
  principal_submatrix(ind,S);
  return S;
}

inline void swap(Symmatrix &A,Symmatrix &B)
{ 
 Real *hm=A.m;A.m=B.m;B.m=hm;
 Integer hi=A.nr;A.nr=B.nr;B.nr=hi;
 hi=A.mem_dim;A.mem_dim=B.mem_dim;B.mem_dim=hi;
#if (CONICBUNDLE_DEBUG>=1)
 bool hb=A.is_init;A.is_init=B.is_init;B.is_init=hb;
#endif
}

inline Symmatrix& xbpeya(Symmatrix& x,const Symmatrix& y,Real alpha=1.,Real beta=0.)
  //returns x= alpha*y+beta*x, where y may be transposed (ytrans=1)
  //if beta==0. then x is initialized to the correct size
{
  if (beta==0.){
    chk_init(y);
    x.init(y,alpha);
  }
  else {
    chk_add(x,y);
    mat_xbpeya((x.rowdim()*(x.rowdim()+1))/2,
               x.get_store(),y.get_store(),alpha,beta);
  }
  return x;
}

inline Symmatrix& xeyapzb(Symmatrix& x,const Symmatrix& y,const Symmatrix& z,Real alpha=1.,Real beta=1.)
  //returns x= alpha*y+beta*z
  //x is initialized to the correct size
{
  chk_add(y,z);
  x.newsize(y.rowdim()); chk_set_init(x,1);
  mat_xeyapzb((x.rowdim()*(x.rowdim()+1))/2,
              x.get_store(),y.get_store(),z.get_store(),alpha,beta);
  return x;
}
  
inline Symmatrix& Symmatrix::operator=(const Symmatrix &A)   
    {return xeya(A);}
inline Symmatrix& Symmatrix::operator+=(const Symmatrix &A) 
    {return xpeya(A);}
inline Symmatrix& Symmatrix::operator-=(const Symmatrix &A) 
    {return xpeya(A,-1.);}
inline Symmatrix& Symmatrix::operator%=(const Symmatrix &A)
    { chk_add(*this,A); mat_xhadey((nr*(nr+1))/2,m,A.m); return *this; } 
inline Symmatrix  Symmatrix::operator-() const    
    {return Symmatrix(*this,-1.);}

inline Symmatrix& Symmatrix::operator*=(Real d)
    {chk_init(*this); mat_xmultea(nr*(nr+1)/2,m,d); return *this;}
inline Symmatrix& Symmatrix::operator/=(Real d)
    {chk_init(*this); mat_xmultea(nr*(nr+1)/2,m,1./d); return *this;}
inline Symmatrix& Symmatrix::operator+=(Real d)
    {chk_init(*this); mat_xpea(nr*(nr+1)/2,m,d); return *this;}
inline Symmatrix& Symmatrix::operator-=(Real d)
    {chk_init(*this); mat_xpea(nr*(nr+1)/2,m,-d); return *this;}

inline Matrix operator*(const Symmatrix &A,const Symmatrix &B)
    {Matrix C;return genmult(Matrix(A),B,C);}

inline Symmatrix operator%(const Symmatrix &A,const Symmatrix &B) 
    {Symmatrix C(A);C%=B;return C;}

inline Symmatrix operator+(const Symmatrix &A,const Symmatrix &B) 
    {Symmatrix C(A);return C.xpeya(B);}

inline Symmatrix operator-(const Symmatrix &A,const Symmatrix &B)
    {Symmatrix C(A);return C.xpeya(B,-1.);}

inline Matrix operator*(const Symmatrix &A,const Matrix &B) 
    {Matrix C; return genmult(A,B,C);}

inline Matrix operator*(const Matrix &A,const Symmatrix &B)
    {Matrix C; return genmult(A,B,C);}

inline Matrix operator+(const Symmatrix &A,const Matrix &B)
    {Matrix C(A);return C.xpeya(B);}

inline Matrix operator+(const Matrix &A,const Symmatrix &B)
    {Matrix C(A);return C.xpeya(B);}

inline Matrix operator-(const Symmatrix &A,const Matrix &B)
    {Matrix C(A);return C.xpeya(B,-1);}

inline Matrix operator-(const Matrix &A,const Symmatrix &B)
    {Matrix C(A);return C.xpeya(B,-1);}

inline Symmatrix operator*(const Symmatrix& A,Real d)
    { return Symmatrix(A,d);}

inline Symmatrix operator*(Real d,const Symmatrix &A)
    { return Symmatrix(A,d);}

inline Symmatrix operator/(const Symmatrix& A,Real d)
    { return Symmatrix(A,1./d);}

inline Symmatrix operator+(const Symmatrix& A,Real d)        
    { Symmatrix B(A); return B+=d;}

inline Symmatrix operator+(Real d,const Symmatrix &A)
    {Symmatrix B(A); return B+=d;}

inline Symmatrix operator-(const Symmatrix& A,Real d)        
    { Symmatrix B(A); return B-=d;}

inline Symmatrix operator-(Real d,const Symmatrix &A)    
    { Symmatrix B(A,-1.);return B+=d;}


inline Matrix svec(const Symmatrix& A) //=(A(11),sqrt(2)A(12),...A(22)..)^t
    {Matrix sv;svec(A,sv);return sv;}

inline Symmatrix skron(const Symmatrix& A,const Symmatrix& B,Real alpha=1.,bool add=false,Integer startindex_S=-1)
{Symmatrix S;skron(A,B,S,alpha,add,startindex_S);return S;}

inline Real norm2(const Symmatrix& A)
    {return ::sqrt(ip(A,A));}

inline Symmatrix transpose(const Symmatrix& A)
    {return Symmatrix(A);}


inline Matrix::Matrix(const Symmatrix &A,Real d)

    {init_to_zero(); xeya(A,d);}

    inline Matrix& Matrix::init(const Symmatrix &A,Real d)

        { return xeya(A,d); }

    inline Matrix& Matrix::operator=(const Symmatrix &A)

        { return xeya(A); }

    inline Matrix& Matrix::operator*=(const Symmatrix &A)

        { Matrix C; return *this=genmult(*this,A,C);}

    inline Matrix& Matrix::operator+=(const Symmatrix &A)

        { return xpeya(A); }

    inline Matrix& Matrix::operator-=(const Symmatrix &A)

        { return xpeya(A,-1.); }


}

#ifndef CH_MATRIX_CLASSES__SPARSMAT_HXX
#include "sparsmat.hxx"
#endif

#endif