AbstractLinAlgPack_DirectSparseSolverMA28.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include <fstream>
00032 #include <algorithm>
00033 
00034 #include "AbstractLinAlgPack_DirectSparseSolverMA28.hpp"
00035 #include "AbstractLinAlgPack_MatrixScaling_Strategy.hpp"
00036 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00037 #include "DenseLinAlgPack_PermVecMat.hpp"
00038 #include "Teuchos_AbstractFactoryStd.hpp"
00039 #include "Teuchos_TestForException.hpp"
00040 #include "Teuchos_Workspace.hpp"
00041 #include "Teuchos_dyn_cast.hpp"
00042 #include "FortranTypes_f_open_file.hpp"
00043 
00044 namespace {
00045 //
00046 template< class T >
00047 inline
00048 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
00049 // A cast to const is needed because the standard does not return a reference from
00050 // valarray<>::operator[]() const.
00051 template <class T>
00052 std::valarray<T>& cva(const std::valarray<T>& va )
00053 {
00054   return const_cast<std::valarray<T>&>(va);
00055 }
00056 }
00057 
00058 namespace AbstractLinAlgPack {
00059 
00060 //
00061 // Implementation of DirectSparseSolver(Imp) interface using MA28.
00062 //
00063 // Here are some little bits of knowledge about MA28 that I need
00064 // to record after many hours of hard work.
00065 //
00066 // * The 1979 paper in ACM TOMS (Vol. 5, No. 1, pages 27), seems
00067 // to suggest that MA28 pivots by column for numerical stability
00068 // but I am not sure about this.
00069 //
00070 // * When factoring a rectangular matrix, you must set 
00071 // LBLOCK = .FALSE. or the row and column permutations
00072 // extracted from IKEEP(:,2) and IKEEP(:,3) respectivly
00073 // are meaningless.
00074 //
00075 // ToDo: Finish this discussion!
00076 //
00077 
00078 // ToDo:
00079 // a) Add an option for printing the values of the common
00080 //    block parameters to out or to a file.  This can
00081 //    be used to get a feel for the performance of
00082 //    ma28
00083 // b) Add provisions for an external client to change
00084 //    the control options of MA28.  Most of these are
00085 //    stored as common block variables.
00086 
00087 // //////////////////////////////////////////////////
00088 // DirectSparseSolverMA28::FactorizationStructureMA28
00089 
00090 DirectSparseSolverMA28::FactorizationStructureMA28::FactorizationStructureMA28()
00091   :m_(0),n_(0),max_n_(0),nz_(0),licn_(0),lirn_(0)
00092    ,u_(0.1),scaling_(NO_SCALING)
00093 {}
00094 
00095 // //////////////////////////////////////////////////
00096 // DirectSparseSolverMA28::BasisMatrixMA28
00097 
00098 // Overridden from BasisMatrixImp
00099 
00100 Teuchos::RefCountPtr<DirectSparseSolverImp::BasisMatrixImp>
00101 DirectSparseSolverMA28::BasisMatrixMA28::create_matrix() const
00102 {
00103   return Teuchos::rcp(new BasisMatrixMA28);
00104 }
00105 
00106 void DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(
00107   VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
00108   ) const 
00109 {
00110   using Teuchos::dyn_cast;
00111   using Teuchos::Workspace;
00112   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00113   size_type k;
00114 
00115   // Get concrete objects
00116   const FactorizationStructureMA28
00117     &fs = dyn_cast<const FactorizationStructureMA28>(*this->get_fact_struc());
00118   const FactorizationNonzerosMA28
00119     &fn = dyn_cast<const FactorizationNonzerosMA28>(*this->get_fact_nonzeros());
00120 
00121   // Validate input
00122 #ifdef TEUCHOS_DEBUG
00123   TEST_FOR_EXCEPTION(
00124     y == NULL, std::invalid_argument
00125     ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
00126 #endif
00127   const size_type y_dim = y->dim(), x_dim = x.dim();
00128 #ifdef TEUCHOS_DEBUG
00129   TEST_FOR_EXCEPTION(
00130     fs.rank_ != y_dim, std::invalid_argument
00131     ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
00132   TEST_FOR_EXCEPTION(
00133     fs.rank_ != x_dim, std::invalid_argument
00134     ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
00135 #endif
00136 
00137   VectorDenseMutableEncap  yd(*y);
00138   VectorDenseEncap         xd(x);
00139 
00140   // Allocate workspace memory
00141   Workspace<value_type>  xfull_s(wss,fs.max_n_,false);
00142   DVectorSlice                 xfull(&xfull_s[0],xfull_s.size());
00143   Workspace<value_type>  ws(wss,fs.max_n_,false);
00144   DVectorSlice                 w(&ws[0],ws.size());
00145 
00146   // Get a context for transpose or no transpose
00147   const IVector
00148     &row_perm = (M_trans == BLAS_Cpp::no_trans) ? fs.row_perm_ : fs.col_perm_,
00149     &col_perm = (M_trans == BLAS_Cpp::no_trans) ? fs.col_perm_ : fs.row_perm_;
00150 
00151   // Copy x into positions in full w
00152   // Here, the rhs vector is set with only those equations that
00153   // are part of the nonsingular set.  It is important that the
00154   // ordering be the same as the original ordering sent to
00155   // MA28AD().
00156   xfull = 0.0;
00157   for( k = 1; k <= x_dim; ++k ) 
00158     xfull(row_perm(k)) = xd()(k);
00159   
00160   // Scale the rhs
00161   if( fs.matrix_scaling_.get() )
00162     fs.matrix_scaling_->scale_rhs( M_trans, xfull.raw_ptr() );
00163 
00164   // Solve for the rhs
00165   FortranTypes::f_int mtype = ( (M_trans == BLAS_Cpp::no_trans) ? 1 : 0 );
00166   fs.ma28_.ma28cd(
00167     fs.max_n_, &cva(fn.a_)[0], fs.licn_, &cva(fs.icn_)[0], &cva(fs.ikeep_)[0]
00168     ,xfull.raw_ptr(), w.raw_ptr(), mtype );
00169 
00170   // Scale the lhs
00171   if( fs.matrix_scaling_.get() )
00172     fs.matrix_scaling_->scale_rhs( M_trans, xfull.raw_ptr() );
00173 
00174   // Copy the solution into y
00175   // Here, the solution vector is set with only those variables that
00176   // are in the basis.  It is important that the
00177   // ordering be the same as the original ordering sent to
00178   // MA28AD().
00179   for( k = 1; k <= y_dim; ++k )
00180     yd()(k) = xfull(col_perm(k));
00181   
00182 }
00183 
00184 // //////////////////////////////////////////////////
00185 // DirectSparseSolverMA28
00186 
00187 // Constructors/initializers
00188 
00189 DirectSparseSolverMA28::DirectSparseSolverMA28(
00190   value_type          estimated_fillin_ratio
00191   ,value_type         u
00192   ,bool               grow
00193   ,value_type         tol
00194   ,index_type         nsrch
00195   ,bool               lbig
00196   ,bool               print_ma28_outputs
00197   ,const std::string& output_file_name
00198   )
00199   :estimated_fillin_ratio_(estimated_fillin_ratio)
00200   ,u_(u)
00201   ,grow_(grow)
00202   ,tol_(tol)
00203   ,nsrch_(nsrch)
00204   ,lbig_(lbig)
00205   ,print_ma28_outputs_(print_ma28_outputs)
00206   ,output_file_name_(output_file_name)
00207   ,file_output_num_(0)
00208 {}
00209 
00210 // Overridden from DirectSparseSolver
00211 
00212 const DirectSparseSolver::basis_matrix_factory_ptr_t
00213 DirectSparseSolverMA28::basis_matrix_factory() const
00214 {
00215   namespace mmp = MemMngPack;
00216   return Teuchos::rcp(new Teuchos::AbstractFactoryStd<BasisMatrix,BasisMatrixMA28>());
00217 }
00218 
00219 void DirectSparseSolverMA28::estimated_fillin_ratio(
00220   value_type estimated_fillin_ratio
00221   )
00222 {
00223   estimated_fillin_ratio_ = estimated_fillin_ratio;
00224 }
00225 
00226 // Overridden from DirectSparseSolverImp
00227 
00228 const Teuchos::RefCountPtr<DirectSparseSolver::FactorizationStructure>
00229 DirectSparseSolverMA28::create_fact_struc() const
00230 {
00231   return Teuchos::rcp(new FactorizationStructureMA28);
00232 }
00233 
00234 const Teuchos::RefCountPtr<DirectSparseSolverImp::FactorizationNonzeros>
00235 DirectSparseSolverMA28::create_fact_nonzeros() const
00236 {
00237   return Teuchos::rcp(new FactorizationNonzerosMA28);
00238 }
00239 
00240 void DirectSparseSolverMA28::imp_analyze_and_factor(
00241   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00242   ,FactorizationStructure                         *fact_struc
00243   ,FactorizationNonzeros                          *fact_nonzeros
00244   ,DenseLinAlgPack::IVector                            *row_perm
00245   ,DenseLinAlgPack::IVector                            *col_perm
00246   ,size_type                                      *rank
00247   ,std::ostream                                   *out
00248   )
00249 {
00250   using Teuchos::dyn_cast;
00251   typedef MatrixConvertToSparse MCTS;
00252   using Teuchos::Workspace;
00253   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00254 
00255   if(out)
00256     *out << "\nUsing MA28 to analyze and factor a new matrix ...\n";
00257 
00258   // Get the concrete factorization and nonzeros objects
00259   FactorizationStructureMA28
00260     &fs = dyn_cast<FactorizationStructureMA28>(*fact_struc);
00261   FactorizationNonzerosMA28
00262     &fn = dyn_cast<FactorizationNonzerosMA28>(*fact_nonzeros);
00263 
00264   // Set MA28 parameters
00265   set_ma28_parameters(&fs);
00266   
00267   // Get the dimensions of things.
00268   const index_type
00269     m = A.rows(),
00270     n = A.cols(),
00271     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00272 
00273   // Validate input
00274   TEST_FOR_EXCEPTION(
00275     n <= 0 || m <= 0 || m > n, std::invalid_argument
00276     ,"DirectSparseSolverMA28::imp_analyze_and_factor(...) : Error!" );
00277 
00278   // Memorize the dimenstions for checks later
00279   fs.m_ = m; fs.n_ = n; fs.nz_ = nz;
00280   fs.max_n_ = my_max(fs.m_,fs.n_);
00281 
00282   // By default set licn and ircn equal to estimated_fillin_ratio * nz.
00283   if( estimated_fillin_ratio_ < 1.0 ) {
00284     if( out ) *out << "Warning, client set estimated_fillin_ratio = " << estimated_fillin_ratio_
00285              << " < 1.0.\nSetting estimated_fillin_ratio = 10.0 ...\n";
00286     estimated_fillin_ratio_ = 10.0;
00287   }
00288   if( fs.licn_ < fs.nz_ ) fs.licn_ = static_cast<index_type>(estimated_fillin_ratio_ * fs.nz_);
00289   if( fs.lirn_ < fs.nz_ ) fs.lirn_ = static_cast<index_type>(estimated_fillin_ratio_ * fs.nz_);
00290 
00291   // Initialize structure storage
00292   fs.ivect_.resize(fs.nz_); // perminatly stores nz row indexes
00293   fs.jvect_.resize(fs.nz_); // perminatly stores nz column indexes
00294 
00295   index_type iflag = 0;
00296   for( int num_fac = 0; num_fac < 5; ++num_fac ) {
00297     
00298     // Initialize matrix factorization storage and temporary storage
00299     fs.icn_.resize(fs.licn_); // first nz entries stores the column indexes
00300      fn.a_.resize(fs.licn_);
00301     fs.ikeep_.resize( fs.ma28_.lblock() ? 5*fs.max_n_ :  4*fs.max_n_ + 1 );
00302     Workspace<index_type>  irn_tmp_(wss,fs.lirn_), iw(wss,8*fs.max_n_);
00303     Workspace<value_type>  w(wss,fs.max_n_);
00304     
00305     // Fill in the matrix information
00306     A.coor_extract_nonzeros(
00307       MCTS::EXTRACT_FULL_MATRIX
00308       ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00309       ,fs.nz_
00310       ,&fn.a_[0]
00311       ,fs.nz_
00312       ,&fs.ivect_[0]
00313       ,&fs.jvect_[0]
00314       );
00315     std::copy( &fs.ivect_[0], &fs.ivect_[0] + fs.nz_, &irn_tmp_[0] );
00316     std::copy( &fs.jvect_[0], &fs.jvect_[0] + fs.nz_, &fs.icn_[0] );
00317     
00318     // Scale the matrix
00319     if( fs.matrix_scaling_.get() )
00320       fs.matrix_scaling_->scale_matrix(
00321         fs.m_, fs.n_, fs.nz_, &fs.ivect_[0], &fs.jvect_[0], true
00322         ,&fn.a_[0]
00323         );
00324     
00325     // Analyze and factor the matrix
00326     if(out)
00327       *out << "\nCalling ma28ad(...) ...\n";
00328     fs.ma28_.ma28ad(
00329       fs.max_n_, fs.nz_, &fn.a_[0], fs.licn_, &irn_tmp_[0], fs.lirn_, &fs.icn_[0], fs.u_
00330       ,&fs.ikeep_[0], &iw[0], &w[0], &iflag
00331       );
00332     
00333     if(iflag != 0 && out)
00334       *out << "\nMA28AD returned iflag = " << iflag << " != 0!\n";
00335 
00336     // Print MA28 outputs
00337     print_ma28_outputs(true,iflag,fs,&w[0],out);
00338 
00339     if( iflag >= 0 ) break;
00340 
00341     switch( iflag ) {
00342       case LICN_AND_LIRN_TOO_SMALL:
00343       case LICN_TOO_SMALL:
00344       case LICN_FAR_TOO_SMALL:
00345       case LIRN_TOO_SMALL:
00346         if(out)
00347           *out << "\nWarning! iflag = " << iflag << ", LIRN and/or LICN are too small!\n"
00348              << "Increasing lirn = " << fs.lirn_ << " amd licn = " << fs.licn_
00349              << " by a factor of 10\n"
00350              << "(try increasing estimated_fillin_ratio = " << estimated_fillin_ratio_
00351              << " to a larger value next time)...\n";
00352         fs.lirn_ = 10 * fs.lirn_;
00353         fs.licn_ = 10 * fs.licn_;
00354         break;
00355     }
00356   }
00357 
00358   // Check for errors and throw exception if you have to.
00359   ThrowIFlagException(iflag);
00360 
00361   // Extract the basis matrix selection
00362 
00363   *rank = fs.ma28_.irank();
00364 
00365   row_perm->resize(fs.m_);
00366   if( *rank < fs.m_ ) {
00367     index_type
00368       *row_perm_ikeep = &fs.ikeep_[fs.max_n_],
00369       *row_perm_itr   = &(*row_perm)(1),
00370       *row_perm_last  = row_perm_itr + fs.m_;
00371     for(; row_perm_itr != row_perm_last;)
00372       *row_perm_itr++ = abs(*row_perm_ikeep++);
00373     // Sort partitions in assending order (required!)
00374     std::sort(&(*row_perm)[0]           , &(*row_perm)[0] + (*rank) );
00375     std::sort(&(*row_perm)[0] + (*rank) , &(*row_perm)[0] + m       );
00376   }
00377   else {
00378     DenseLinAlgPack::identity_perm( row_perm );
00379   }
00380 
00381   col_perm->resize(fs.n_);
00382   if( *rank < fs.n_ ) {
00383     index_type
00384       *col_perm_ikeep = &fs.ikeep_[2*fs.max_n_],
00385       *col_perm_itr   = &(*col_perm)(1),
00386       *col_perm_last  = col_perm_itr + fs.n_;
00387     for(; col_perm_itr != col_perm_last;)
00388       *col_perm_itr++ = abs(*col_perm_ikeep++);
00389     // Sort partitions in assending order (required!)
00390     std::sort(&(*col_perm)[0]           , &(*col_perm)[0] + (*rank) );
00391     std::sort(&(*col_perm)[0] + (*rank) , &(*col_perm)[0] + n       );
00392   }
00393   else {
00394     DenseLinAlgPack::identity_perm( col_perm );
00395   }
00396 
00397   // Set internal copy of basis selection
00398   fs.row_perm_ = *row_perm;
00399   fs.col_perm_ = *col_perm;
00400   fs.rank_     = *rank;
00401 
00402 }
00403 
00404 void DirectSparseSolverMA28::imp_factor(
00405   const AbstractLinAlgPack::MatrixConvertToSparse   &A
00406   ,const FactorizationStructure                   &fact_struc
00407   ,FactorizationNonzeros                          *fact_nonzeros
00408   ,std::ostream                                   *out
00409   )
00410 {
00411   using Teuchos::dyn_cast;
00412   typedef MatrixConvertToSparse MCTS;
00413   using Teuchos::Workspace;
00414   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00415 
00416   if(out)
00417     *out << "\nUsing MA28 to factor a new matrix ...\n";
00418 
00419   // Get the concrete factorization and nonzeros objects
00420   const FactorizationStructureMA28
00421     &fs = dyn_cast<const FactorizationStructureMA28>(fact_struc);
00422   FactorizationNonzerosMA28
00423     &fn = dyn_cast<FactorizationNonzerosMA28>(*fact_nonzeros);
00424 
00425   // Get the dimensions of things.
00426   const index_type
00427     m = A.rows(),
00428     n = A.cols(),
00429     nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
00430 
00431   // Validate input
00432 #ifdef TEUCHOS_DEBUG
00433   TEST_FOR_EXCEPTION(
00434     m != fs.m_ || n != fs.n_ || fs.nz_ != nz, std::invalid_argument
00435     ,"DirectSparseSolverMA28::imp_factor(...) : Error, "
00436     "A is not compatible with matrix passed to imp_analyze_and_factor()!" );
00437 #endif
00438 
00439   // Initialize matrix factorization storage and temporary storage
00440   if(fn.a_.size() < fs.licn_)  fn.a_.resize(fs.licn_);
00441   Workspace<index_type>   iw(wss,5*fs.max_n_);
00442   Workspace<value_type>   w(wss,fs.max_n_);
00443 
00444   // Fill in the matrix nonzeros (we already have the structure)
00445   A.coor_extract_nonzeros(
00446     MCTS::EXTRACT_FULL_MATRIX
00447     ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
00448     ,fs.nz_
00449     ,&fn.a_[0]
00450     ,0
00451     ,NULL
00452     ,NULL
00453     );
00454 
00455   // Scale the matrix
00456   if( fs.matrix_scaling_.get() )
00457     fs.matrix_scaling_->scale_matrix(
00458       fs.m_, fs.n_, fs.nz_, &cva(fs.ivect_)[0], &cva(fs.jvect_)[0], false
00459       ,&fn.a_[0]
00460       );
00461 
00462   // Factor the matrix
00463   index_type iflag = 0;
00464   fs.ma28_.ma28bd(
00465     fs.max_n_, fs.nz_, &fn.a_[0], fs.licn_, &cva(fs.ivect_)[0], &cva(fs.jvect_)[0], &cva(fs.icn_)[0]
00466     ,&cva(fs.ikeep_)[0], &iw[0], &w[0], &iflag
00467     );
00468 
00469   if(iflag != 0 && out)
00470     *out << "\nMA28BD returned iflag = " << iflag << " != 0!\n";
00471 
00472   // Print MA28 outputs
00473   print_ma28_outputs(false,iflag,fs,&w[0],out);
00474 
00475   // Check for errors and throw exception if you have to.
00476   ThrowIFlagException(iflag);
00477 
00478 }
00479 
00480 // private
00481 
00482 void DirectSparseSolverMA28::set_ma28_parameters( FactorizationStructureMA28* fs )
00483 {
00484   // Set common block parameters
00485   fs->ma28_.lblock( FortranTypes::F_FALSE ); // Do not permute to block triangular form (*** This is critical!)
00486   fs->u_ = u_;
00487   fs->ma28_.grow( grow_ ? FortranTypes::F_TRUE : FortranTypes::F_FALSE );
00488   fs->ma28_.tol(tol_);
00489   fs->ma28_.nsrch(nsrch_);
00490   fs->ma28_.lbig( lbig_ ? FortranTypes::F_TRUE : FortranTypes::F_FALSE );
00491   // Setup output file
00492   if( output_file_name_.length() > 0 && fs->ma28_.lp() == 0 ) {
00493     // Open a fortran file
00494     index_type iout = 25; // Unique?
00495     FortranTypes::f_open_file( iout, output_file_name_.c_str() );
00496     fs->ma28_.mp(iout);
00497     fs->ma28_.lp(iout);
00498   }
00499   else if( output_file_name_.length() == 0 && fs->ma28_.lp() ) {
00500     fs->ma28_.mp(0);
00501     fs->ma28_.lp(0);
00502   } 
00503 }
00504 
00505 void DirectSparseSolverMA28::print_ma28_outputs(
00506   bool                               ma28ad_bd
00507   ,index_type                        iflag
00508   ,const FactorizationStructureMA28  &fs
00509   ,const value_type                  w[]
00510   ,std::ostream                      *out
00511   )
00512 {
00513   if( print_ma28_outputs_ && out ) {
00514     *out << "\nReturn parameters from MA28 (call number = " << ++file_output_num_ << ")\n";
00515     if( fs.ma28_.grow() == FortranTypes::F_TRUE )
00516       *out << "w(1)   = " << w[0] << std::endl;
00517     *out << "rmin   = " << fs.ma28_.rmin() << std::endl;
00518     *out << "irncp  = " << fs.ma28_.irncp() << std::endl;
00519     *out << "icncp  = " << fs.ma28_.icncp() << std::endl;
00520     *out << "minirn = " << fs.ma28_.minirn() << std::endl;
00521     *out << "minicn = " << fs.ma28_.minicn() << std::endl;
00522     *out << "irank  = " << fs.ma28_.irank() << std::endl;
00523     *out << "themax = " << fs.ma28_.themax() << std::endl;
00524     if( fs.ma28_.lbig() == FortranTypes::F_TRUE )
00525       *out << "big    = " << fs.ma28_.big() << std::endl;
00526     *out << "ndrop  = " << fs.ma28_.ndrop() << std::endl;
00527     if( iflag >= 0 ) {
00528       *out << "\nAnalysis:\n"
00529          << "estimated_fillin_ratio can be reduced to max(minirn,minicn)/nz = "
00530          << "max(" << fs.ma28_.minirn() << "," << fs.ma28_.minicn() << ")/" << fs.nz_
00531          << " = " << my_max( fs.ma28_.minirn(), fs.ma28_.minicn() ) / (double)fs.nz_
00532          << std::endl;
00533     }
00534   }
00535 }
00536 
00537 void DirectSparseSolverMA28::ThrowIFlagException(index_type iflag)
00538 {
00539   E_IFlag e_iflag = static_cast<E_IFlag>(iflag);
00540   const char msg_err_head[] = "DirectSparseSolverMA28::ThrowIFlagException(iflag) : Error";
00541   switch(e_iflag) {
00542     case SLOW_ITER_CONV :
00543       TEST_FOR_EXCEPTION(
00544         true, std::runtime_error
00545         ,msg_err_head << ", Convergence to slow" );
00546     case MAXIT_REACHED :
00547       TEST_FOR_EXCEPTION(
00548         true, std::runtime_error
00549         ,msg_err_head << ", Maximum iterations exceeded");
00550     case MA28BD_CALLED_WITH_DROPPED :
00551       TEST_FOR_EXCEPTION(
00552         true, std::logic_error
00553         ,msg_err_head << ", ma28bd called with elements dropped in ma28ad");
00554     case DUPLICATE_ELEMENTS :
00555       TEST_FOR_EXCEPTION(
00556         true, FactorizationFailure
00557         ,msg_err_head << ", Duplicate elements have been detected");
00558     case NEW_NONZERO_ELEMENT :
00559       TEST_FOR_EXCEPTION(
00560         true, FactorizationFailure
00561         ,msg_err_head << ", A new non-zero element has be passed to ma28bd that was not ot ma28ad");
00562     case N_OUT_OF_RANGE :
00563       TEST_FOR_EXCEPTION(
00564         true, FactorizationFailure
00565         ,msg_err_head << ", 1 <=max(n,m) <= 32767 has been violated");
00566     case NZ_LE_ZERO :
00567       TEST_FOR_EXCEPTION(
00568         true, std::logic_error
00569         ,msg_err_head << ", nz <= 0 has been violated");
00570     case LICN_LE_NZ :
00571       TEST_FOR_EXCEPTION(
00572         true, std::logic_error
00573         ,msg_err_head << ", licn <= nz has been violated");
00574     case LIRN_LE_NZ :
00575       TEST_FOR_EXCEPTION(
00576         true, std::logic_error
00577         ,msg_err_head << ", lirn <= nz has been violated");
00578     case ERROR_DURRING_BLOCK_TRI :
00579       TEST_FOR_EXCEPTION(
00580         true, FactorizationFailure
00581         ,msg_err_head << ", An error has occured durring block triangularization");
00582     case LICN_AND_LIRN_TOO_SMALL :
00583       TEST_FOR_EXCEPTION(
00584         true, FactorizationFailure
00585         ,msg_err_head << ", licn and lirn are to small to hold matrix factorization");
00586     case LICN_TOO_SMALL :
00587       TEST_FOR_EXCEPTION(
00588         true, FactorizationFailure
00589         ,msg_err_head << ", licn is to small to hold matrix factorization");
00590     case LICN_FAR_TOO_SMALL :
00591       TEST_FOR_EXCEPTION(
00592         true, FactorizationFailure
00593         ,msg_err_head << ", licn is to far small to hold matrix factorization");
00594     case LIRN_TOO_SMALL :
00595       TEST_FOR_EXCEPTION(
00596         true, FactorizationFailure
00597         ,msg_err_head << ", lirn is to small to hold matrix factorization");
00598     case NUMERICALLY_SINGULAR :
00599       TEST_FOR_EXCEPTION(
00600         true, FactorizationFailure
00601         ,msg_err_head << ", matrix is numerically singular, see \'abort2\'");
00602     case STRUCTURALLY_SINGULAR :
00603       TEST_FOR_EXCEPTION(
00604         true, FactorizationFailure
00605         ,msg_err_head << ", matrix is structurally singular, see \'abort1\'");
00606     default:
00607       return; // We don't throw exceptions for other values of iflag.
00608   }
00609 }
00610 
00611 } // end namespace AbstractLinAlgPack 

Generated on Thu Sep 18 12:35:11 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1