CMgramsparse_withoutdiag.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_CMGRAMSPARSE_WITHOUTDIAG_HXX
#define CONICBUNDLE_CMGRAMSPARSE_WITHOUTDIAG_HXX

#include "Coeffmat.hxx"
#include "CMsymdense.hxx"

namespace ConicBundle {


  class CMgramsparse_withoutdiag: public Coeffmat
  {
  private:
    CH_Matrix_Classes::Sparsemat A;  
    CH_Matrix_Classes::Sparsesym di; 
    bool positive; 
  public:
    CMgramsparse_withoutdiag(const CH_Matrix_Classes::Sparsemat& Ain,bool pos=true,CoeffmatInfo* cip=0)
    {A=Ain;di=sparseDiag(rowsip(A));positive=pos;CM_type=CM_gramsparsewd;infop=cip;}
    virtual ~CMgramsparse_withoutdiag(){}

    
    virtual Coeffmat* clone() const 
    {return new CMgramsparse_withoutdiag(A,positive,ConicBundle::clone(infop));}

    virtual CH_Matrix_Classes::Integer dim() const  { return A.rowdim(); }
    
    virtual CH_Matrix_Classes::Real operator()(CH_Matrix_Classes::Integer i,CH_Matrix_Classes::Integer j) const 
    {if (i==j) return 0.; return  (positive? 1.:-1.)*CH_Matrix_Classes::ip(A.row(i),A.row(j));}

