ROOT_528-00b_version: math/quadp/src/TQpLinSolverBase.cxx Source File

00001 // @(#)root/quadp:$Id: TQpLinSolverBase.cxx 36024 2010-10-01 16:03:30Z moneta $
00002 // Author: Eddy Offermann   May 2004
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 /*************************************************************************
00013  * Parts of this file are copied from the OOQP distribution and          *
00014  * are subject to the following license:                                 *
00015  *                                                                       *
00016  * COPYRIGHT 2001 UNIVERSITY OF CHICAGO                                  *
00017  *                                                                       *
00018  * The copyright holder hereby grants you royalty-free rights to use,    *
00019  * reproduce, prepare derivative works, and to redistribute this software*
00020  * to others, provided that any changes are clearly documented. This     *
00021  * software was authored by:                                             *
00022  *                                                                       *
00023  *   E. MICHAEL GERTZ      gertz@mcs.anl.gov                             *
00024  *   Mathematics and Computer Science Division                           *
00025  *   Argonne National Laboratory                                         *
00026  *   9700 S. Cass Avenue                                                 *
00027  *   Argonne, IL 60439-4844                                              *
00028  *                                                                       *
00029  *   STEPHEN J. WRIGHT     swright@cs.wisc.edu                           *
00030  *   Computer Sciences Department                                        *
00031  *   University of Wisconsin                                             *
00032  *   1210 West Dayton Street                                             *
00033  *   Madison, WI 53706   FAX: (608)262-9777                              *
00034  *                                                                       *
00035  * Any questions or comments may be directed to one of the authors.      *
00036  *                                                                       *
00037  * ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES OF   *
00038  * ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT, AND     *
00039  * OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A CONTRACT   *
00040  * WITH THE DEPARTMENT OF ENERGY.                                        *
00041  *************************************************************************/
00042 
00043 //////////////////////////////////////////////////////////////////////////
00044 //                                                                      //
00045 // TQpLinSolverBase                                                     //
00046 //                                                                      //
00047 // Implementation of main solver for linear systems                     //
00048 //                                                                      //
00049 //////////////////////////////////////////////////////////////////////////
00050 
00051 #include "Riostream.h"
00052 #include "TQpLinSolverBase.h"
00053 #include "TMatrixD.h"
00054 
00055 ClassImp(TQpLinSolverBase)
00056 
00057 //______________________________________________________________________________
00058 TQpLinSolverBase::TQpLinSolverBase()
00059 {
00060 // Default constructor
00061 
00062    fNx   = 0;
00063    fMy   = 0;
00064    fMz   = 0;
00065    fNxup = 0;
00066    fNxlo = 0;
00067    fMcup = 0;
00068    fMclo = 0;
00069    fFactory = 0;
00070 }
00071 
00072 
00073 //______________________________________________________________________________
00074 TQpLinSolverBase::TQpLinSolverBase(TQpProbBase *factory,TQpDataBase *data)
00075 {
00076 // Constructor
00077 
00078    fFactory = factory;
00079 
00080    fNx = data->fNx;
00081    fMy = data->fMy;
00082    fMz = data->fMz;
00083 
00084    fXloIndex.ResizeTo(data->fXloIndex); fXloIndex = data->fXloIndex;
00085    fXupIndex.ResizeTo(data->fXupIndex); fXupIndex = data->fXupIndex;
00086    fCloIndex.ResizeTo(data->fCloIndex); fCloIndex = data->fCloIndex;
00087    fCupIndex.ResizeTo(data->fCupIndex); fCupIndex = data->fCupIndex;
00088 
00089    fNxlo = fXloIndex.NonZeros();
00090    fNxup = fXupIndex.NonZeros();
00091    fMclo = fCloIndex.NonZeros();
00092    fMcup = fCupIndex.NonZeros();
00093 
00094    if (fNxup+fNxlo > 0) {
00095       fDd.ResizeTo(fNx);
00096       fDq.ResizeTo(fNx);
00097       data->GetDiagonalOfQ(fDq);
00098    }
00099    fNomegaInv.ResizeTo(fMz);
00100    fRhs      .ResizeTo(fNx+fMy+fMz);
00101 }
00102 
00103 
00104 //______________________________________________________________________________
00105 TQpLinSolverBase::TQpLinSolverBase(const TQpLinSolverBase &another) : TObject(another), 
00106                                                                       fFactory(another.fFactory)
00107 {
00108 // Copy constructor
00109 
00110    *this = another;
00111 }
00112 
00113 
00114 //______________________________________________________________________________
00115 void TQpLinSolverBase::Factor(TQpDataBase * /* prob */,TQpVar *vars)
00116 {
00117 // Sets up the matrix for the main linear system in "augmented system" form. The
00118 // actual factorization is performed by a routine specific to either the sparse
00119 // or dense case
00120 
00121    R__ASSERT(vars->ValidNonZeroPattern());
00122 
00123    if (fNxlo+fNxup > 0) {
00124       fDd.ResizeTo(fDq);
00125       fDd = fDq;
00126    }
00127    this->ComputeDiagonals(fDd,fNomegaInv,
00128       vars->fT,vars->fLambda,vars->fU,vars->fPi,
00129       vars->fV,vars->fGamma,vars->fW,vars->fPhi);
00130    if (fNxlo+fNxup > 0) this->PutXDiagonal(fDd);
00131 
00132    fNomegaInv.Invert();
00133    fNomegaInv *= -1.;
00134 
00135    if (fMclo+fMcup > 0) this->PutZDiagonal(fNomegaInv);
00136 }
00137 
00138 
00139 //______________________________________________________________________________
00140 void TQpLinSolverBase::ComputeDiagonals(TVectorD &dd,TVectorD &omega,
00141                                         TVectorD &t, TVectorD &lambda,
00142                                         TVectorD &u, TVectorD &pi,
00143                                         TVectorD &v, TVectorD &gamma,
00144                                         TVectorD &w, TVectorD &phi)
00145 {
00146 // Computes the diagonal matrices in the augmented system from the current set of variables
00147 
00148    if (fNxup+fNxlo > 0) {
00149       if (fNxlo > 0) AddElemDiv(dd,1.0,gamma,v,fXloIndex);
00150       if (fNxup > 0) AddElemDiv(dd,1.0,phi  ,w,fXupIndex);
00151    }
00152    omega.Zero();
00153    if (fMclo > 0) AddElemDiv(omega,1.0,lambda,t,fCloIndex);
00154    if (fMcup > 0) AddElemDiv(omega,1.0,pi,    u,fCupIndex);
00155 }
00156 
00157 
00158 //______________________________________________________________________________
00159 void TQpLinSolverBase::Solve(TQpDataBase *prob,TQpVar *vars,TQpResidual *res,TQpVar *step)
00160 {
00161 // Solves the system for a given set of residuals. Assembles the right-hand side appropriate
00162 // to the matrix factored in factor, solves the system using the factorization produced there,
00163 // partitions the solution vector into step components, then recovers the step components
00164 // eliminated during the block elimination that produced the augmented system form .
00165 
00166    R__ASSERT(vars->ValidNonZeroPattern());
00167    R__ASSERT(res ->ValidNonZeroPattern());
00168 
00169    (step->fX).ResizeTo(res->fRQ); step->fX = res->fRQ;
00170    if (fNxlo > 0) {
00171       TVectorD &vInvGamma = step->fV;
00172       vInvGamma.ResizeTo(vars->fGamma); vInvGamma = vars->fGamma;
00173       ElementDiv(vInvGamma,vars->fV,fXloIndex);
00174 
00175       AddElemMult(step->fX,1.0,vInvGamma,res->fRv);
00176       AddElemDiv (step->fX,1.0,res->fRgamma,vars->fV,fXloIndex);
00177    }
00178 
00179    if (fNxup > 0) {
00180       TVectorD &wInvPhi = step->fW;
00181       wInvPhi.ResizeTo(vars->fPhi); wInvPhi = vars->fPhi;
00182       ElementDiv(wInvPhi,vars->fW,fXupIndex);
00183 
00184       AddElemMult(step->fX,1.0,wInvPhi,res->fRw);
00185       AddElemDiv (step->fX,-1.0,res->fRphi,vars->fW,fXupIndex);
00186    }
00187 
00188    // start by partially computing step->fS
00189    (step->fS).ResizeTo(res->fRz); step->fS = res->fRz;
00190    if (fMclo > 0) {
00191       TVectorD &tInvLambda = step->fT;
00192       tInvLambda.ResizeTo(vars->fLambda); tInvLambda = vars->fLambda;
00193       ElementDiv(tInvLambda,vars->fT,fCloIndex);
00194 
00195       AddElemMult(step->fS,1.0,tInvLambda,res->fRt);
00196       AddElemDiv (step->fS,1.0,res->fRlambda,vars->fT,fCloIndex);
00197    }
00198 
00199    if (fMcup > 0) {
00200       TVectorD &uInvPi = step->fU;
00201 
00202       uInvPi.ResizeTo(vars->fPi); uInvPi = vars->fPi;
00203       ElementDiv(uInvPi,vars->fU,fCupIndex);
00204 
00205       AddElemMult(step->fS,1.0,uInvPi,res->fRu);
00206       AddElemDiv (step->fS,-1.0,res->fRpi,vars->fU,fCupIndex);
00207    }
00208 
00209    (step->fY).ResizeTo(res->fRA); step->fY = res->fRA;
00210    (step->fZ).ResizeTo(res->fRC); step->fZ = res->fRC;
00211 
00212    if (fMclo > 0)
00213       this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fLambda,prob);
00214    else
00215       this->SolveXYZS(step->fX,step->fY,step->fZ,step->fS,step->fPi,prob);
00216 
00217    if (fMclo > 0) {
00218       (step->fT).ResizeTo(step->fS); step->fT = step->fS;
00219       Add(step->fT,-1.0,res->fRt);
00220       (step->fT).SelectNonZeros(fCloIndex);
00221 
00222       (step->fLambda).ResizeTo(res->fRlambda); step->fLambda = res->fRlambda;
00223       AddElemMult(step->fLambda,-1.0,vars->fLambda,step->fT);
00224       ElementDiv(step->fLambda,vars->fT,fCloIndex);
00225    }
00226 
00227    if (fMcup > 0) {
00228       (step->fU).ResizeTo(res->fRu); step->fU = res->fRu;
00229       Add(step->fU,-1.0,step->fS);
00230       (step->fU).SelectNonZeros(fCupIndex);
00231 
00232       (step->fPi).ResizeTo(res->fRpi); step->fPi = res->fRpi;
00233       AddElemMult(step->fPi,-1.0,vars->fPi,step->fU);
00234       ElementDiv(step->fPi,vars->fU,fCupIndex);
00235    }
00236 
00237    if (fNxlo > 0) {
00238       (step->fV).ResizeTo(step->fX); step->fV = step->fX;
00239       Add(step->fV,-1.0,res->fRv);
00240       (step->fV).SelectNonZeros(fXloIndex);
00241 
00242       (step->fGamma).ResizeTo(res->fRgamma); step->fGamma = res->fRgamma;
00243       AddElemMult(step->fGamma,-1.0,vars->fGamma,step->fV);
00244       ElementDiv(step->fGamma,vars->fV,fXloIndex);
00245    }
00246 
00247    if (fNxup > 0) {
00248       (step->fW).ResizeTo(res->fRw); step->fW = res->fRw;
00249       Add(step->fW,-1.0,step->fX);
00250       (step->fW).SelectNonZeros(fXupIndex);
00251 
00252       (step->fPhi).ResizeTo(res->fRphi); step->fPhi = res->fRphi;
00253       AddElemMult(step->fPhi,-1.0,vars->fPhi,step->fW);
00254       ElementDiv(step->fPhi,vars->fW,fXupIndex);
00255    }
00256    R__ASSERT(step->ValidNonZeroPattern());
00257 }
00258 
00259 
00260 //______________________________________________________________________________
00261 void TQpLinSolverBase::SolveXYZS(TVectorD &stepx,TVectorD &stepy,
00262                                  TVectorD &stepz,TVectorD &steps,
00263                                  TVectorD & /* ztemp */, TQpDataBase * /* prob */ )
00264 {
00265 // Assemble right-hand side of augmented system and call SolveCompressed to solve it
00266 
00267    AddElemMult(stepz,-1.0,fNomegaInv,steps);
00268    this->JoinRHS(fRhs,stepx,stepy,stepz);
00269 
00270    this->SolveCompressed(fRhs);
00271 
00272    this->SeparateVars(stepx,stepy,stepz,fRhs);
00273 
00274    stepy *= -1.;
00275    stepz *= -1.;
00276 
00277    Add(steps,-1.0,stepz);
00278    ElementMult(steps,fNomegaInv);
00279    steps *= -1.;
00280 }
00281 
00282 
00283 //______________________________________________________________________________
00284 void TQpLinSolverBase::JoinRHS(TVectorD &rhs_out, TVectorD &rhs1_in,
00285                                TVectorD &rhs2_in,TVectorD &rhs3_in)
00286 {
00287 // Assembles a single vector object from three given vectors .
00288 //     rhs_out (output) final joined vector
00289 //     rhs1_in (input) first part of rhs
00290 //     rhs2_in (input) middle part of rhs
00291 //     rhs3_in (input) last part of rhs .
00292 
00293    fFactory->JoinRHS(rhs_out,rhs1_in,rhs2_in,rhs3_in);
00294 }
00295 
00296 
00297 //______________________________________________________________________________
00298 void TQpLinSolverBase::SeparateVars(TVectorD &x_in,TVectorD &y_in,
00299                                     TVectorD &z_in,TVectorD &vars_in)
00300 {
00301 // Extracts three component vectors from a given aggregated vector.
00302 //     vars_in  (input) aggregated vector
00303 //     x_in (output) first part of vars
00304 //     y_in (output) middle part of vars
00305 //     z_in (output) last part of vars
00306 
00307    fFactory->SeparateVars(x_in,y_in,z_in,vars_in);
00308 }
00309 
00310 
00311 //______________________________________________________________________________
00312 TQpLinSolverBase &TQpLinSolverBase::operator=(const TQpLinSolverBase &source)
00313 {
00314 // Assignment opeartor
00315 
00316    if (this != &source) {
00317       TObject::operator=(source);
00318 
00319       fNx   = source.fNx;
00320       fMy   = source.fMy;
00321       fMz   = source.fMz;
00322       fNxup = source.fNxup;
00323       fNxlo = source.fNxlo;
00324       fMcup = source.fMcup;
00325       fMclo = source.fMclo;
00326 
00327       fNomegaInv.ResizeTo(source.fNomegaInv); fNomegaInv = source.fNomegaInv;
00328       fRhs      .ResizeTo(source.fRhs);       fRhs       = source.fRhs;
00329 
00330       fDd       .ResizeTo(source.fDd);        fDd        = source.fDd;
00331       fDq       .ResizeTo(source.fDq);        fDq        = source.fDq;
00332 
00333       fXupIndex .ResizeTo(source.fXupIndex);  fXupIndex  = source.fXupIndex;
00334       fCupIndex .ResizeTo(source.fCupIndex);  fCupIndex  = source.fCupIndex;
00335       fXloIndex .ResizeTo(source.fXloIndex);  fXloIndex  = source.fXloIndex;
00336       fCloIndex .ResizeTo(source.fCloIndex);  fCloIndex  = source.fCloIndex;
00337 
00338       // LM : copy also pointer data member
00339       fFactory = source.fFactory;
00340    }
00341    return *this;
00342 }