Epetra Package Browser (Single Doxygen Collection) Development
Epetra_SerialSpdDenseSolver.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_SerialSpdDenseSolver.h"
00044 #include "Epetra_SerialSymDenseMatrix.h"
00045 #include "Epetra_SerialDenseMatrix.h"
00046 
00047 //=============================================================================
00048 Epetra_SerialSpdDenseSolver::Epetra_SerialSpdDenseSolver(void)
00049   : Epetra_SerialDenseSolver(),
00050     SCOND_(-1.0),
00051     SymMatrix_(0),
00052     SymFactor_(0)
00053 {
00054 }
00055 //=============================================================================
00056 Epetra_SerialSpdDenseSolver::~Epetra_SerialSpdDenseSolver()
00057 {
00058   if (SymFactor_ != SymMatrix_ && SymFactor_ != 0) {
00059     delete SymFactor_; SymFactor_ = 0; Factor_ = 0;
00060   }
00061 }
00062 //=============================================================================
00063 int Epetra_SerialSpdDenseSolver::SetMatrix(Epetra_SerialSymDenseMatrix & A_in) {
00064   
00065   SymMatrix_=&A_in;
00066   SymFactor_=&A_in;
00067   SCOND_ = -1.0;
00068   // Also call SerialDensematrix set method
00069   return(Epetra_SerialDenseSolver::SetMatrix( (Epetra_SerialDenseMatrix &) A_in));
00070 }
00071 //=============================================================================
00072 int Epetra_SerialSpdDenseSolver::Factor(void) {
00073   if (Factored()) return(0); // Already factored
00074   if (Inverted()) EPETRA_CHK_ERR(-100); // Cannot factor inverted matrix
00075   int ierr = 0;
00076 
00077   ANORM_ = SymMatrix_->OneNorm();
00078 
00079 
00080   // If we want to refine the solution, then the factor must
00081   // be stored separatedly from the original matrix
00082 
00083   if (A_ == AF_)
00084     if (RefineSolution_ ) {
00085       SymFactor_ = new Epetra_SerialSymDenseMatrix(*SymMatrix_);
00086       Factor_ = SymFactor_;
00087       AF_ = SymFactor_->A();
00088       LDAF_ = SymFactor_->LDA();
00089     }
00090   if (Equilibrate_) ierr = EquilibrateMatrix();
00091 
00092   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00093   
00094   POTRF (SymMatrix_->UPLO(), N_, AF_, LDAF_, &INFO_);
00095   Factored_ = true;
00096   double DN = N_;
00097   UpdateFlops((DN*DN*DN)/3.0);
00098 
00099   EPETRA_CHK_ERR(INFO_);
00100   return(0);
00101 
00102 }
00103 
00104 //=============================================================================
00105 int Epetra_SerialSpdDenseSolver::Solve(void) {
00106   int ierr = 0;
00107 
00108   // We will call one of four routines depending on what services the user wants and 
00109   // whether or not the matrix has been inverted or factored already.
00110   //
00111   // If the matrix has been inverted, use DGEMM to compute solution.
00112   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00113   // call the X interface.
00114   // Otherwise, if the matrix is already factored we will call the TRS interface.
00115   // Otherwise, if the matrix is unfactored we will call the SV interface.
00116 
00117   if (Equilibrate_) {
00118     ierr = Epetra_SerialDenseSolver::EquilibrateRHS();
00119     B_Equilibrated_ = true;
00120   }
00121   EPETRA_CHK_ERR(ierr);
00122   if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
00123   if (!A_Equilibrated_ && B_Equilibrated_) EPETRA_CHK_ERR(-2);
00124   if (B_==0) EPETRA_CHK_ERR(-3); // No B
00125   if (X_==0) EPETRA_CHK_ERR(-4); // No B
00126 
00127   if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
00128 
00129   double DN = N_;
00130   double DNRHS = NRHS_;
00131   if (Inverted()) {
00132 
00133     if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
00134     GEMM('N', 'N', N_, NRHS_, N_, 1.0, AF_, LDAF_,
00135     B_, LDB_, 0.0, X_, LDX_);
00136     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00137     UpdateFlops(2.0*DN*DN*DNRHS);
00138     Solved_ = true;
00139   }
00140   else {
00141 
00142     if (!Factored()) Factor(); // Matrix must be factored
00143     if (B_!=X_) {
00144        *LHS_ = *RHS_; // Copy B to X if needed 
00145        X_ = LHS_->A(); LDX_ = LHS_->LDA();
00146     }
00147     
00148     POTRS(SymMatrix_->UPLO(), N_, NRHS_, AF_, LDAF_, X_, LDX_, &INFO_);
00149     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00150     UpdateFlops(2.0*DN*DN*DNRHS);
00151     Solved_ = true;
00152 
00153   }
00154   int ierr1=0;
00155   if (RefineSolution_) ierr1 = ApplyRefinement();
00156   if (ierr1!=0) {
00157     EPETRA_CHK_ERR(ierr1);
00158   }
00159   else {
00160     EPETRA_CHK_ERR(ierr);
00161   }
00162   if (Equilibrate_) ierr1 = Epetra_SerialDenseSolver::UnequilibrateLHS();
00163   EPETRA_CHK_ERR(ierr1);
00164   return(0);
00165 }
00166 //=============================================================================
00167 int Epetra_SerialSpdDenseSolver::ApplyRefinement(void)
00168 {
00169   double DN = N_;
00170   double DNRHS = NRHS_;
00171   if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
00172   if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
00173 
00174   if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
00175   FERR_ = new double[NRHS_];
00176   if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
00177   BERR_ = new double[NRHS_];
00178   AllocateWORK();
00179   AllocateIWORK();
00180   
00181   PORFS(SymMatrix_->UPLO(), N_, NRHS_, A_, LDA_, AF_, LDAF_,
00182          B_, LDB_, X_, LDX_, FERR_, BERR_, 
00183          WORK_, IWORK_, &INFO_);
00184   
00185   
00186   SolutionErrorsEstimated_ = true;
00187   ReciprocalConditionEstimated_ = true;
00188   SolutionRefined_ = true;
00189   
00190   UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
00191   
00192   EPETRA_CHK_ERR(INFO_);
00193   return(0);
00194   
00195 }
00196 
00197 //=============================================================================
00198 int Epetra_SerialSpdDenseSolver::ComputeEquilibrateScaling(void) {
00199   if (R_!=0) return(0); // Already computed
00200  
00201   double DN = N_;
00202   R_ = new double[N_];
00203   C_ = R_;
00204   
00205   POEQU (N_, AF_, LDAF_, R_, &SCOND_, &AMAX_, &INFO_);
00206   if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
00207 
00208   if (SCOND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
00209 
00210   C_ = R_; // Set column scaling pointer so we can use EquilibrateRHS and UnequilibrateLHS from base class
00211   UpdateFlops(2.0*DN*DN);
00212   
00213   return(0);
00214 }
00215 
00216 //=============================================================================
00217 int Epetra_SerialSpdDenseSolver::EquilibrateMatrix(void) {
00218   int i, j;
00219   int ierr = 0;
00220 
00221   if (A_Equilibrated_) return(0); // Already done
00222   if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute S if needed
00223   if (ierr!=0) EPETRA_CHK_ERR(ierr);
00224   if (SymMatrix_->Upper()) {
00225     if (A_==AF_) {
00226       double * ptr;
00227       for (j=0; j<N_; j++) {
00228   ptr = A_ + j*LDA_;
00229   double s1 = R_[j];
00230   for (i=0; i<=j; i++) {
00231     *ptr = *ptr*s1*R_[i];
00232     ptr++;
00233   }
00234       }
00235     }
00236     else {
00237       double * ptr;
00238       double * ptr1;
00239       for (j=0; j<N_; j++) {
00240   ptr = A_ + j*LDA_;
00241   ptr1 = AF_ + j*LDAF_;
00242   double s1 = R_[j];
00243   for (i=0; i<=j; i++) {
00244     *ptr = *ptr*s1*R_[i];
00245     ptr++;
00246     *ptr1 = *ptr1*s1*R_[i];
00247     ptr1++;
00248   }
00249       }
00250     }
00251   }
00252   else {
00253     if (A_==AF_) {
00254       double * ptr;
00255       for (j=0; j<N_; j++) {
00256   ptr = A_ + j + j*LDA_;
00257   double s1 = R_[j];
00258   for (i=j; i<N_; i++) {
00259     *ptr = *ptr*s1*R_[i];
00260     ptr++;
00261   }
00262       }
00263     }
00264     else {
00265       double * ptr;
00266       double * ptr1;
00267       for (j=0; j<N_; j++) {
00268   ptr = A_ + j + j*LDA_;
00269   ptr1 = AF_ + j + j*LDAF_;
00270   double s1 = R_[j];
00271   for (i=j; i<N_; i++) {
00272     *ptr = *ptr*s1*R_[i];
00273     ptr++;
00274     *ptr1 = *ptr1*s1*R_[i];
00275     ptr1++;
00276   }
00277       }
00278     }
00279   }
00280   A_Equilibrated_ = true;
00281   double NumFlops = (double) ((N_+1)*N_/2);
00282   if (A_==AF_) NumFlops += NumFlops;
00283   UpdateFlops(NumFlops);
00284   
00285   return(0);
00286 }
00287 
00288 //=============================================================================
00289 int Epetra_SerialSpdDenseSolver::Invert(void)
00290 {
00291   if (!Factored()) Factor(); // Need matrix factored.
00292   POTRI ( SymMatrix_->UPLO(), N_, AF_, LDAF_, &INFO_);
00293   // Copy lower/upper triangle to upper/lower triangle: make full inverse
00294   SymMatrix_->CopyUPLOMat(SymMatrix_->Upper(), AF_, LDAF_, N_);
00295   double DN = N_;
00296   UpdateFlops((DN*DN*DN));
00297   Inverted_ = true;
00298   Factored_ = false;
00299   
00300   EPETRA_CHK_ERR(INFO_);
00301   return(0);
00302 }
00303 
00304 //=============================================================================
00305 int Epetra_SerialSpdDenseSolver::ReciprocalConditionEstimate(double & Value)
00306 {
00307   int ierr = 0;
00308   if (ReciprocalConditionEstimated()) {
00309     Value = RCOND_;
00310     return(0); // Already computed, just return it.
00311   }
00312 
00313   if (ANORM_<0.0) ANORM_ = SymMatrix_->OneNorm();
00314   if (ierr!=0) EPETRA_CHK_ERR(ierr-1);
00315   if (!Factored()) ierr = Factor(); // Need matrix factored.
00316   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00317 
00318   AllocateWORK();
00319   AllocateIWORK();
00320   POCON( SymMatrix_->UPLO(), N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
00321   ReciprocalConditionEstimated_ = true;
00322   Value = RCOND_;
00323   UpdateFlops(2*N_*N_); // Not sure of count
00324   EPETRA_CHK_ERR(INFO_);
00325   return(0);
00326 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines