CMlowrankdd.hxx

Go to the documentation of this file.

#ifndef CONICBUNDLE_CMLOWRANKDD_HXX
#define CONICBUNDLE_CMLOWRANKDD_HXX

#include "Coeffmat.hxx"

namespace ConicBundle {


  class CMlowrankdd: public Coeffmat
  {
  private:
    CH_Matrix_Classes::Matrix A; 
    CH_Matrix_Classes::Matrix B; 
  public:
    CMlowrankdd(const CH_Matrix_Classes::Matrix& Ain,const CH_Matrix_Classes::Matrix& Bin,CoeffmatInfo* cip=0)
    {A=Ain;B=Bin;CM_type=CM_lowrankdd;infop=cip;}
    virtual ~CMlowrankdd(){}

    
    virtual Coeffmat* clone() const 
    {return new CMlowrankdd(A,B,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 
    {
      return 
        CH_Matrix_Classes::mat_ip(A.coldim(),A.get_store()+i,A.rowdim(),B.get_store()+j,B.rowdim())
       +CH_Matrix_Classes::mat_ip(A.coldim(),B.get_store()+i,B.rowdim(),A.get_store()+j,A.rowdim());
    }

    virtual void make_symmatrix(CH_Matrix_Classes::Symmatrix& S) const
    { CH_Matrix_Classes::rank2add(A,B,S,2.);}

    virtual CH_Matrix_Classes::Real norm(void) const
    { CH_Matrix_Classes::Matrix C,D; CH_Matrix_Classes::genmult(A,B,C,1.,0.,1); CH_Matrix_Classes::genmult(C,C,D); CH_Matrix_Classes::Real d=2.*CH_Matrix_Classes::trace(D);
      CH_Matrix_Classes::genmult(A,A,C,1.,0.,1); CH_Matrix_Classes::genmult(B,B,D,1.,0.,1);
      return sqrt(2.*CH_Matrix_Classes::ip(C,D)+d);}

   virtual Coeffmat* subspace(const CH_Matrix_Classes::Matrix& P) const
    {CH_Matrix_Classes::Matrix C,D; CH_Matrix_Classes::genmult(P,A,C,1.,0.,1); CH_Matrix_Classes::genmult(P,B,D,1.,0.,1);
      return new CMlowrankdd(C,D,ConicBundle::clone(infop)); }
    
    virtual void multiply(CH_Matrix_Classes::Real d)
    { A*=d; if (infop) infop->multiply(d);}

    virtual CH_Matrix_Classes::Real ip(const CH_Matrix_Classes::Symmatrix& S) const
    { CH_Matrix_Classes::Matrix C; return 2.*CH_Matrix_Classes::ip(CH_Matrix_Classes::genmult(S,A,C),B); }
    
    virtual CH_Matrix_Classes::Real gramip(const CH_Matrix_Classes::Matrix& P) const
    { CH_Matrix_Classes::Matrix C,D; CH_Matrix_Classes::genmult(P,A,C,1.,0.,1); CH_Matrix_Classes::genmult(P,B,D,1.,0.,1);
      return 2.*CH_Matrix_Classes::ip(C,D); }

    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::Matrix C(A.coldim(),P.coldim());
      CH_Matrix_Classes::Real *cp=C.get_store();
      for(CH_Matrix_Classes::Integer j=0;j<C.coldim();j++){
        const CH_Matrix_Classes::Real *pp=P.get_store()+start_row+j*P.rowdim();
        const CH_Matrix_Classes::Real *ap=A.get_store();
        CH_Matrix_Classes::Integer ard=A.rowdim();
        for(CH_Matrix_Classes::Integer i=C.rowdim();--i>=0;){
          (*cp++)=CH_Matrix_Classes::mat_ip(ard,ap,pp);
          ap+=ard;
        }
      } 
      chk_set_init(C,1);
      CH_Matrix_Classes::Matrix D(B.coldim(),P.coldim());
      cp=D.get_store();
      {for(CH_Matrix_Classes::Integer j=0;j<D.coldim();j++){
        const CH_Matrix_Classes::Real *pp=P.get_store()+start_row+j*P.rowdim();
        const CH_Matrix_Classes::Real *ap=B.get_store();
        CH_Matrix_Classes::Integer ard=B.rowdim();
        for(CH_Matrix_Classes::Integer i=D.rowdim();--i>=0;){
          (*cp++)=CH_Matrix_Classes::mat_ip(ard,ap,pp);
          ap+=ard;
        }
      }}        
      chk_set_init(D,1);
      CH_Matrix_Classes::Real trval=0;
      if (Lam==0) {
        trval=CH_Matrix_Classes::ip(C,D);
      }
      else {
        assert(Lam->dim()==P.coldim());
        const CH_Matrix_Classes::Real *cp=C.get_store();
        const CH_Matrix_Classes::Real *dp=D.get_store();
        const CH_Matrix_Classes::Real *lp=Lam->get_store();
        for (CH_Matrix_Classes::Integer i=0;i<C.coldim();i++,cp+=C.rowdim(),dp+=D.rowdim())
          trval+=(*lp++)*CH_Matrix_Classes::mat_ip(C.rowdim(),cp,dp);
      }
      return 2.*trval; 
    }

    virtual void addmeto(CH_Matrix_Classes::Symmatrix& S,CH_Matrix_Classes::Real d=1.) const
    { CH_Matrix_Classes::rank2add(A,B,S,2.*d,1.); }

    virtual void addprodto(CH_Matrix_Classes::Matrix& D,const CH_Matrix_Classes::Matrix&C ,CH_Matrix_Classes::Real d=1.) const
    {CH_Matrix_Classes::Matrix E; CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(B,C,E,1.,0.,1),D,d,1.);
     CH_Matrix_Classes::genmult(B,CH_Matrix_Classes::genmult(A,C,E,1.,0.,1),D,d,1.);}

