// @(#)root/quadp:$Id$ // 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() { // Default constructor fNx = 0; fMy = 0; fMz = 0; fNxup = 0; fNxlo = 0; fMcup = 0; fMclo = 0; fFactory = 0; } //______________________________________________________________________________ TQpLinSolverBase::TQpLinSolverBase(TQpProbBase *factory,TQpDataBase *data) { // Constructor 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), fFactory(another.fFactory) { // Copy constructor *this = another; } //______________________________________________________________________________ void TQpLinSolverBase::Factor(TQpDataBase * /* prob */,TQpVar *vars) { // Sets up the matrix for the main linear system in "augmented system" form. The // actual factorization is performed by a routine specific to either the sparse // or dense case R__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) { // Computes the diagonal matrices in the augmented system from the current set of variables 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) { // Solves the system for a given set of residuals. Assembles the right-hand side appropriate // to the matrix factored in factor, solves the system using the factorization produced there, // partitions the solution vector into step components, then recovers the step components // eliminated during the block elimination that produced the augmented system form . R__ASSERT(vars->ValidNonZeroPattern()); R__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); } R__ASSERT(step->ValidNonZeroPattern()); } //______________________________________________________________________________ void TQpLinSolverBase::SolveXYZS(TVectorD &stepx,TVectorD &stepy, TVectorD &stepz,TVectorD &steps, TVectorD & /* ztemp */, TQpDataBase * /* prob */ ) { // Assemble right-hand side of augmented system and call SolveCompressed to solve it 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_out, TVectorD &rhs1_in, TVectorD &rhs2_in,TVectorD &rhs3_in) { // Assembles a single vector object from three given vectors . // rhs_out (output) final joined vector // rhs1_in (input) first part of rhs // rhs2_in (input) middle part of rhs // rhs3_in (input) last part of rhs . fFactory->JoinRHS(rhs_out,rhs1_in,rhs2_in,rhs3_in); } //______________________________________________________________________________ void TQpLinSolverBase::SeparateVars(TVectorD &x_in,TVectorD &y_in, TVectorD &z_in,TVectorD &vars_in) { // Extracts three component vectors from a given aggregated vector. // vars_in (input) aggregated vector // x_in (output) first part of vars // y_in (output) middle part of vars // z_in (output) last part of vars fFactory->SeparateVars(x_in,y_in,z_in,vars_in); } //______________________________________________________________________________ TQpLinSolverBase &TQpLinSolverBase::operator=(const TQpLinSolverBase &source) { // Assignment opeartor 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; // LM : copy also pointer data member fFactory = source.fFactory; } return *this; }