Go to the documentation of this file.
13 #include "Coeffmat.hxx"
14 #include "CMsymdense.hxx"
16 namespace ConicBundle {
29  {
30  private:
33  bool positive;
34  public:
37  {A=Ain;di=sparseDiag(rowsip(A));positive=pos;CM_type=CM_gramsparsewd;infop=cip;}
39  virtual ~CMgramsparse_withoutdiag(){}
43  virtual Coeffmat* clone() const
44  {return new CMgramsparse_withoutdiag(A,positive,ConicBundle::clone(infop));}
47  virtual CH_Matrix_Classes::Integer dim() const { return A.rowdim(); }
51  {if (i==j) return 0.; return (positive? 1.:-1.)*CH_Matrix_Classes::ip(A.row(i),A.row(j));}
55  {if (positive) {
57  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
59  S(j,j)-=di.get_colval()(i);
60  }
61  }
62  else {
64  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
66  S(j,j)+=di.get_colval()(i);
67  }
68  }
69  }
72  virtual CH_Matrix_Classes::Real norm(void) const
74  for(CH_Matrix_Classes::Integer i=0;i<A.get_rowinfo().rowdim();i++){
75  CH_Matrix_Classes::genmult(A,A.row(A.get_rowinfo()(i,0)),vec,1.,0.,0,1);
76  sum+=CH_Matrix_Classes::ip(vec,vec)-di(A.get_rowinfo()(i,0),A.get_rowinfo()(i,0));
77  }
78  return std::sqrt(sum);
79  }
82  virtual Coeffmat* subspace(const CH_Matrix_Classes::Matrix& P) const
84  CH_Matrix_Classes::genmult(P,A,B,1.,0.,1);
87  CH_Matrix_Classes::genmult(P,di,B,1.,0.,1);
89  if (!positive) S*=-1;
90  return new CMsymdense(S,ConicBundle::clone(infop));
91  }
95  { di*=std::abs(d);
96  if (d<0.) {A*=sqrt(-d); positive=!positive;} else {A*=sqrt(d);}
97  if (infop) infop->multiply(d); }
106  {
108  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
110  sum+=CH_Matrix_Classes::ip(vec,vec)*di.get_colval()(i);
111  }
113  if (positive) return CH_Matrix_Classes::ip(B,B)-sum;
114  else return -CH_Matrix_Classes::ip(B,B)+sum;
115  }
119  {
120  CH_Matrix_Classes::Real trval=0.;
121  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
122  CH_Matrix_Classes::Matrix vec=P.row(di.get_colindex()(i)+start_row);
123  trval-=CH_Matrix_Classes::rowip(P,start_row+di.get_colindex()(i),Lam)*di.get_colval()(i);
124  }
127  const CH_Matrix_Classes::Integer brd=B.rowdim();
128  const CH_Matrix_Classes::Integer pcd=P.coldim();
129  const CH_Matrix_Classes::Integer prd=P.rowdim();
130  const CH_Matrix_Classes::Real *pp=P.get_store()+start_row;
131  const CH_Matrix_Classes::Indexmatrix &colinfo=A.get_colinfo();
132  const CH_Matrix_Classes::Real *vp=(A.get_colval()).get_store();
133  const CH_Matrix_Classes::Integer *ip=(A.get_colindex()).get_store();
134  for(CH_Matrix_Classes::Integer j=0;j<colinfo.rowdim();j++){
135  CH_Matrix_Classes::Integer indj=colinfo(j,0);
136  for(CH_Matrix_Classes::Integer i=colinfo(j,1);--i>=0;){
137  CH_Matrix_Classes::mat_xpeya(pcd,bp+indj,brd,pp+(*ip++),prd,(*vp++));
138  }
139  }
140  if (Lam==0) {
141  trval+=CH_Matrix_Classes::ip(B,B);
142  }
143  else {
144  assert(Lam->dim()==P.coldim());
145  const CH_Matrix_Classes::Real *bp=B.get_store();
146  const CH_Matrix_Classes::Real *lp=Lam->get_store();
147  for (CH_Matrix_Classes::Integer i=0;i<B.coldim();i++,bp+=B.rowdim())
148  trval+=(*lp++)*CH_Matrix_Classes::mat_ip(B.rowdim(),bp);
149  }
150  return (positive?trval:-trval);
151  }
155  {if (!positive) d*=-1;
156  CH_Matrix_Classes::rankadd(A,S,d,1.); S-=di*d;}
160  {CH_Matrix_Classes::Matrix D;if (!positive) d*=-1;
162  CH_Matrix_Classes::genmult(di,C,B,-d,1.,0);
163  }
167  { if (!positive) d*=-1;
170  CH_Matrix_Classes::genmult(di,C,B,-d,1.,0);}
174  {
175  R.init(P.coldim(),Q.coldim(),0.);
176  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
177  genmult(P.row(di.get_colindex()(i)),Q.row(di.get_colindex()(i)),R,di.get_colval()(i),1.,1);
178  }
179  CH_Matrix_Classes::Matrix tmp1; CH_Matrix_Classes::genmult(P,A,tmp1,1.,0.,1,0);
180  CH_Matrix_Classes::Matrix tmp2; CH_Matrix_Classes::genmult(A,Q,tmp2,1.,0.,1,0);
181  if (positive) CH_Matrix_Classes::genmult(tmp1,tmp2,R,1.,-1.,0,0);
182  else CH_Matrix_Classes::genmult(tmp1,tmp2,R,-1.,1.,0,0);
183  }
187  { return 4*A.nonzeros()+2*di.get_colindex().dim(); }
190  virtual int dense() const
191  {return 0;}
194  virtual int sparse() const
195  { return 0;}
200  CH_Matrix_Classes::Matrix& /* val */,
201  CH_Matrix_Classes::Real /* d=1. */)const
202  {return 0;}
205  virtual int support_in(const CH_Matrix_Classes::Sparsesym& /* S */) const
206  {return 0;}
210  {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);}
215  if (positive) CH_Matrix_Classes::rankadd(B,S); else CH_Matrix_Classes::rankadd(B,S,-1.);
216  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
217  rankadd(P.row(di.get_colindex()(i)),S,(positive?-1.:1.)*di.get_colval()(i),1.,1);
218  }
219  }
223  {
226  const CH_Matrix_Classes::Integer brd=B.rowdim();
227  const CH_Matrix_Classes::Integer pcd=P.coldim();
228  const CH_Matrix_Classes::Integer prd=P.rowdim();
229  const CH_Matrix_Classes::Real *pp=P.get_store()+start_row;
230  const CH_Matrix_Classes::Indexmatrix &colinfo=A.get_colinfo();
231  const CH_Matrix_Classes::Real *vp=(A.get_colval()).get_store();
232  const CH_Matrix_Classes::Integer *ip=(A.get_colindex()).get_store();
233  for(CH_Matrix_Classes::Integer j=0;j<colinfo.rowdim();j++){
234  CH_Matrix_Classes::Integer indj=colinfo(j,0);
235  for(CH_Matrix_Classes::Integer i=colinfo(j,1);--i>=0;){
236  CH_Matrix_Classes::mat_xpeya(pcd,bp+indj,brd,pp+(*ip++),prd,(*vp++));
237  }
238  }
239  if (positive) CH_Matrix_Classes::rankadd(B,S,alpha,1.,1); else CH_Matrix_Classes::rankadd(B,S,-alpha,1.,1);
240  for(CH_Matrix_Classes::Integer i=0;i<di.get_colindex().dim();i++){
241  rankadd(P.row(di.get_colindex()(i)+start_row),S,(positive?-alpha:alpha)*di.get_colval()(i),1.,1);
242  }
243  }
247  CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int btrans=0) const
248  {
250  CH_Matrix_Classes::genmult(A,CH_Matrix_Classes::genmult(A,B,D,1.,0.,1,btrans),C,(positive?1.:-1.)*alpha,beta);
251  return CH_Matrix_Classes::genmult(di,B,C,(positive?-1.:1.)*alpha,1.,btrans);
252  }
256  CH_Matrix_Classes::Real alpha=1.,CH_Matrix_Classes::Real beta=0.,int btrans=0) const
257  {
259  return CH_Matrix_Classes::genmult(CH_Matrix_Classes::genmult(B,A,D,1.,0.,btrans),A,C,(positive?1.:-1.)*alpha,beta,0,1);
260  return CH_Matrix_Classes::genmult(B,di,C,(positive?-1.:1.)*alpha,1.,btrans);
261  }
264  virtual int equal(const Coeffmat* p,double tol=1e-6) const
265  {
266  const CMgramsparse_withoutdiag *pp=dynamic_cast<const CMgramsparse_withoutdiag *>(p);
267  if (pp==0)
268  return 0;
269  return (positive==pp->positive) && CH_Matrix_Classes::equal(A,pp->A,tol);
270  }
273  virtual std::ostream& display(std::ostream& o) const
274  {o<<"CMgramsparse_withoutdiag\n"; A.display(o); return o;}
277  virtual std::ostream& out(std::ostream& o) const
278  {return o<<"GRAM_SPARSE_WITHOUTDIAG\n"<<positive<<"\n"<<A;}
281  virtual std::istream& in(std::istream& i)
282  { i>>positive>>A; di=sparseDiag(rowsip(A)); return i;}
285  CMgramsparse_withoutdiag(std::istream& is,CoeffmatInfo* cip=0)
286  {CM_type=CM_gramsparse;infop=cip;in(is);}
288  //--- specific routines
290  const CH_Matrix_Classes::Sparsemat& get_A() const {return A;}
292  bool get_positive() const {return positive;}
293 };
297 }
298 #endif
int Integer
all integer numbers in calculations and indexing are of this type
Definition: matop.hxx:40
Header declaring the classes ConicBundle::Coeffmat, ConicBundle::CoeffmatPointer, ConicBundle::Coeffm...
const Indexmatrix & get_colindex() const
returns the index vector of the column representation holding the row index for each element ...
Definition: sparssym.hxx:246
Integer rowdim() const
returns the row dimension
Definition: matrix.hxx:215
double Real
all real numbers in calculations are of this type
Definition: matop.hxx:50
CH_Matrix_Classes::Matrix & genmult(const MinorantBundle &A, const CH_Matrix_Classes::Matrix &B, CH_Matrix_Classes::Matrix &C, CH_Matrix_Classes::Real alpha=1., CH_Matrix_Classes::Real beta=0., int Atrans=0, int Btrans=0, CH_Matrix_Classes::Matrix *Coffset=0)
computes and returns C=alpha*A*B+beta*C where A and B may be transposed and A is considered to have t...
virtual int equal(const Coeffmat *p, double tol=1e-6) const
returns 1, if p is the same derived class and entries differ by less than tol, otherwise zero ...
Definition: CMgramsparse_withoutdiag.hxx:264
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
returns ip(*this,QQ^T)=trace Q^T(*this)Q for Q=P.rows(start_row,start_row+dim-1)
Definition: CMgramsparse_withoutdiag.hxx:118
allows to memorize the scalings applied to a Coeffmat and offers the basis for storing further user d...
Definition: Coeffmat.hxx:52
Matrix rowsip(const Matrix &A)
returns the row vector of the squared Frobenius norm of all rowd i of A, i.e., the sum of A(i...
Definition: matrix.hxx:1221
bool get_positive() const
returns the flag on whether the Gram matrix is used in positive or negative form
Definition: CMgramsparse_withoutdiag.hxx:292
CMgramsparse_withoutdiag(std::istream &is, CoeffmatInfo *cip=0)
constructor with istream and possibly additional user information
Definition: CMgramsparse_withoutdiag.hxx:285
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
computes C= alpha*B^(T if btrans)*(*this) + beta*C, C is also returned
Definition: CMgramsparse_withoutdiag.hxx:255
const Indexmatrix & get_colinfo() const
returns information on nonzero columns, k by 3, listing: index, %nonzeros, first index in colindex/co...
Definition: sparsmat.hxx:265
defines a base class for coefficient matrices in semidefinite programming, in particular for use with...
Definition: Coeffmat.hxx:125
void display(std::ostream &out, int precision=0, int width=0, int screenwidth=0) const
displays a matrix in a pretty way for bounded screen widths; for variables of value zero default valu...
Matrix class for integral values of type Integer
Definition: indexmat.hxx:195
virtual void addmeto(CH_Matrix_Classes::Symmatrix &S, CH_Matrix_Classes::Real d=1.) const
computes S+=d*(*this);
Definition: CMgramsparse_withoutdiag.hxx:154
Matrix & init(const Matrix &A, Real d=1., int atrans=0)
initialize to *this=A*d where A may be transposed
Definition: matrix.hxx:1035
Sparsesym sparseDiag(const Matrix &A, Real tol=SPARSE_ZERO_TOL)
forms a sparse symmetrix matrix having vector A on its diagonal
virtual CH_Matrix_Classes::Integer prodvec_flops() const
returns an estimate of number of flops to compute addprodto for a vector
Definition: CMgramsparse_withoutdiag.hxx:186
virtual std::ostream & display(std::ostream &o) const
display constraint information
Definition: CMgramsparse_withoutdiag.hxx:273
CMgramsparse_withoutdiag(const CH_Matrix_Classes::Sparsemat &Ain, bool pos=true, CoeffmatInfo *cip=0)
copy Ain, the flag for positive/negative and store the user information
Definition: CMgramsparse_withoutdiag.hxx:36
bool positive
if true use A*A^T-Diag, if false use -A*A^T+Diag
Definition: CMgramsparse_withoutdiag.hxx:33
virtual Coeffmat * clone() const
makes an explicit copy of itself and returns a pointer to it
Definition: CMgramsparse_withoutdiag.hxx:43
virtual CH_Matrix_Classes::Real operator()(CH_Matrix_Classes::Integer i, CH_Matrix_Classes::Integer j) const
returns the value of the matrix element (i,j)
Definition: CMgramsparse_withoutdiag.hxx:50
const Matrix & get_colval() const
returns the value vector of the column representation holding the value for each element ...
Definition: sparsmat.hxx:269
Real rowip(const Matrix &A, Integer i, const Matrix *scaling=0)
returns the squared Frobenius norm of row i of A, i.e., the sum of A(i,j)*A(i,j) over all j with poss...
Definition: matrix.hxx:1192
Real * get_store()
returns the current address of the internal value array; use cautiously, do not use delete! ...
Definition: matrix.hxx:326
Matrix class of symmetric matrices with real values of type Real
Definition: symmat.hxx:43
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
computes C= alpha*(*this)*B^(T if btrans) + beta*C, C is also returned
Definition: CMgramsparse_withoutdiag.hxx:246
virtual void addprodto(CH_Matrix_Classes::Matrix &B, const CH_Matrix_Classes::Sparsemat &C, CH_Matrix_Classes::Real d=1.) const
computes B+=d*(*this)*C
Definition: CMgramsparse_withoutdiag.hxx:166
conic bundle method solver for sum of convex functions. See the ConicBundle_Manual for a quick introd...
Definition: CBSolver.hxx:22
Matrix row(Integer i) const
returns row i copied to a new matrix
Coeffmattype CM_type
in order to enable type identification
Definition: Coeffmat.hxx:134
CoeffmatInfo * infop
allows the user to specify and output additional information
Definition: Coeffmat.hxx:135
virtual CH_Matrix_Classes::Integer dim() const
returns the order of the represented symmetric matrix
Definition: CMgramsparse_withoutdiag.hxx:47
implements a Gram matrix with zero diagonal as Coeffmat for a sparse rectangular CH_Matrix_Classes::...
Definition: CMgramsparse_withoutdiag.hxx:28
virtual void project(CH_Matrix_Classes::Symmatrix &S, const CH_Matrix_Classes::Matrix &P) const
computes S=P^T*(*this)*P
Definition: CMgramsparse_withoutdiag.hxx:213
CH_Matrix_Classes::Sparsesym di
this is the precomputed diagonal di=diag(A*A^T)
Definition: CMgramsparse_withoutdiag.hxx:32
implements a general dense symmetric Coeffmat based on CH_Matrix_Classes::Symmatrix (for use with Mat...
Definition: CMsymdense.hxx:27
Symmatrix & rankadd(const Matrix &A, Symmatrix &C, Real alpha=1., Real beta=0., int trans=0)
returns C=beta*C+alpha* A*A^T, where A may be transposed. If beta==0. then C is initiliazed to the co...
Matrix class of symmetric matrices with real values of type Real
Definition: sparssym.hxx:89
void multiply(CH_Matrix_Classes::Real sf)
scales the scale factor
Definition: Coeffmat.hxx:70
virtual int support_in(const CH_Matrix_Classes::Sparsesym &) const
returns 0 if the support of the costraint matrix is not contained in the support of the sparse symmet...
Definition: CMgramsparse_withoutdiag.hxx:205
Integer coldim() const
returns the column dimension
Definition: matrix.hxx:218
virtual int sparse() const
returns 0 if not sparse, otherwise 1
Definition: CMgramsparse_withoutdiag.hxx:194
for CMgramsparse
Definition: Coeffmat.hxx:37
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
computes S+=Q^T(*this)Q for Q=P.rows(start_row,start_row+dim-1)
Definition: CMgramsparse_withoutdiag.hxx:222
Matrix class for real values of type Real
Definition: matrix.hxx:74
virtual CH_Matrix_Classes::Real gramip(const CH_Matrix_Classes::Matrix &P) const
returns ip(*this,PP^T)=trace P^T(*this)P
Definition: CMgramsparse_withoutdiag.hxx:105
const Matrix & get_colval() const
returns the value vector of the column representation holding the value for each element ...
Definition: sparssym.hxx:248
Val mat_ip(Integer len, const Val *x, const Val *y, const Val *d=0)
return sum(x[i]*y[i]) summing over len elements of the arrays x and y.
Definition: matop.hxx:1096
Matrix class of sparse matrices with real values of type Real
Definition: sparsmat.hxx:74
void mat_xpeya(Integer len, Val *x, const Val *y, const Val a)
Set x[i]+=a*y[i] for len elements of the arrays x and y.
Definition: matop.hxx:543
CH_Matrix_Classes::Sparsemat A
the Coeffmat acts like A*A^T-Diag(A*A^T) or its negative
Definition: CMgramsparse_withoutdiag.hxx:31
virtual CH_Matrix_Classes::Real ip(const CH_Matrix_Classes::Symmatrix &S) const
returns ip(*this,S)=trace(*this*S), the trace inner product
Definition: CMgramsparse_withoutdiag.hxx:100
const Indexmatrix & get_rowinfo() const
returns information on nonzero rows, k by 3, listing: index, %nonzeros, first index in rowindex/rowva...
Definition: sparsmat.hxx:272
virtual void addprodto(CH_Matrix_Classes::Matrix &B, const CH_Matrix_Classes::Matrix &C, CH_Matrix_Classes::Real d=1.) const
computes B+=d*(*this)*C
Definition: CMgramsparse_withoutdiag.hxx:159
virtual void make_symmatrix(CH_Matrix_Classes::Symmatrix &S) const
returns a dense symmetric constraint matrix
Definition: CMgramsparse_withoutdiag.hxx:54
virtual Coeffmat * subspace(const CH_Matrix_Classes::Matrix &P) const
delivers a new object on the heap corresponding to the matrix P^T(*this)P, the caller is responsible ...
Definition: CMgramsparse_withoutdiag.hxx:82
const Indexmatrix & get_colindex() const
returns the index vector of the column representation holding the row index for each element ...
Definition: sparsmat.hxx:267
Integer rowdim() const
returns the row dimension
Definition: indexmat.hxx:321
virtual CH_Matrix_Classes::Real norm(void) const
returns the Frobenius norm of the matrix
Definition: CMgramsparse_withoutdiag.hxx:72
virtual int sparse(CH_Matrix_Classes::Indexmatrix &, CH_Matrix_Classes::Indexmatrix &, CH_Matrix_Classes::Matrix &, CH_Matrix_Classes::Real) const
returns 0 if not sparse. If it is sparse it returns 1 and the nonzero structure in I...
Definition: CMgramsparse_withoutdiag.hxx:198
Indexmatrix & genmult(const Indexmatrix &A, const Indexmatrix &B, Indexmatrix &C, Integer alpha=1, Integer beta=0, int atrans=0, int btrans=0)
returns C=beta*C+alpha*A*B, where A and B may be transposed; C must not be equal to A and B; if beta=...
virtual CH_Matrix_Classes::Real ip(const CH_Matrix_Classes::Sparsesym &S) const
returns the inner product of the constraint matrix with S
Definition: CMgramsparse_withoutdiag.hxx:209
Sparsemat row(Integer i) const
returns row i copied to a new sparse matrix
double abs(double d)
absolute value of a double
Definition: mymath.hxx:25
virtual int dense() const
returns 1 if its structure is as bad as its dense symmetric representation, otherwise 0 ...
Definition: CMgramsparse_withoutdiag.hxx:190
for CMgramsparse_withoutdiag
Definition: Coeffmat.hxx:39
Integer sum(const Indexmatrix &A)
returns the sum over all elements of A, i.e., (1 1 ... 1)*A*(1 1 ... 1)^T
virtual void left_right_prod(const CH_Matrix_Classes::Matrix &P, const CH_Matrix_Classes::Matrix &Q, CH_Matrix_Classes::Matrix &R) const
computes R=P^T*(*this)*Q
Definition: CMgramsparse_withoutdiag.hxx:173
void dim(Integer &_nr, Integer &_nc) const
returns the number of rows in _nr and the number of columns in _nc
Definition: indexmat.hxx:315
bool equal(const Matrix &A, const Matrix &B)
returns true if both matrices have the same size and the same elements
Definition: matrix.hxx:1433
Integer coldim() const
returns the column dimension
Definition: sparsmat.hxx:202
Integer rowdim() const
returns the row dimension
Definition: sparsmat.hxx:199
double sqrt(int a)
return sqrt for int a
Definition: mymath.hxx:121
Header declaring the class ConicBundle::CMsymdense (needed for ConicBundle::AffineMatrixFunction) ...
CoeffmatInfo * clone(const CoeffmatInfo *cip)
if cip is not zero, it calls and returns cip->clone() and 0 otherwise
Definition: Coeffmat.hxx:106
virtual std::ostream & out(std::ostream &o) const
put entire contents onto ostream with the class type in the beginning so that the derived class can b...
Definition: CMgramsparse_withoutdiag.hxx:277
Real ip(const Matrix &A, const Matrix &B)
returns the usual inner product of A and B, i.e., the sum of A(i,j)*B(i,j) over all i...
Definition: matrix.hxx:1165
virtual void multiply(CH_Matrix_Classes::Real d)
multiply constraint permanentely by d; this is to allow scaling or sign changes in the constraints ...
Definition: CMgramsparse_withoutdiag.hxx:94
Integer nonzeros() const
returns the number of nonzeros
Definition: sparsmat.hxx:205
const CH_Matrix_Classes::Sparsemat & get_A() const
returns the const reference to the internal matrix A forming the Gram matrix
Definition: CMgramsparse_withoutdiag.hxx:290
virtual std::istream & in(std::istream &i)
counterpart to out(), does not read the class type, though. This is assumed to have been read in orde...
Definition: CMgramsparse_withoutdiag.hxx:281