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