    virtual void addprodto(CH_Matrix_Classes::Matrix& D,const CH_Matrix_Classes::Sparsemat&C ,CH_Matrix_Classes::Real d=1.) const
    {CH_Matrix_Classes::Matrix E; CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(B,C,E,1.,0.,1),D,d,1.);
     CH_Matrix_Classes::genmult(B,CH_Matrix_Classes::genmult(A,C,E,1.,0.,1),D,d,1.);}

    virtual void left_right_prod(const CH_Matrix_Classes::Matrix& P,const CH_Matrix_Classes::Matrix& Q,CH_Matrix_Classes::Matrix& R) const 
    {
      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(B,Q,tmp2,1.,0.,1,0);
      CH_Matrix_Classes::genmult(tmp1,tmp2,R,1.,0.,0,0);
      CH_Matrix_Classes::genmult(P,B,tmp1,1.,0.,1,0);
      CH_Matrix_Classes::genmult(A,Q,tmp2,1.,0.,1,0);
      CH_Matrix_Classes::genmult(tmp1,tmp2,R,1.,1.,0,0);
    }

    virtual CH_Matrix_Classes::Integer prodvec_flops() const 
    { return 4*A.rowdim()*A.coldim()+4*B.rowdim()*B.coldim(); }

    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
    {return 2.*CH_Matrix_Classes::ip(A,S*B);}
    
