Epetra Package Browser (Single Doxygen Collection) Development
Epetra_SerialDenseSVD.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 
00044 #include "Epetra_SerialDenseSVD.h"
00045 #include "Epetra_SerialDenseMatrix.h"
00046 
00047 //=============================================================================
00048 Epetra_SerialDenseSVD::Epetra_SerialDenseSVD()
00049   : Epetra_Object("Epetra::SerialDenseSVD"),
00050     Epetra_CompObject(),
00051     Epetra_BLAS(),
00052     Epetra_LAPACK(),
00053     Transpose_(false),
00054     Factored_(false),
00055     Solved_(false),
00056     Inverted_(false),
00057     TRANS_('N'),
00058     M_(0),
00059    N_(0),
00060    Min_MN_(0),
00061    NRHS_(0),
00062    LDA_(0),
00063    LDAI_(0),
00064    LDB_(0),
00065    LDX_(0),
00066    INFO_(0),
00067    LWORK_(0),
00068    IWORK_(0),
00069    ANORM_(0.0),
00070   Matrix_(0),
00071   LHS_(0),
00072   RHS_(0),
00073   Inverse_(0),
00074   A_(0),
00075   AI_(0),
00076   WORK_(0),
00077   U_(0),
00078   S_(0),
00079   Vt_(0),
00080   B_(0),
00081   X_(0),
00082     UseTranspose_(false)
00083 {
00084   InitPointers();
00085   ResetMatrix();
00086   ResetVectors();
00087 }
00088 //=============================================================================
00089 Epetra_SerialDenseSVD::~Epetra_SerialDenseSVD()
00090 {
00091   DeleteArrays();
00092 }
00093 //=============================================================================
00094 void Epetra_SerialDenseSVD::InitPointers()
00095 {
00096   IWORK_ = 0;
00097 //  FERR_ = 0;
00098 //  BERR_ = 0;
00099 //  Factor_ =0;
00100   Inverse_ =0;
00101 //  AF_ = 0;
00102   AI_ = 0;
00103 //  IPIV_ = 0;
00104   WORK_ = 0;
00105 //  R_ = 0;
00106 //  C_ = 0;
00107   U_ = 0;
00108   S_ = 0;
00109   Vt_ = 0;
00110   INFO_ = 0;
00111   LWORK_ = 0;    
00112 }
00113 //=============================================================================
00114 void Epetra_SerialDenseSVD::DeleteArrays()
00115 {
00116   if (U_ !=0) {delete[] U_; U_ = 0;}
00117   if (S_ !=0) {delete[] S_; S_ = 0;}
00118   if (Vt_ !=0) {delete[] Vt_; Vt_ = 0;}
00119   if (IWORK_ != 0) {delete [] IWORK_; IWORK_ = 0;}
00120 //  if (FERR_ != 0)  {delete [] FERR_; FERR_ = 0;}
00121 //  if (BERR_ != 0)  {delete [] BERR_; BERR_ = 0;}
00122 //  if (Factor_ != Matrix_ && Factor_ != 0)   {delete Factor_; Factor_ = 0;}
00123 //  if (Factor_ !=0) Factor_ = 0;
00124   if (Inverse_ != 0)   {delete Inverse_; Inverse_ = 0;}
00125 //  if (AF_ !=0) AF_ = 0;
00126   if (AI_ !=0) AI_ = 0;
00127 //  if (IPIV_ != 0)  {delete [] IPIV_;IPIV_ = 0;}
00128   if (WORK_ != 0)  {delete [] WORK_;WORK_ = 0;}
00129 //  if (R_ != 0 && R_ != C_)     {delete [] R_;R_ = 0;}
00130 //  if (R_ != 0) R_ = 0;
00131 //  if (C_ != 0)     {delete [] C_;C_ = 0;}
00132   INFO_ = 0;
00133   LWORK_ = 0;    
00134 }
00135 //=============================================================================
00136 void Epetra_SerialDenseSVD::ResetMatrix()
00137 {
00138   DeleteArrays();
00139   ResetVectors();
00140   Matrix_ = 0;
00141   Inverse_ = 0;
00142 //  Factor_ = 0;
00143 //  A_Equilibrated_ = false;
00144   Factored_ = false;
00145   Inverted_ = false;
00146   M_ = 0;
00147   N_ = 0;
00148   Min_MN_ = 0;
00149   LDA_ = 0;
00150 //  LDAF_ = 0;
00151   LDAI_ = 0;
00152   ANORM_ = -1.0;
00153 //  RCOND_ = -1.0;
00154 //  ROWCND_ = -1.0;
00155 //  COLCND_ = -1.0;
00156 //  AMAX_ = -1.0;
00157   A_ = 0;
00158 
00159   if( U_ ) { delete [] U_; U_ = 0; }
00160   if( S_ ) { delete [] S_; S_ = 0; }
00161   if( Vt_ ) { delete [] Vt_; Vt_ = 0; }
00162 }
00163 //=============================================================================
00164 int Epetra_SerialDenseSVD::SetMatrix(Epetra_SerialDenseMatrix & A_in) {
00165   ResetMatrix();
00166   Matrix_ = &A_in;
00167 //  Factor_ = &A_in;
00168   M_ = A_in.M();
00169   N_ = A_in.N();
00170   Min_MN_ = EPETRA_MIN(M_,N_);
00171   LDA_ = A_in.LDA();
00172 //  LDAF_ = LDA_;
00173   A_ = A_in.A();
00174 //  AF_ = A_in.A();
00175   return(0);
00176 }
00177 //=============================================================================
00178 void Epetra_SerialDenseSVD::ResetVectors()
00179 {
00180   LHS_ = 0;
00181   RHS_ = 0;
00182   B_ = 0;
00183   X_ = 0;
00184 //  ReciprocalConditionEstimated_ = false;
00185 //  SolutionRefined_ = false;
00186   Solved_ = false;
00187 //  SolutionErrorsEstimated_ = false;
00188 //  B_Equilibrated_ = false;
00189   NRHS_ = 0;
00190   LDB_ = 0;
00191   LDX_ = 0;
00192 }
00193 //=============================================================================
00194 int Epetra_SerialDenseSVD::SetVectors(Epetra_SerialDenseMatrix & X_in, Epetra_SerialDenseMatrix & B_in)
00195 {
00196   if (B_in.M()!=X_in.M() || B_in.N() != X_in.N()) EPETRA_CHK_ERR(-1);
00197   if (B_in.A()==0) EPETRA_CHK_ERR(-2);
00198   if (B_in.LDA()<1) EPETRA_CHK_ERR(-3);
00199   if (X_in.A()==0) EPETRA_CHK_ERR(-4);
00200   if (X_in.LDA()<1) EPETRA_CHK_ERR(-5);
00201 
00202   ResetVectors(); 
00203   LHS_ = &X_in;
00204   RHS_ = &B_in;
00205   NRHS_ = B_in.N();
00206 
00207   B_ = B_in.A();
00208   LDB_ = B_in.LDA();
00209   X_ = X_in.A();
00210   LDX_ = X_in.LDA();
00211   return(0);
00212 }
00213 //=============================================================================
00214 int Epetra_SerialDenseSVD::Factor() {
00215 
00216   int ierr = 0;
00217 
00218   ANORM_ = Matrix_->OneNorm(); // Compute 1-Norm of A
00219 
00220   //allocate U_, S_, and Vt_ if not already done
00221   if(U_==0)
00222   {
00223     U_ = new double[M_*N_];
00224     S_ = new double[M_];
00225     Vt_ = new double[M_*N_];
00226   }
00227   else //zero them out
00228   {
00229     for( int i = 0; i < M_; ++i ) S_[i]=0.0;
00230     for( int i = 0; i < M_*N_; ++i )
00231     {
00232       U_[i]=0.0;
00233       Vt_[i]=0.0;
00234     }
00235   }
00236 
00237 //  if (Equilibrate_) ierr = EquilibrateMatrix();
00238 
00239   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00240   
00241   //allocate temp work space
00242   int lwork = 5*M_;
00243   double *work = new double[lwork];
00244   char job = 'A';
00245 
00246   //create temporary matrix to avoid writeover of original
00247   Epetra_SerialDenseMatrix tempMat( *Matrix_ );
00248   GESVD( job, job, M_, N_, tempMat.A(), LDA_, S_, U_, N_, Vt_, M_, work, &lwork, &INFO_ );
00249 
00250   delete [] work;
00251 
00252   Factored_ = true;
00253   double DN = N_;
00254   UpdateFlops(2.0*(DN*DN*DN)/3.0);
00255 
00256   EPETRA_CHK_ERR(INFO_);
00257   return(0);
00258 }
00259 
00260 //=============================================================================
00261 int Epetra_SerialDenseSVD::Solve(void) {
00262 
00263   //FOR NOW, ONLY ALLOW SOLVES IF INVERTED!!!!
00264   //NO REFINEMENT!!!
00265   //NO EQUILIBRATION!!!
00266 
00267   // We will call one of four routines depending on what services the user wants and 
00268   // whether or not the matrix has been inverted or factored already.
00269   //
00270   // If the matrix has been inverted, use DGEMM to compute solution.
00271   // Otherwise, if the user want the matrix to be equilibrated or wants a refined solution, we will
00272   // call the X interface.
00273   // Otherwise, if the matrix is already factored we will call the TRS interface.
00274   // Otherwise, if the matrix is unfactored we will call the SV interface.
00275 
00276 
00277 /*
00278   if (Equilibrate_) {
00279     ierr = EquilibrateRHS();
00280     B_Equilibrated_ = true;
00281   }
00282   EPETRA_CHK_ERR(ierr);
00283   if (A_Equilibrated_ && !B_Equilibrated_) EPETRA_CHK_ERR(-1); // Matrix and vectors must be similarly scaled
00284   if (!A_Equilibrated_ && B_Equilibrated_) EPETRA_CHK_ERR(-2);
00285   if (B_==0) EPETRA_CHK_ERR(-3); // No B
00286   if (X_==0) EPETRA_CHK_ERR(-4); // No B
00287 
00288   if (ShouldEquilibrate() && !A_Equilibrated_) ierr = 1; // Warn that the system should be equilibrated.
00289 */
00290 
00291   double DN = N_;
00292   double DNRHS = NRHS_;
00293   if (Inverted()) {
00294 
00295     if (B_==X_) EPETRA_CHK_ERR(-100); // B and X must be different for this case
00296 
00297     GEMM(TRANS_, 'N', N_, NRHS_, N_, 1.0, AI_, LDAI_, B_, LDB_, 0.0, X_, LDX_);
00298     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00299     UpdateFlops(2.0*DN*DN*DNRHS);
00300     Solved_ = true;
00301   }
00302   else EPETRA_CHK_ERR(-101); //Currently, for solve must have inverse
00303 /*
00304   else {
00305 
00306     if (!Factored()) Factor(); // Matrix must be factored
00307     
00308     if (B_!=X_) *LHS_ = *RHS_; // Copy B to X if needed
00309     GETRS(TRANS_, N_, NRHS_, AF_, LDAF_, IPIV_, X_, LDX_, &INFO_);
00310     if (INFO_!=0) EPETRA_CHK_ERR(INFO_);
00311     UpdateFlops(2.0*DN*DN*DNRHS);
00312     Solved_ = true;
00313 
00314   }
00315 
00316   int ierr1=0;
00317   if (RefineSolution_ && !Inverted()) ierr1 = ApplyRefinement();
00318   if (ierr1!=0) EPETRA_CHK_ERR(ierr1)
00319   else
00320     EPETRA_CHK_ERR(ierr);
00321   
00322   if (Equilibrate_) ierr1 = UnequilibrateLHS();
00323   EPETRA_CHK_ERR(ierr1);
00324 */
00325   return(0);
00326 }
00327 /*
00328 //=============================================================================
00329 int Epetra_SerialDenseSVD::ApplyRefinement(void)
00330 {
00331   double DN = N_;
00332   double DNRHS = NRHS_;
00333   if (!Solved()) EPETRA_CHK_ERR(-100); // Must have an existing solution
00334   if (A_==AF_) EPETRA_CHK_ERR(-101); // Cannot apply refine if no original copy of A.
00335 
00336   if (FERR_ != 0) delete [] FERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
00337   FERR_ = new double[NRHS_];
00338   if (BERR_ != 0) delete [] BERR_; // Always start with a fresh copy of FERR_, since NRHS_ may change
00339   BERR_ = new double[NRHS_];
00340   AllocateWORK();
00341   AllocateIWORK();
00342   
00343   GERFS(TRANS_, N_, NRHS_, A_, LDA_, AF_, LDAF_, IPIV_,
00344          B_, LDB_, X_, LDX_, FERR_, BERR_, 
00345          WORK_, IWORK_, &INFO_);
00346   
00347   
00348   SolutionErrorsEstimated_ = true;
00349   ReciprocalConditionEstimated_ = true;
00350   SolutionRefined_ = true;
00351   
00352   UpdateFlops(2.0*DN*DN*DNRHS); // Not sure of count
00353   
00354   EPETRA_CHK_ERR(INFO_);
00355   return(0);
00356 
00357 }
00358 
00359 //=============================================================================
00360 int Epetra_SerialDenseSVD::ComputeEquilibrateScaling(void) {
00361   if (R_!=0) return(0); // Already computed
00362  
00363   double DM = M_;
00364   double DN = N_;
00365   R_ = new double[M_];
00366   C_ = new double[N_];
00367   
00368   GEEQU (M_, N_, AF_, LDAF_, R_, C_, &ROWCND_, &COLCND_, &AMAX_, &INFO_);
00369   if (INFO_ != 0) EPETRA_CHK_ERR(INFO_);
00370 
00371   if (COLCND_<0.1 || ROWCND_<0.1 || AMAX_ < Epetra_Underflow || AMAX_ > Epetra_Overflow) ShouldEquilibrate_ = true;
00372 
00373   UpdateFlops(4.0*DM*DN);
00374   
00375   return(0);
00376 }
00377 
00378 //=============================================================================
00379 int Epetra_SerialDenseSVD::EquilibrateMatrix(void)
00380 {
00381   int i, j;
00382   int ierr = 0;
00383 
00384   double DN = N_;
00385   double DM = M_;
00386 
00387   if (A_Equilibrated_) return(0); // Already done
00388   if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
00389   if (ierr!=0) EPETRA_CHK_ERR(ierr);
00390   if (A_==AF_) {
00391     double * ptr;
00392     for (j=0; j<N_; j++) {
00393       ptr = A_ + j*LDA_;
00394       double s1 = C_[j];
00395       for (i=0; i<M_; i++) {
00396   *ptr = *ptr*s1*R_[i];
00397   ptr++;
00398       }
00399     }
00400     UpdateFlops(2.0*DM*DN);
00401   }
00402   else {
00403     double * ptr;
00404     double * ptr1;
00405     for (j=0; j<N_; j++) {
00406       ptr = A_ + j*LDA_;
00407       ptr1 = AF_ + j*LDAF_;
00408       double s1 = C_[j];
00409       for (i=0; i<M_; i++) {
00410   *ptr = *ptr*s1*R_[i];
00411   ptr++;
00412   *ptr1 = *ptr1*s1*R_[i];
00413   ptr1++;
00414       }
00415     }
00416     UpdateFlops(4.0*DM*DN);
00417   }
00418   
00419   A_Equilibrated_ = true;
00420   
00421   return(0);
00422 }
00423 
00424 //=============================================================================
00425 int Epetra_SerialDenseSVD::EquilibrateRHS(void)
00426 {
00427   int i, j;
00428   int ierr = 0;
00429 
00430   if (B_Equilibrated_) return(0); // Already done
00431   if (R_==0) ierr = ComputeEquilibrateScaling(); // Compute R and C if needed
00432   if (ierr!=0) EPETRA_CHK_ERR(ierr);
00433 
00434   double * R = R_;
00435   if (Transpose_) R = C_;
00436 
00437   double * ptr;
00438   for (j=0; j<NRHS_; j++) {
00439     ptr = B_ + j*LDB_;
00440     for (i=0; i<M_; i++) {
00441       *ptr = *ptr*R[i];
00442       ptr++;
00443     }
00444   }
00445 
00446   
00447   B_Equilibrated_ = true;
00448   UpdateFlops((double) N_*(double) NRHS_);
00449   
00450   return(0);
00451 }
00452 
00453 //=============================================================================
00454 int Epetra_SerialDenseSVD::UnequilibrateLHS(void)
00455 {
00456   int i, j;
00457 
00458   if (!B_Equilibrated_) return(0); // Nothing to do
00459 
00460   double * C = C_;
00461   if (Transpose_) C = R_;
00462 
00463   double * ptr;
00464   for (j=0; j<NRHS_; j++) {
00465     ptr = X_ + j*LDX_;
00466     for (i=0; i<N_; i++) {
00467       *ptr = *ptr*C[i];
00468       ptr++;
00469     }
00470   }
00471 
00472   
00473   UpdateFlops((double) N_ *(double) NRHS_);
00474   
00475   return(0);
00476 }
00477 */
00478 
00479 //=============================================================================
00480 int Epetra_SerialDenseSVD::Invert( double rthresh, double athresh )
00481 {
00482   if (!Factored()) Factor(); // Need matrix factored.
00483 
00484   //apply threshold
00485   double thresh = S_[0]*rthresh + athresh;
00486   int num_replaced = 0;
00487   for( int i = 0; i < M_; ++i )
00488     if( S_[i] < thresh )
00489     {
00490 //cout <<  num_replaced << thresh << " " << S_[0] << " " << S_[i] << endl;
00491 //      S_[i] = thresh;
00492       S_[i] = 0.0;
00493       ++num_replaced;
00494     }
00495 
00496   //scale cols of U_ with reciprocal singular values
00497   double *p = U_;
00498   for( int i = 0; i < N_; ++i )
00499   {
00500     double scale = 0.0;
00501     if( S_[i] ) scale = 1./S_[i];
00502     for( int j = 0; j < M_; ++j ) *p++ *= scale;
00503   }
00504 
00505   //create new Inverse_ if necessary
00506   if( Inverse_ == 0 )
00507   {
00508     Inverse_ = new Epetra_SerialDenseMatrix();
00509     Inverse_->Shape( N_, M_ );
00510     AI_ = Inverse_->A();
00511     LDAI_ = Inverse_->LDA();
00512   }
00513 /*
00514   else //zero it out
00515   {
00516     for( int i = 0; i < Inverse_->M(); ++i )
00517       for( int j = 0; j < Inverse_->N(); ++j )
00518         (*Inverse_)(i,j) = 0.0;
00519   }
00520 */
00521 
00522   GEMM( 'T', 'T', M_, M_, M_, 1.0, Vt_, M_, U_, M_, 0.0, AI_, M_ );
00523 
00524   double DN = N_;
00525   UpdateFlops((DN*DN*DN));
00526   Inverted_ = true;
00527   Factored_ = false;
00528   
00529   EPETRA_CHK_ERR(INFO_);
00530   return(num_replaced);
00531 }
00532 
00533 /*
00534 //=============================================================================
00535 int Epetra_SerialDenseSVD::ReciprocalConditionEstimate(double & Value)
00536 {
00537   int ierr = 0;
00538   if (ReciprocalConditionEstimated()) {
00539     Value = RCOND_;
00540     return(0); // Already computed, just return it.
00541   }
00542 
00543   if (ANORM_<0.0) ANORM_ = Matrix_->OneNorm();
00544   if (!Factored()) ierr = Factor(); // Need matrix factored.
00545   if (ierr!=0) EPETRA_CHK_ERR(ierr-2);
00546 
00547   AllocateWORK();
00548   AllocateIWORK();
00549   // We will assume a one-norm condition number
00550   GECON( '1', N_, AF_, LDAF_, ANORM_, &RCOND_, WORK_, IWORK_, &INFO_);
00551   ReciprocalConditionEstimated_ = true;
00552   Value = RCOND_;
00553   UpdateFlops(2*N_*N_); // Not sure of count
00554   EPETRA_CHK_ERR(INFO_);
00555   return(0);
00556 }
00557 */
00558 
00559 //=============================================================================
00560 void Epetra_SerialDenseSVD::Print(ostream& os) const {
00561 
00562   if (Matrix_!=0) os << *Matrix_;
00563 //  if (Factor_!=0) os << *Factor_;
00564   if (S_!=0) for( int i = 0; i < M_; ++i ) cout << "(" << i << "," << S_[i] << ")\n";
00565   if (Inverse_!=0) os << *Inverse_;
00566   if (LHS_!=0) os << *LHS_;
00567   if (RHS_!=0) os << *RHS_;
00568 
00569 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines