    // @(#)root/quadp:$Name:  $:$Id: TQpDataDens.cxx,v 1.4 2004/06/09 12:23:16 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;
    }