    virtual void make_symmatrix(CH_Matrix_Classes::Symmatrix& S) const
    {if (positive) { 
        CH_Matrix_Classes::rankadd(A,S);
        for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
          CH_Matrix_Classes::Integer j=di.get_colindex()(i);
          S(j,j)-=di.get_colval()(i);
        }
      }
      else { 
        CH_Matrix_Classes::rankadd(A,S,-1.);
        for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
          CH_Matrix_Classes::Integer j=di.get_colindex()(i);
          S(j,j)+=di.get_colval()(i);
        }
      }
    }

    virtual CH_Matrix_Classes::Real norm(void) const 
    { CH_Matrix_Classes::Real sum=0.; CH_Matrix_Classes::Matrix vec; 
      for(CH_Matrix_Classes::Integer i=0;i<A.get_rowinfo().rowdim();i++){
        CH_Matrix_Classes::genmult(A,A.row(A.get_rowinfo()(i,0)),vec,1.,0.,0,1);
        sum+=CH_Matrix_Classes::ip(vec,vec)-di(A.get_rowinfo()(i,0),A.get_rowinfo()(i,0));
      }
      return std::sqrt(sum);
    }

    virtual Coeffmat* subspace(const CH_Matrix_Classes::Matrix& P) const
    { CH_Matrix_Classes::Matrix B;
      CH_Matrix_Classes::genmult(P,A,B,1.,0.,1);
      CH_Matrix_Classes::Symmatrix S;
      CH_Matrix_Classes::rankadd(B,S);
      CH_Matrix_Classes::genmult(P,di,B,1.,0.,1);
      S-=CH_Matrix_Classes::Symmatrix(B*P);
      if (!positive) S*=-1;
      return new CMsymdense(S,ConicBundle::clone(infop)); 
    }
    
    virtual void multiply(CH_Matrix_Classes::Real d)
    { di*=std::abs(d);
      if (d<0.) {A*=sqrt(-d); positive=!positive;} else {A*=sqrt(d);} 
      if (infop) infop->multiply(d); }

    virtual CH_Matrix_Classes::Real ip(const CH_Matrix_Classes::Symmatrix& S) const
    {CH_Matrix_Classes::Matrix B; if (positive) return CH_Matrix_Classes::ip(CH_Matrix_Classes::genmult(S,A,B),A)-CH_Matrix_Classes::ip(S,di);
     else return -CH_Matrix_Classes::ip(CH_Matrix_Classes::genmult(S,A,B),A)+CH_Matrix_Classes::ip(S,di); }
    
    virtual CH_Matrix_Classes::Real gramip(const CH_Matrix_Classes::Matrix& P) const
    {
      CH_Matrix_Classes::Real sum=0.;
      for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
        CH_Matrix_Classes::Matrix vec=P.row(di.get_colindex()(i));
        sum+=CH_Matrix_Classes::ip(vec,vec)*di.get_colval()(i);
      }
      CH_Matrix_Classes::Matrix B; CH_Matrix_Classes::genmult(P,A,B,1.,0.,1);
      if (positive) return CH_Matrix_Classes::ip(B,B)-sum; 
      else return -CH_Matrix_Classes::ip(B,B)+sum;
    }

    virtual CH_Matrix_Classes::Real gramip(const CH_Matrix_Classes::Matrix& P,CH_Matrix_Classes::Integer start_row,const CH_Matrix_Classes::Matrix* Lam=0) const
    {
      CH_Matrix_Classes::Real trval=0.;
      for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
        CH_Matrix_Classes::Matrix vec=P.row(di.get_colindex()(i)+start_row);
        trval-=CH_Matrix_Classes::rowip(P,start_row+di.get_colindex()(i),Lam)*di.get_colval()(i);
      }
      CH_Matrix_Classes::Matrix B(A.coldim(),P.coldim(),0.);
      CH_Matrix_Classes::Real *bp=B.get_store();
      const CH_Matrix_Classes::Integer brd=B.rowdim();
      const CH_Matrix_Classes::Integer pcd=P.coldim();
      const CH_Matrix_Classes::Integer prd=P.rowdim();
      const CH_Matrix_Classes::Real *pp=P.get_store()+start_row;
      const CH_Matrix_Classes::Indexmatrix &colinfo=A.get_colinfo();
      const CH_Matrix_Classes::Real *vp=(A.get_colval()).get_store();
      const CH_Matrix_Classes::Integer *ip=(A.get_colindex()).get_store();
      for(CH_Matrix_Classes::Integer j=0;j<colinfo.rowdim();j++){
        CH_Matrix_Classes::Integer indj=colinfo(j,0);
        for(CH_Matrix_Classes::Integer i=colinfo(j,1);--i>=0;){
          CH_Matrix_Classes::mat_xpeya(pcd,bp+indj,brd,pp+(*ip++),prd,(*vp++));
        }
      }
      if (Lam==0) {
        trval+=CH_Matrix_Classes::ip(B,B);
      }
      else {
        assert(Lam->dim()==P.coldim());
        const CH_Matrix_Classes::Real *bp=B.get_store();
        const CH_Matrix_Classes::Real *lp=Lam->get_store();
        for (CH_Matrix_Classes::Integer i=0;i<B.coldim();i++,bp+=B.rowdim())
          trval+=(*lp++)*CH_Matrix_Classes::mat_ip(B.rowdim(),bp);
      }
      return (positive?trval:-trval);
    }

    virtual void addmeto(CH_Matrix_Classes::Symmatrix& S,CH_Matrix_Classes::Real d=1.) const
    {if (!positive) d*=-1;
     CH_Matrix_Classes::rankadd(A,S,d,1.); S-=di*d;}

    virtual void addprodto(CH_Matrix_Classes::Matrix& B,const CH_Matrix_Classes::Matrix&C ,CH_Matrix_Classes::Real d=1.) const
    {CH_Matrix_Classes::Matrix D;if (!positive) d*=-1;
     CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(A,C,D,1.,0.,1),B,d,1.);
     CH_Matrix_Classes::genmult(di,C,B,-d,1.,0);
     }

    virtual void addprodto(CH_Matrix_Classes::Matrix& B,const CH_Matrix_Classes::Sparsemat&C ,CH_Matrix_Classes::Real d=1.) const
    { if (!positive) d*=-1;
      CH_Matrix_Classes::Matrix D;
     CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(A,C,D,1.,0.,1),B,d,1.);
     CH_Matrix_Classes::genmult(di,C,B,-d,1.,0);}

    virtual void left_right_prod(const CH_Matrix_Classes::Matrix& P,const CH_Matrix_Classes::Matrix& Q,CH_Matrix_Classes::Matrix& R) const
    {
      R.init(P.coldim(),Q.coldim(),0.);
      for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
        genmult(P.row(di.get_colindex()(i)),Q.row(di.get_colindex()(i)),R,di.get_colval()(i),1.,1);
      }
      CH_Matrix_Classes::Matrix tmp1; CH_Matrix_Classes::genmult(P,A,tmp1,1.,0.,1,0);
      CH_Matrix_Classes::Matrix tmp2; CH_Matrix_Classes::genmult(A,Q,tmp2,1.,0.,1,0);
      if (positive) CH_Matrix_Classes::genmult(tmp1,tmp2,R,1.,-1.,0,0);
      else CH_Matrix_Classes::genmult(tmp1,tmp2,R,-1.,1.,0,0);
    }

    virtual CH_Matrix_Classes::Integer prodvec_flops() const 
    { return 4*A.nonzeros()+2*di.get_colindex().dim(); }

    virtual int dense() const
    {return 0;}
    
    virtual int sparse() const
    { return 0;}
    
   virtual int sparse(CH_Matrix_Classes::Indexmatrix& /* I */,
                       CH_Matrix_Classes::Indexmatrix& /* J */,
                       CH_Matrix_Classes::Matrix& /* val */,
                       CH_Matrix_Classes::Real /* d=1. */)const
    {return 0;}
    
    virtual int support_in(const CH_Matrix_Classes::Sparsesym& /* S */) const
    {return 0;}

    virtual CH_Matrix_Classes::Real ip(const CH_Matrix_Classes::Sparsesym& S) const
    {if (positive) return CH_Matrix_Classes::ip(A,S*A)-CH_Matrix_Classes::ip(di,S); else return -CH_Matrix_Classes::ip(A,S*A)+CH_Matrix_Classes::ip(di,S);}
    
    virtual void project(CH_Matrix_Classes::Symmatrix& S,const CH_Matrix_Classes::Matrix& P) const
    {CH_Matrix_Classes::Matrix B; CH_Matrix_Classes::genmult(P,A,B,1.,0.,1);
     if (positive) CH_Matrix_Classes::rankadd(B,S); else CH_Matrix_Classes::rankadd(B,S,-1.);
     for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
       rankadd(P.row(di.get_colindex()(i)),S,(positive?-1.:1.)*di.get_colval()(i),1.,1);
     } 
    }

    virtual void add_projection(CH_Matrix_Classes::Symmatrix& S,const CH_Matrix_Classes::Matrix& P,CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Integer start_row=0) const
    {
      CH_Matrix_Classes::Matrix B(A.coldim(),P.coldim(),0.);
      CH_Matrix_Classes::Real *bp=B.get_store();
      const CH_Matrix_Classes::Integer brd=B.rowdim();
      const CH_Matrix_Classes::Integer pcd=P.coldim();
      const CH_Matrix_Classes::Integer prd=P.rowdim();
      const CH_Matrix_Classes::Real *pp=P.get_store()+start_row;
      const CH_Matrix_Classes::Indexmatrix &colinfo=A.get_colinfo();
      const CH_Matrix_Classes::Real *vp=(A.get_colval()).get_store();
      const CH_Matrix_Classes::Integer *ip=(A.get_colindex()).get_store();
      for(CH_Matrix_Classes::Integer j=0;j<colinfo.rowdim();j++){
        CH_Matrix_Classes::Integer indj=colinfo(j,0);
        for(CH_Matrix_Classes::Integer i=colinfo(j,1);--i>=0;){
          CH_Matrix_Classes::mat_xpeya(pcd,bp+indj,brd,pp+(*ip++),prd,(*vp++));
        }
      }
      if (positive) CH_Matrix_Classes::rankadd(B,S,alpha,1.,1); else CH_Matrix_Classes::rankadd(B,S,-alpha,1.,1);      
      for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
        rankadd(P.row(di.get_colindex()(i)+start_row),S,(positive?-alpha:alpha)*di.get_colval()(i),1.,1);
     } 
    }

    virtual const CH_Matrix_Classes::Matrix& postgenmult(const CH_Matrix_Classes::Matrix& B,CH_Matrix_Classes::Matrix& C,
                             CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int btrans=0) const
    { 
      CH_Matrix_Classes::Matrix D; 
      CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(A,B,D,1.,0.,1,btrans),C,(positive?1.:-1.)*alpha,beta);
      return CH_Matrix_Classes::genmult(di,B,C,(positive?-1.:1.)*alpha,1.,btrans);
    }

    virtual const CH_Matrix_Classes::Matrix& pregenmult(const CH_Matrix_Classes::Matrix& B,CH_Matrix_Classes::Matrix& C,
                             CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int btrans=0) const
    { 
      CH_Matrix_Classes::Matrix D; 
      return CH_Matrix_Classes::genmult(CH_Matrix_Classes::genmult(B,A,D,1.,0.,btrans),A,C,(positive?1.:-1.)*alpha,beta,0,1);
      return CH_Matrix_Classes::genmult(B,di,C,(positive?-1.:1.)*alpha,1.,btrans);
    }

    virtual int equal(const Coeffmat* p,double tol=1e-6) const
    {
      const CMgramsparse_withoutdiag *pp=dynamic_cast<const CMgramsparse_withoutdiag *>(p);
      if (pp==0) 
        return 0;
      return (positive==pp->positive) && CH_Matrix_Classes::equal(A,pp->A,tol);
    }

    virtual std::ostream& display(std::ostream& o) const 
    {o<<"CMgramsparse_withoutdiag\n"; A.display(o); return o;}
    
    virtual std::ostream& out(std::ostream& o) const
    {return o<<"GRAM_SPARSE_WITHOUTDIAG\n"<<positive<<"\n"<<A;}

    virtual std::istream& in(std::istream& i)
    { i>>positive>>A; di=sparseDiag(rowsip(A)); return i;}

    CMgramsparse_withoutdiag(std::istream& is,CoeffmatInfo* cip=0)
    {CM_type=CM_gramsparse;infop=cip;in(is);}

    //--- specific routines
    const CH_Matrix_Classes::Sparsemat& get_A() const {return A;}
    bool get_positive() const {return positive;}
};


}
#endif