Epetra_SerialDenseSVD.cpp

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

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