    virtual void project(CH_Matrix_Classes::Symmatrix& S,const CH_Matrix_Classes::Matrix& P) const
    {CH_Matrix_Classes::Matrix C,D; CH_Matrix_Classes::genmult(P,A,C,1.,0.,1); CH_Matrix_Classes::genmult(P,B,D,1.,0.,1);
     CH_Matrix_Classes::rank2add(C,D,S,2.);}

    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 C(A.coldim(),P.coldim());
      CH_Matrix_Classes::Real *cp=C.get_store();
      for(CH_Matrix_Classes::Integer j=0;j<C.coldim();j++){
        const CH_Matrix_Classes::Real *pp=P.get_store()+start_row+j*P.rowdim();
        const CH_Matrix_Classes::Real *ap=A.get_store();
        CH_Matrix_Classes::Integer ard=A.rowdim();
        for(CH_Matrix_Classes::Integer i=C.rowdim();--i>=0;){
          (*cp++)=CH_Matrix_Classes::mat_ip(ard,ap,pp);
          ap+=ard;
        }
      } 
      chk_set_init(C,1);
      CH_Matrix_Classes::Matrix D(B.coldim(),P.coldim());
      cp=D.get_store();
      {for(CH_Matrix_Classes::Integer j=0;j<D.coldim();j++){
        const CH_Matrix_Classes::Real *pp=P.get_store()+start_row+j*P.rowdim();
        const CH_Matrix_Classes::Real *ap=B.get_store();
        CH_Matrix_Classes::Integer ard=B.rowdim();
        for(CH_Matrix_Classes::Integer i=D.rowdim();--i>=0;){
          (*cp++)=CH_Matrix_Classes::mat_ip(ard,ap,pp);
          ap+=ard;
        }
      }}        
      chk_set_init(D,1);
      CH_Matrix_Classes::rank2add(C,D,S,2.*alpha,1.,1); 
    }

   virtual const CH_Matrix_Classes::Matrix& postgenmult(const CH_Matrix_Classes::Matrix& D,CH_Matrix_Classes::Matrix& C,
                             CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int dtrans=0) const
    { 
      CH_Matrix_Classes::Matrix E; 
      CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(B,D,E,1.,0.,1,dtrans),C,alpha,beta);
      return CH_Matrix_Classes::genmult(B,CH_Matrix_Classes::genmult(A,D,E,1.,0.,1,dtrans),C,alpha,1.);
    }

    virtual const CH_Matrix_Classes::Matrix& pregenmult(const CH_Matrix_Classes::Matrix& D,CH_Matrix_Classes::Matrix& C,
                             CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int dtrans=0) const
    { 
      CH_Matrix_Classes::Matrix E; 
      CH_Matrix_Classes::genmult(CH_Matrix_Classes::genmult(D,A,E,1.,0.,dtrans),B,C,alpha,beta,0,1);
      return CH_Matrix_Classes::genmult(CH_Matrix_Classes::genmult(D,B,E,1.,0.,dtrans),A,C,alpha,1.,0,1);
    }

    virtual int equal(const Coeffmat* p,double tol=1e-6) const
    {
      const CMlowrankdd *pp=dynamic_cast<const CMlowrankdd *>(p);
      if (pp==0) 
        return 0;
      if ((A.rowdim()==(pp->A).rowdim())&&(A.coldim()==(pp->A).coldim()) &&
          (B.rowdim()==(pp->B).rowdim())&&(B.coldim()==(pp->B).coldim()) &&
          (CH_Matrix_Classes::norm2(A-pp->A)<tol) && (CH_Matrix_Classes::norm2(B-pp->B)<tol))
        return 1;
      if ((A.rowdim()==(pp->B).rowdim())&&(A.coldim()==(pp->B).coldim()) &&
          (B.rowdim()==(pp->A).rowdim())&&(B.coldim()==(pp->A).coldim()) &&
          (CH_Matrix_Classes::norm2(A-pp->B)<tol) && (CH_Matrix_Classes::norm2(B-pp->A)<tol))
        return 1;
      return 0;
    }

    virtual std::ostream& display(std::ostream& o) const 
    {o<<"CMlowrankdd\n";A.display(o);B.display(o);return o;}
    
    virtual std::ostream& out(std::ostream& o) const
    {return o<<"LOWRANK_DENSE_DENSE\n"<<A<<B;}

    virtual std::istream& in(std::istream& i)
    {
     i>>A>>B;
     if((A.rowdim()!=B.rowdim())||(A.coldim()!=B.coldim())){
         i.clear(std::ios::failbit);
         if (CH_Matrix_Classes::materrout) (*CH_Matrix_Classes::materrout)<<"*** ERROR: CMlowrankdd::in(): dimensions of A and B do not match"<<std::endl;
     }
     return i;
    }
    
    CMlowrankdd(std::istream& is,CoeffmatInfo* cip=0)
    {CM_type=CM_lowrankdd;infop=cip;in(is);}


};


}
#endif