Epetra Package Browser (Single Doxygen Collection) Development
Epetra_SerialDenseSolver.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include "Epetra_SerialDenseSolver.h"
00044 #include "Epetra_SerialDenseMatrix.h"
00045 
00046 //=============================================================================
00047 Epetra_SerialDenseSolver::Epetra_SerialDenseSolver()
00048   : Epetra_CompObject(),
00049     Epetra_BLAS(),
00050     Epetra_LAPACK(),
00051     Equilibrate_(false),
00052     ShouldEquilibrate_(false),
00053     A_Equilibrated_(false),
00054     B_Equilibrated_(false),
00055     Transpose_(false),
00056     Factored_(false),
00057     EstimateSolutionErrors_(false),
00058     SolutionErrorsEstimated_(false),
00059     Solved_(false),
00060     Inverted_(false),
00061     ReciprocalConditionEstimated_(false),
00062     RefineSolution_(false),
00063     SolutionRefined_(false),
00064     TRANS_('N'),
00065     M_(0),
00066     N_(0),
00067     Min_MN_(0),
00068     NRHS_(0),
00069     LDA_(0),
00070     LDAF_(0),
00071     LDB_(0),
00072     LDX_(0),
00073     INFO_(0),
00074     LWORK_(0),
00075     IPIV_(0),
00076     IWORK_(0),
00077     ANORM_(0.0),
00078     RCOND_(0.0),
00079     ROWCND_(0.0),
00080     COLCND_(0.0),
00081     AMAX_(0.0),
00082     Matrix_(0),
00083     LHS_(0),
00084     RHS_(0),
00085     Factor_(0),
00086     A_(0),
00087     FERR_(0),
00088     BERR_(0),
00089     AF_(0),
00090     WORK_(0),
00091     R_(0),
00092     C_(0),
00093     B_(0),
00094     X_(0)
00095 {
00096   InitPointers();
00097   ResetMatrix();
00098   ResetVectors();
00099 }
00100 //=============================================================================
00101 Epetra_SerialDenseSolver::~Epetra_SerialDenseSolver()
00102 {
00103   DeleteArrays();
00104 }
00105 //=============================================================================
00106 void Epetra_SerialDenseSolver::InitPointers()
00107 {
00108   IWORK_ = 0;
00109   FERR_ = 0;
00110   BERR_ = 0;
00111   Factor_ =0;
00112   Matrix_ =0;
00113   AF_ = 0;
00114   IPIV_ = 0;
00115   WORK_ = 0;
00116   R_ = 0;
00117   C_ = 0;
00118   INFO_ = 0;
00119   LWORK_ = 0;    
00120 }
00121 //=============================================================================
00122 void Epetra_SerialDenseSolver::DeleteArrays()
00123 {
00124   if (IWORK_ != 0) {delete [] IWORK_; IWORK_ = 0;}
00125   if (FERR_ != 0)  {delete [] FERR_; FERR_ = 0;}
00126   if (BERR_ != 0)  {delete [] BERR_; BERR_ = 0;}
00127   if (Factor_ != Matrix_ && Factor_ != 0)   {delete Factor_; Factor_ = 0;}
00128   if (Factor_ !=0) Factor_ = 0;
00129   if (AF_ !=0) AF_ = 0;
00130   if (IPIV_ != 0)  {delete [] IPIV_;IPIV_ = 0;}
00131   if (WORK_ != 0)  {delete [] WORK_;WORK_ = 0;}
00132   if (R_ != 0 && R_ != C_)     {delete [] R_;R_ = 0;}
00133   if (R_ != 0) R_ = 0;
00134   if (C_ != 0)     {delete [] C_;C_ = 0;}
00135   INFO_ = 0;
00136   LWORK_ = 0;    
00137 }
00138 //=============================================================================
00139 void Epetra_SerialDenseSolver::ResetMatrix()
00140 {
00141   DeleteArrays();
00142   ResetVectors();
00143   Matrix_ = 0;
00144   Factor_ = 0;
00145   A_Equilibrated_ = false;
00146   Factored_ = false;
00147   Inverted_ = false;
00148   M_ = 0;
00149   N_ = 0;
00150   Min_MN_ = 0;
00151   LDA_ = 0;
00152   LDAF_ = 0;
00153   ANORM_ = -1.0;
00154   RCOND_ = -1.0;
00155   ROWCND_ = -1.0;
00156   COLCND_ = -1.0;
00157   AMAX_ = -1.0;
00158   A_ = 0;
00159 
00160 }
00161 //=============================================================================
00162 int Epetra_SerialDenseSolver::SetMatrix(Epetra_SerialDenseMatrix & A_in) {
00163   ResetMatrix();
00164   Matrix_ = &A_in;
00165   Factor_ = &A_in;
00166   M_ = A_in.M();
00167   N_ = A_in.N();
00168   Min_MN_ = EPETRA_MIN(M_,N_);
00169   LDA_ = A_in.LDA();
00170   LDAF_ = LDA_;
00171   A_ = A_in.A();
00172   AF_ = A_in.A();
00173   return(0);
00174 }
00175 //=============================================================================
00176 void Epetra_SerialDenseSolver::ResetVectors()
00177 {
00178   LHS_ = 0;
00179   RHS_ = 0;
00180   B_ = 0;
00181   X_ = 0;
00182   ReciprocalConditionEstimated_ = false;
00183   SolutionRefined_ = false;
00184   Solved_ = false;
00185   SolutionErrorsEstimated_ = false;
00186   B_Equilibrated_ = false;
00187   NRHS_ = 0;
00188   LDB_ = 0;
00189   LDX_ = 0;
00190 }
00191 //=============================================================================
00192 int Epetra_SerialDenseSolver::SetVectors(Epetra_SerialDenseMatrix & X_in, Epetra_SerialDenseMatrix & B_in)
00193 {
00194   if (B_in.M()!=X_in.M() || B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
00195   if (B_in.A()==0) EPETRA_CHK_ERR(-2);
00196   if (B_in.LDA()<1) EPETRA_CHK_ERR(-3);
00197   if (X_in.A()==0) EPETRA_CHK_ERR(-4);
00198   if (X_in.LDA()<1) EPETRA_CHK_ERR(-5);
00199 
00200   ResetVectors(); 
00201   LHS_ = &X_in;
00202   RHS_ = &B_in;
00203   NRHS_ = B_in.N();
00204 
00205   B_ = B_in.A();
00206   LDB_ = B_in.LDA();
00207   X_ = X_in.A();
00208   LDX_ = X_in.LDA();
00209   return(0);
00210 }
00211 //=============================================================================
00212 void Epetra_SerialDenseSolver::EstimateSolutionErrors(bool Flag) {
00213   EstimateSolutionErrors_ = Flag;
00214   // If the errors are estimated, this implies that the solution must be refined
00215   RefineSolution_ = RefineSolution_ || Flag;
00216   return;
00217 }
00218 //=============================================================================
00219 int Epetra_SerialDenseSolver::Factor(void) {
00220   if (Factored()) return(0); // Already factored
00221   if (Inverted()) EPETRA_CHK_ERR(-100); // Cannot factor inverted matrix
00222   int ierr = 0;
00223 
00224   ANORM_ = Matrix_->OneNorm(); // Compute 1-Norm of A
00225 
00226 
00227   // If we want to refine the solution, then the factor must
00228   // be stored separatedly from the original matrix
00229 
00230   if (A_ == AF_)
00231     if (RefineSolution_ ) {
00232       Factor_ = new Epetra_SerialDenseMatrix(*Matrix_);
00233       AF_ = Factor_->A();
00234       LDAF_ = Factor_->LDA();
00235     }
00236   
00237   if (Equilibrate_) ierr = EquilibrateMatrix();
00238 
00239   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00240   
00241   if (IPIV_==0) IPIV_ = new int[Min_MN_]; // Allocated Pivot vector if not already done.
00242 
00243   GETRF (M_, N_, AF_, LDAF_, IPIV_, &INFO_);
00244 
00245   Factored_ = true;
00246   double DN = N_;
00247   UpdateFlops(2.0*(DN*DN*DN)/3.0);
00248 
00249   EPETRA_CHK_ERR(INFO_);
00250   return(0);
00251 
00252 }
00253 
00254 //=============================================================================
00255 int Epetra_SerialDenseSolver::Solve(void) {
00256   int ierr = 0;
00257 
00258   // We will call one of four routines depending on what services the user wants and 
00259   // whether or not the matrix has been inverted or factored already.
00260   //
00261   // If the matrix has been inverted, use DGEMM to compute solution.
00262   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00263   // call the X interface.
00264   // Otherwise, if the matrix is already factored we will call the TRS interface.
00265   // Otherwise, if the matrix is unfactored we will call the SV interface.
00266 
00267 
00268 
00269   if (Equilibrate_) {
00270     ierr = EquilibrateRHS();
00271     B_Equilibrated_ = true;
00272   }
00273   EPETRA_CHK_ERR(ierr);
00274   if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
00275   if (!A_Equilibrated_ && B_Equilibrated_) EPETRA_CHK_ERR(-2);
00276   if (B_==0) EPETRA_CHK_ERR(-3); // No B
00277   if (X_==0) EPETRA_CHK_ERR(-4); // No X
00278 
00279   if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
00280 
00281   double DN = N_;
00282   double DNRHS = NRHS_;
00283   if (Inverted()) {
00284 
00285     if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
00286 
00287     GEMM(TRANS_, 'N', N_, NRHS_, N_, 1.0, AF_, LDAF_,
00288     B_, LDB_, 0.0, X_, LDX_);
00289     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00290     UpdateFlops(2.0*DN*DN*DNRHS);
00291     Solved_ = true;
00292   }
00293   else {
00294 
00295     if (!Factored()) Factor(); // Matrix must be factored
00296     
00297     if (B_!=X_) {
00298        *LHS_ = *RHS_; // Copy B to X if needed
00299        X_ = LHS_->A(); LDX_ = LHS_->LDA();
00300     }
00301     GETRS(TRANS_, N_, NRHS_, AF_, LDAF_, IPIV_, X_, LDX_, &INFO_);
00302     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00303     UpdateFlops(2.0*DN*DN*DNRHS);
00304     Solved_ = true;
00305 
00306   }
00307   int ierr1=0;
00308   if (RefineSolution_ && !Inverted()) ierr1 = ApplyRefinement();
00309   if (ierr1!=0) EPETRA_CHK_ERR(ierr1)
00310   else
00311     EPETRA_CHK_ERR(ierr);
00312   
00313   if (Equilibrate_) ierr1 = UnequilibrateLHS();
00314   EPETRA_CHK_ERR(ierr1);
00315   return(0);
00316 }
00317 //=============================================================================
00318 int Epetra_SerialDenseSolver::ApplyRefinement(void)
00319 {
00320   double DN = N_;
00321   double DNRHS = NRHS_;
00322   if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
00323   if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
00324 
00325   if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
00326   FERR_ = new double[NRHS_];
00327   if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of BERR_, since NRHS_ may change
00328   BERR_ = new double[NRHS_];
00329   AllocateWORK();
00330   AllocateIWORK();
00331   
00332   GERFS(TRANS_, N_, NRHS_, A_, LDA_, AF_, LDAF_, IPIV_,
00333          B_, LDB_, X_, LDX_, FERR_, BERR_, 
00334          WORK_, IWORK_, &INFO_);
00335   
00336   
00337   SolutionErrorsEstimated_ = true;
00338   ReciprocalConditionEstimated_ = true;
00339   SolutionRefined_ = true;
00340   
00341   UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
00342   
00343   EPETRA_CHK_ERR(INFO_);
00344   return(0);
00345 
00346 }
00347 
00348 //=============================================================================
00349 int Epetra_SerialDenseSolver::ComputeEquilibrateScaling(void) {
00350   if (R_!=0) return(0); // Already computed
00351  
00352   double DM = M_;
00353   double DN = N_;
00354   R_ = new double[M_];
00355   C_ = new double[N_];
00356   
00357   GEEQU (M_, N_, AF_, LDAF_, R_, C_, &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00358   if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
00359 
00360   if (COLCND_<0.1 || ROWCND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
00361 
00362   UpdateFlops(4.0*DM*DN);
00363   
00364   return(0);
00365 }
00366 
00367 //=============================================================================
00368 int Epetra_SerialDenseSolver::EquilibrateMatrix(void)
00369 {
00370   int i, j;
00371   int ierr = 0;
00372 
00373   double DN = N_;
00374   double DM = M_;
00375 
00376   if (A_Equilibrated_) return(0); // Already done
00377   if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
00378   if (ierr!=0) EPETRA_CHK_ERR(ierr);
00379   if (A_==AF_) {
00380     double * ptr;
00381     for (j=0; j<N_; j++) {
00382       ptr = A_ + j*LDA_;
00383       double s1 = C_[j];
00384       for (i=0; i<M_; i++) {
00385   *ptr = *ptr*s1*R_[i];
00386   ptr++;
00387       }
00388     }
00389     UpdateFlops(2.0*DM*DN);
00390   }
00391   else {
00392     double * ptr;
00393     double * ptr1;
00394     for (j=0; j<N_; j++) {
00395       ptr = A_ + j*LDA_;
00396       ptr1 = AF_ + j*LDAF_;
00397       double s1 = C_[j];
00398       for (i=0; i<M_; i++) {
00399   *ptr = *ptr*s1*R_[i];
00400   ptr++;
00401   *ptr1 = *ptr1*s1*R_[i];
00402   ptr1++;
00403       }
00404     }
00405     UpdateFlops(4.0*DM*DN);
00406   }
00407   
00408   A_Equilibrated_ = true;
00409   
00410   return(0);
00411 }
00412 
00413 //=============================================================================
00414 int Epetra_SerialDenseSolver::EquilibrateRHS(void)
00415 {
00416   int i, j;
00417   int ierr = 0;
00418 
00419   if (B_Equilibrated_) return(0); // Already done
00420   if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
00421   if (ierr!=0) EPETRA_CHK_ERR(ierr);
00422 
00423   double * R_tmp = R_;
00424   if (Transpose_) R_tmp = C_;
00425 
00426   double * ptr;
00427   for (j=0; j<NRHS_; j++) {
00428     ptr = B_ + j*LDB_;
00429     for (i=0; i<M_; i++) {
00430       *ptr = *ptr*R_tmp[i];
00431       ptr++;
00432     }
00433   }
00434 
00435   
00436   B_Equilibrated_ = true;
00437   UpdateFlops((double) N_*(double) NRHS_);
00438   
00439   return(0);
00440 }
00441 
00442 //=============================================================================
00443 int Epetra_SerialDenseSolver::UnequilibrateLHS(void)
00444 {
00445   int i, j;
00446 
00447   if (!B_Equilibrated_) return(0); // Nothing to do
00448 
00449   double * C_tmp = C_;
00450   if (Transpose_) C_tmp = R_;
00451 
00452   double * ptr;
00453   for (j=0; j<NRHS_; j++) {
00454     ptr = X_ + j*LDX_;
00455     for (i=0; i<N_; i++) {
00456       *ptr = *ptr*C_tmp[i];
00457       ptr++;
00458     }
00459   }
00460 
00461   
00462   UpdateFlops((double) N_ *(double) NRHS_);
00463   
00464   return(0);
00465 }
00466 
00467 //=============================================================================
00468 int Epetra_SerialDenseSolver::Invert(void)
00469 {
00470   if (!Factored()) Factor(); // Need matrix factored.
00471 
00472   /* This section work with LAPACK Version 3.0 only 
00473   // Setting LWORK = -1 and calling GETRI will return optimal work space size in WORK_TMP
00474   int LWORK_TMP = -1;
00475   double WORK_TMP;
00476   GETRI ( N_, AF_, LDAF_, IPIV_, &WORK_TMP, &LWORK_TMP, &INFO_);
00477   LWORK_TMP = WORK_TMP; // Convert to integer
00478   if (LWORK_TMP>LWORK_) {
00479   if (WORK_!=0) delete [] WORK_;
00480   LWORK_ = LWORK_TMP;
00481   WORK_ = new double[LWORK_];
00482   }
00483   */
00484   // This section will work with any version of LAPACK 
00485   AllocateWORK();
00486 
00487   GETRI ( N_, AF_, LDAF_, IPIV_, WORK_, &LWORK_, &INFO_);
00488 
00489   double DN = N_;
00490   UpdateFlops((DN*DN*DN));
00491   Inverted_ = true;
00492   Factored_ = false;
00493   
00494   EPETRA_CHK_ERR(INFO_);
00495   return(0);
00496 }
00497 
00498 //=============================================================================
00499 int Epetra_SerialDenseSolver::ReciprocalConditionEstimate(double & Value)
00500 {
00501   int ierr = 0;
00502   if (ReciprocalConditionEstimated()) {
00503     Value = RCOND_;
00504     return(0); // Already computed, just return it.
00505   }
00506 
00507   if (ANORM_<0.0) ANORM_ = Matrix_->OneNorm();
00508   if (!Factored()) ierr = Factor(); // Need matrix factored.
00509   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00510 
00511   AllocateWORK();
00512   AllocateIWORK();
00513   // We will assume a one-norm condition number
00514   GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
00515   ReciprocalConditionEstimated_ = true;
00516   Value = RCOND_;
00517   UpdateFlops(2*N_*N_); // Not sure of count
00518   EPETRA_CHK_ERR(INFO_);
00519   return(0);
00520 }
00521 //=============================================================================
00522 void Epetra_SerialDenseSolver::Print(std::ostream& os) const {
00523 
00524   if (Matrix_!=0) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
00525   if (Factor_!=0) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
00526   if (LHS_   !=0) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
00527   if (RHS_   !=0) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;
00528 
00529 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines