Epetra_SerialDenseSolver.cpp

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

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7