// @(#)root/quadp:$Name: $:$Id: TQpDataDens.cxx,v 1.3 2004/06/02 06:45:34 brun Exp $
// Author: Eddy Offermann May 2004
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
/*************************************************************************
* Parts of this file are copied from the OOQP distribution and *
* are subject to the following license: *
* *
* COPYRIGHT 2001 UNIVERSITY OF CHICAGO *
* *
* The copyright holder hereby grants you royalty-free rights to use, *
* reproduce, prepare derivative works, and to redistribute this software*
* to others, provided that any changes are clearly documented. This *
* software was authored by: *
* *
* E. MICHAEL GERTZ gertz@mcs.anl.gov *
* Mathematics and Computer Science Division *
* Argonne National Laboratory *
* 9700 S. Cass Avenue *
* Argonne, IL 60439-4844 *
* *
* STEPHEN J. WRIGHT swright@cs.wisc.edu *
* Computer Sciences Department *
* University of Wisconsin *
* 1210 West Dayton Street *
* Madison, WI 53706 FAX: (608)262-9777 *
* *
* Any questions or comments may be directed to one of the authors. *
* *
* ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES OF *
* ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT, AND *
* OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A CONTRACT *
* WITH THE DEPARTMENT OF ENERGY. *
*************************************************************************/
#include "Riostream.h"
#include "TQpDataDens.h"
//////////////////////////////////////////////////////////////////////////
// //
// TQpDataDens //
// //
// Data for the dense QP formulation //
// //
//////////////////////////////////////////////////////////////////////////
ClassImp(TQpDataDens)
//______________________________________________________________________________
TQpDataDens::TQpDataDens(Int_t nx,Int_t my,Int_t mz)
: TQpDataBase(nx,my,mz)
{
fQ.ResizeTo(fNx,fNx);
fA.ResizeTo(fMy,fNx);
fC.ResizeTo(fMz,fNx);
}
//______________________________________________________________________________
TQpDataDens::TQpDataDens(TVectorD &c_in, TMatrixDSym &Q_in,
TVectorD &xlow_in,TVectorD &ixlow_in,
TVectorD &xupp_in,TVectorD &ixupp_in,
TMatrixD &A_in, TVectorD &bA_in,
TMatrixD &C_in,
TVectorD &clow_in,TVectorD &iclow_in,
TVectorD &cupp_in,TVectorD &icupp_in)
{
fG .ResizeTo(c_in) ; fG = c_in;
fBa .ResizeTo(bA_in) ; fBa = bA_in;
fXloBound.ResizeTo(xlow_in) ; fXloBound = xlow_in;
fXloIndex.ResizeTo(ixlow_in); fXloIndex = ixlow_in;
fXupBound.ResizeTo(xupp_in) ; fXupBound = xupp_in;
fXupIndex.ResizeTo(ixupp_in); fXupIndex = ixupp_in;
fCloBound.ResizeTo(clow_in) ; fCloBound = clow_in;
fCloIndex.ResizeTo(iclow_in); fCloIndex = iclow_in;
fCupBound.ResizeTo(cupp_in) ; fCupBound = cupp_in;
fCupIndex.ResizeTo(icupp_in); fCupIndex = icupp_in;
fNx = fG.GetNrows();
fQ.Use(Q_in);
if (A_in.GetNrows() > 0) {
fA.Use(A_in);
fMy = fA.GetNrows();
} else
fMy = 0;
if (C_in.GetNrows() > 0) {
fC.Use(C_in);
fMz = fC.GetNrows();
} else
fMz = 0;
}
//______________________________________________________________________________
TQpDataDens::TQpDataDens(const TQpDataDens &another) : TQpDataBase(another)
{
*this = another;
}
//______________________________________________________________________________
void TQpDataDens::Qmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x )
{
y *= beta;
if (fQ.GetNoElements() > 0)
y += alpha*(fQ*x);
}
//______________________________________________________________________________
void TQpDataDens::Amult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
{
y *= beta;
if (fA.GetNoElements() > 0)
y += alpha*(fA*x);
}
//______________________________________________________________________________
void TQpDataDens::Cmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
{
y *= beta;
if (fC.GetNoElements() > 0)
y += alpha*(fC*x);
}
//______________________________________________________________________________
void TQpDataDens::ATransmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
{
y *= beta;
if (fA.GetNoElements() > 0)
y += alpha*(TMatrixD(TMatrixD::kTransposed,fA)*x);
}
//______________________________________________________________________________
void TQpDataDens::CTransmult(Double_t beta,TVectorD &y,Double_t alpha,const TVectorD &x)
{
y *= beta;
if (fC.GetNoElements() > 0)
y += alpha*(TMatrixD(TMatrixD::kTransposed,fC)*x);
}
//______________________________________________________________________________
Double_t TQpDataDens::DataNorm()
{
Double_t norm = 0.0;
Double_t componentNorm = fG.NormInf();
if (componentNorm > norm) norm = componentNorm;
TMatrixDSym fQ_abs(fQ);
componentNorm = (fQ_abs.Abs()).Max();
if (componentNorm > norm) norm = componentNorm;
componentNorm = fBa.NormInf();
if (componentNorm > norm) norm = componentNorm;
TMatrixD fA_abs(fQ);
componentNorm = (fA_abs.Abs()).Max();
if (componentNorm > norm) norm = componentNorm;
TMatrixD fC_abs(fQ);
componentNorm = (fC_abs.Abs()).Max();
if (componentNorm > norm) norm = componentNorm;
Assert(fXloBound.MatchesNonZeroPattern(fXloIndex));
componentNorm = fXloBound.NormInf();
if (componentNorm > norm) norm = componentNorm;
Assert(fXupBound.MatchesNonZeroPattern(fXupIndex));
componentNorm = fXupBound.NormInf();
if (componentNorm > norm) norm = componentNorm;
Assert(fCloBound.MatchesNonZeroPattern(fCloIndex));
componentNorm = fCloBound.NormInf();
if (componentNorm > norm) norm = componentNorm;
Assert(fCupBound.MatchesNonZeroPattern(fCupIndex));
componentNorm = fCupBound.NormInf();
if (componentNorm > norm) norm = componentNorm;
return norm;
}
//______________________________________________________________________________
void TQpDataDens::Print(Option_t * /*opt*/) const
{
fQ.Print("Q");
fG.Print("c");
fXloBound.Print("xlow");
fXloIndex.Print("ixlow");
fXupBound.Print("xupp");
fXupIndex.Print("ixupp");
fA.Print("A");
fBa.Print("b");
fC.Print("C");
fCloBound.Print("clow");
fCloIndex.Print("iclow");
fCupBound.Print("cupp");
fCupIndex.Print("icupp");
}
//______________________________________________________________________________
void TQpDataDens::PutQIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
{
m.SetSub(row,col,fQ);
}
//______________________________________________________________________________
void TQpDataDens::PutAIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
{
m.SetSub(row,col,fA);
}
//______________________________________________________________________________
void TQpDataDens::PutCIntoAt(TMatrixDBase &m,Int_t row,Int_t col)
{
m.SetSub(row,col,fC);
}
//______________________________________________________________________________
void TQpDataDens::GetDiagonalOfQ(TVectorD &dq)
{
const Int_t n = TMath::Min(fQ.GetNrows(),fQ.GetNcols());
dq.ResizeTo(n);
dq = TMatrixDDiag(fQ);
}
//______________________________________________________________________________
Double_t TQpDataDens::ObjectiveValue(TQpVar *vars)
{
TVectorD tmp(fG);
this->Qmult(1.0,tmp,0.5,vars->fX);
return tmp*vars->fX;
}
//______________________________________________________________________________
void TQpDataDens::DataRandom(TVectorD &x,TVectorD &y,TVectorD &z,TVectorD &s)
{
Double_t ix = 3074.20374;
TVectorD xdual(fNx);
this->RandomlyChooseBoundedVariables(x,xdual,fXloBound,fXloIndex,fXupBound,fXupIndex,ix,.25,.25,.25);
TVectorD sprime(fMz);
this->RandomlyChooseBoundedVariables(sprime,z,fCloBound,fCloIndex,fCupBound,fCupIndex,ix,.25,.25,.5);
fQ.RandomizePD(0.0,1.0,ix);
fA.Randomize(-10.0,10.0,ix);
fC.Randomize(-10.0,10.0,ix);
y .Randomize(-10.0,10.0,ix);
fG = xdual;
fG -= fQ*x;
fG += TMatrixD(TMatrixD::kTransposed,fA)*y;
fG += TMatrixD(TMatrixD::kTransposed,fC)*z;
fBa = fA*x;
s = fC*x;
// Now compute the real q = s-sprime
const TVectorD q = s-sprime;
// Adjust fCloBound and fCupBound appropriately
Add(fCloBound,1.0,q);
Add(fCupBound,1.0,q);
fCloBound.SelectNonZeros(fCloIndex);
fCupBound.SelectNonZeros(fCupIndex);
}
//______________________________________________________________________________
TQpDataDens &TQpDataDens::operator=(const TQpDataDens &source)
{
if (this != &source) {
TQpDataBase::operator=(source);
fQ.ResizeTo(source.fQ); fQ = source.fQ;
fA.ResizeTo(source.fA); fA = source.fA;
fC.ResizeTo(source.fC); fC = source.fC;
}
return *this;
}