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