    // @(#)root/quadp:$Name:  $:$Id: TQpLinSolverBase.cxx,v 1.2 2004/05/24 12:45:40 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.                                        *
     *************************************************************************/

    //////////////////////////////////////////////////////////////////////////
    //                                                                      //
    // TQpLinSolverBase                                                     //
    //                                                                      //
    // Implementation of main solver for linear systems                     //
    //                                                                      //
    //////////////////////////////////////////////////////////////////////////

    #include "Riostream.h"
    #include "TQpLinSolverBase.h"
    #include "TMatrixD.h"

    ClassImp(TQpLinSolverBase)

    //______________________________________________________________________________
    TQpLinSolverBase::TQpLinSolverBase()
    {
      fNx   = 0;
      fMy   = 0;
      fMz   = 0;
      fNxup = 0;
      fNxlo = 0;
      fMcup = 0;
      fMclo = 0;
    }

    //______________________________________________________________________________
    TQpLinSolverBase::TQpLinSolverBase(TQpProbBase *factory,TQpDataBase *data)
    {
      fFactory = factory;

      fNx = data->fNx;
      fMy = data->fMy;
      fMz = data->fMz;

      fXloIndex.ResizeTo(data->fXloIndex); fXloIndex = data->fXloIndex;
      fXupIndex.ResizeTo(data->fXupIndex); fXupIndex = data->fXupIndex;
      fCloIndex.ResizeTo(data->fCloIndex); fCloIndex = data->fCloIndex;
      fCupIndex.ResizeTo(data->fCupIndex); fCupIndex = data->fCupIndex;

      fNxlo = fXloIndex.NonZeros();
      fNxup = fXupIndex.NonZeros();
      fMclo = fCloIndex.NonZeros();
      fMcup = fCupIndex.NonZeros();

      if (fNxup+fNxlo > 0) {
        fDd.ResizeTo(fNx);
        fDq.ResizeTo(fNx);
        data->GetDiagonalOfQ(fDq);
      }
      fNomegaInv.ResizeTo(fMz);
      fRhs      .ResizeTo(fNx+fMy+fMz);
    }

    //______________________________________________________________________________
    TQpLinSolverBase::TQpLinSolverBase(const TQpLinSolverBase &another) : TObject(another)
    {
      *this = another;
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::Factor(TQpDataBase * /* prob */,TQpVar *vars)
    {
      Assert(vars->ValidNonZeroPattern());

      if (fNxlo+fNxup > 0) {
        fDd.ResizeTo(fDq);
        fDd = fDq;
      }
      this->ComputeDiagonals(fDd,fNomegaInv,
    			 vars->fT,vars->fLambda,vars->fU,vars->fPi,
    			 vars->fV,vars->fGamma,vars->fW,vars->fPhi);
      if (fNxlo+fNxup > 0) this->PutXDiagonal(fDd);

      fNomegaInv.Invert();
      fNomegaInv *= -1.;

      if (fMclo+fMcup > 0) this->PutZDiagonal(fNomegaInv);
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::ComputeDiagonals(TVectorD &dd,TVectorD &omega,
                                            TVectorD &t, TVectorD &lambda,
                                            TVectorD &u, TVectorD &pi,
                                            TVectorD &v, TVectorD &gamma,
                                            TVectorD &w, TVectorD &phi)
    {
      if (fNxup+fNxlo > 0) {
        if (fNxlo > 0) AddElemDiv(dd,1.0,gamma,v,fXloIndex);
        if (fNxup > 0) AddElemDiv(dd,1.0,phi  ,w,fXupIndex);
      }
      omega.Zero();
      if (fMclo > 0) AddElemDiv(omega,1.0,lambda,t,fCloIndex);
      if (fMcup > 0) AddElemDiv(omega,1.0,pi,    u,fCupIndex);
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::Solve(TQpDataBase *prob,TQpVar *vars,TQpResidual *res,TQpVar *step)
    {
      Assert(vars->ValidNonZeroPattern());
      Assert(res ->ValidNonZeroPattern());

      (step->fX).ResizeTo(res->fRQ); step->fX = res->fRQ;
      if (fNxlo > 0) {
        TVectorD &vInvGamma = step->fV;
        vInvGamma.ResizeTo(vars->fGamma); vInvGamma = vars->fGamma;
        ElementDiv(vInvGamma,vars->fV,fXloIndex);

        AddElemMult(step->fX,1.0,vInvGamma,res->fRv);
        AddElemDiv (step->fX,1.0,res->fRgamma,vars->fV,fXloIndex);
      }

      if (fNxup > 0) {
        TVectorD &wInvPhi = step->fW;
        wInvPhi.ResizeTo(vars->fPhi); wInvPhi = vars->fPhi;
        ElementDiv(wInvPhi,vars->fW,fXupIndex);

        AddElemMult(step->fX,1.0,wInvPhi,res->fRw);
        AddElemDiv (step->fX,-1.0,res->fRphi,vars->fW,fXupIndex);
      }

      // start by partially computing step->fS
      (step->fS).ResizeTo(res->fRz); step->fS = res->fRz;
      if (fMclo > 0) {
        TVectorD &tInvLambda = step->fT;
        tInvLambda.ResizeTo(vars->fLambda); tInvLambda = vars->fLambda;
        ElementDiv(tInvLambda,vars->fT,fCloIndex);

        AddElemMult(step->fS,1.0,tInvLambda,res->fRt);
        AddElemDiv (step->fS,1.0,res->fRlambda,vars->fT,fCloIndex);
      }

      if (fMcup > 0) {
        TVectorD &uInvPi = step->fU;

        uInvPi.ResizeTo(vars->fPi); uInvPi = vars->fPi;
        ElementDiv(uInvPi,vars->fU,fCupIndex);

        AddElemMult(step->fS,1.0,uInvPi,res->fRu);
        AddElemDiv (step->fS,-1.0,res->fRpi,vars->fU,fCupIndex);
      }

      (step->fY).ResizeTo(res->fRA); step->fY = res->fRA;
      (step->fZ).ResizeTo(res->fRC); step->fZ = res->fRC;

      if (fMclo > 0)
        this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fLambda,prob);
      else
        this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fPi,prob);

      if (fMclo > 0) {
        (step->fT).ResizeTo(step->fS); step->fT = step->fS;
        Add(step->fT,-1.0,res->fRt);
        (step->fT).SelectNonZeros(fCloIndex);

        (step->fLambda).ResizeTo(res->fRlambda); step->fLambda = res->fRlambda;
        AddElemMult(step->fLambda,-1.0,vars->fLambda,step->fT);
        ElementDiv(step->fLambda,vars->fT,fCloIndex);
      }

      if (fMcup > 0) {
        (step->fU).ResizeTo(res->fRu); step->fU = res->fRu;
        Add(step->fU,-1.0,step->fS);
        (step->fU).SelectNonZeros(fCupIndex);

        (step->fPi).ResizeTo(res->fRpi); step->fPi = res->fRpi;
        AddElemMult(step->fPi,-1.0,vars->fPi,step->fU);
        ElementDiv(step->fPi,vars->fU,fCupIndex);
      }

      if (fNxlo > 0) {
        (step->fV).ResizeTo(step->fX); step->fV = step->fX;
        Add(step->fV,-1.0,res->fRv);
        (step->fV).SelectNonZeros(fXloIndex);

        (step->fGamma).ResizeTo(res->fRgamma); step->fGamma = res->fRgamma;
        AddElemMult(step->fGamma,-1.0,vars->fGamma,step->fV);
        ElementDiv(step->fGamma,vars->fV,fXloIndex);
      }

      if (fNxup > 0) {
        (step->fW).ResizeTo(res->fRw); step->fW = res->fRw;
        Add(step->fW,-1.0,step->fX);
        (step->fW).SelectNonZeros(fXupIndex);

        (step->fPhi).ResizeTo(res->fRphi); step->fPhi = res->fRphi;
        AddElemMult(step->fPhi,-1.0,vars->fPhi,step->fW);
        ElementDiv(step->fPhi,vars->fW,fXupIndex);
      }
      Assert(step->ValidNonZeroPattern());
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::SolveXYZS(TVectorD &stepx,TVectorD &stepy,
                                     TVectorD &stepz,TVectorD &steps,
                                     TVectorD & /* ztemp */, TQpDataBase * /* prob */ )
    {
      AddElemMult(stepz,-1.0,fNomegaInv,steps);
      this->JoinRHS(fRhs,stepx,stepy,stepz);

      this->SolveCompressed(fRhs);

      this->SeparateVars(stepx,stepy,stepz,fRhs);

      stepy *= -1.;
      stepz *= -1.;

      Add(steps,-1.0,stepz);
      ElementMult(steps,fNomegaInv);
      steps *= -1.;
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::JoinRHS(TVectorD &rhs_in, TVectorD &rhs1_in,
                                   TVectorD &rhs2_in,TVectorD &rhs3_in)
    {
      // joinRHS has to be delegated to the factory. This is true because
      // the rhs may be distributed across processors, so the factory is the
      // only object that knows with certainly how to scatter the elements.
      fFactory->JoinRHS(rhs_in,rhs1_in,rhs2_in,rhs3_in);
    }

    //______________________________________________________________________________
     void TQpLinSolverBase::SeparateVars(TVectorD &x_in,TVectorD &y_in,
                                        TVectorD &z_in,TVectorD &vars_in)
    {
      // separateVars has to be delegated to the factory. This is true because
      // the rhs may be distributed across processors, so the factory is the
      // only object that knows with certainly how to scatter the elements.
      fFactory->SeparateVars(x_in,y_in,z_in,vars_in);
    }

    //______________________________________________________________________________
    TQpLinSolverBase &TQpLinSolverBase::operator=(const TQpLinSolverBase &source)
    {
      if (this != &source) {
        TObject::operator=(source);

        fNx   = source.fNx;
        fMy   = source.fMy;
        fMz   = source.fMz;
        fNxup = source.fNxup;
        fNxlo = source.fNxlo;
        fMcup = source.fMcup;
        fMclo = source.fMclo;

        fNomegaInv.ResizeTo(source.fNomegaInv); fNomegaInv = source.fNomegaInv;
        fRhs      .ResizeTo(source.fRhs);       fRhs       = source.fRhs;

        fDd       .ResizeTo(source.fDd);        fDd        = source.fDd;
        fDq       .ResizeTo(source.fDq);        fDq        = source.fDq;

        fXupIndex .ResizeTo(source.fXupIndex);  fXupIndex  = source.fXupIndex;
        fCupIndex .ResizeTo(source.fCupIndex);  fCupIndex  = source.fCupIndex;
        fXloIndex .ResizeTo(source.fXloIndex);  fXloIndex  = source.fXloIndex;
        fCloIndex .ResizeTo(source.fCloIndex);  fCloIndex  = source.fCloIndex;
      }
      return *this;
    }