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

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