AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_DirectSparseSolverMA28.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include <fstream>
45 #include <algorithm>
46 
47 #include "Moocho_ConfigDefs.hpp"
48 
49 
50 #ifdef HAVE_MOOCHO_MA28
51 
52 
53 #include "AbstractLinAlgPack_DirectSparseSolverMA28.hpp"
54 #include "AbstractLinAlgPack_MatrixScaling_Strategy.hpp"
55 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
56 #include "DenseLinAlgPack_PermVecMat.hpp"
57 #include "Teuchos_AbstractFactoryStd.hpp"
58 #include "Teuchos_Assert.hpp"
59 #include "Teuchos_Workspace.hpp"
60 #include "Teuchos_dyn_cast.hpp"
61 #include "FortranTypes_f_open_file.hpp"
62 
63 namespace {
64 //
65 template< class T >
66 inline
67 T my_max( const T& v1, const T& v2 ) { return v1 > v2 ? v1 : v2; }
68 // A cast to const is needed because the standard does not return a reference from
69 // valarray<>::operator[]() const.
70 template <class T>
71 std::valarray<T>& cva(const std::valarray<T>& va )
72 {
73  return const_cast<std::valarray<T>&>(va);
74 }
75 }
76 
77 namespace AbstractLinAlgPack {
78 
79 //
80 // Implementation of DirectSparseSolver(Imp) interface using MA28.
81 //
82 // Here are some little bits of knowledge about MA28 that I need
83 // to record after many hours of hard work.
84 //
85 // * The 1979 paper in ACM TOMS (Vol. 5, No. 1, pages 27), seems
86 // to suggest that MA28 pivots by column for numerical stability
87 // but I am not sure about this.
88 //
89 // * When factoring a rectangular matrix, you must set
90 // LBLOCK = .FALSE. or the row and column permutations
91 // extracted from IKEEP(:,2) and IKEEP(:,3) respectivly
92 // are meaningless.
93 //
94 // ToDo: Finish this discussion!
95 //
96 
97 // ToDo:
98 // a) Add an option for printing the values of the common
99 // block parameters to out or to a file. This can
100 // be used to get a feel for the performance of
101 // ma28
102 // b) Add provisions for an external client to change
103 // the control options of MA28. Most of these are
104 // stored as common block variables.
105 
106 // //////////////////////////////////////////////////
107 // DirectSparseSolverMA28::FactorizationStructureMA28
108 
109 DirectSparseSolverMA28::FactorizationStructureMA28::FactorizationStructureMA28()
110  :m_(0),n_(0),max_n_(0),nz_(0),licn_(0),lirn_(0)
111  ,u_(0.1),scaling_(NO_SCALING)
112 {}
113 
114 // //////////////////////////////////////////////////
115 // DirectSparseSolverMA28::BasisMatrixMA28
116 
117 // Overridden from BasisMatrixImp
118 
121 {
122  return Teuchos::rcp(new BasisMatrixMA28);
123 }
124 
126  VectorMutable* y, BLAS_Cpp::Transp M_trans, const Vector& x
127  ) const
128 {
129  using Teuchos::dyn_cast;
130  using Teuchos::Workspace;
132  size_type k;
133 
134  // Get concrete objects
135  const FactorizationStructureMA28
136  &fs = dyn_cast<const FactorizationStructureMA28>(*this->get_fact_struc());
137  const FactorizationNonzerosMA28
138  &fn = dyn_cast<const FactorizationNonzerosMA28>(*this->get_fact_nonzeros());
139 
140  // Validate input
141 #ifdef TEUCHOS_DEBUG
143  y == NULL, std::invalid_argument
144  ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
145 #endif
146  const size_type y_dim = y->dim(), x_dim = x.dim();
147 #ifdef TEUCHOS_DEBUG
149  fs.rank_ != y_dim, std::invalid_argument
150  ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
152  fs.rank_ != x_dim, std::invalid_argument
153  ,"DirectSparseSolverMA28::BasisMatrixMA28::V_InvMtV(...) : Error! " );
154 #endif
155 
156  VectorDenseMutableEncap yd(*y);
157  VectorDenseEncap xd(x);
158 
159  // Allocate workspace memory
160  Workspace<value_type> xfull_s(wss,fs.max_n_,false);
161  DVectorSlice xfull(&xfull_s[0],xfull_s.size());
162  Workspace<value_type> ws(wss,fs.max_n_,false);
163  DVectorSlice w(&ws[0],ws.size());
164 
165  // Get a context for transpose or no transpose
166  const IVector
167  &row_perm = (M_trans == BLAS_Cpp::no_trans) ? fs.row_perm_ : fs.col_perm_,
168  &col_perm = (M_trans == BLAS_Cpp::no_trans) ? fs.col_perm_ : fs.row_perm_;
169 
170  // Copy x into positions in full w
171  // Here, the rhs vector is set with only those equations that
172  // are part of the nonsingular set. It is important that the
173  // ordering be the same as the original ordering sent to
174  // MA28AD().
175  xfull = 0.0;
176  for( k = 1; k <= x_dim; ++k )
177  xfull(row_perm(k)) = xd()(k);
178 
179  // Scale the rhs
180  if( fs.matrix_scaling_.get() )
181  fs.matrix_scaling_->scale_rhs( M_trans, xfull.raw_ptr() );
182 
183  // Solve for the rhs
184  FortranTypes::f_int mtype = ( (M_trans == BLAS_Cpp::no_trans) ? 1 : 0 );
185  fs.ma28_.ma28cd(
186  fs.max_n_, &cva(fn.a_)[0], fs.licn_, &cva(fs.icn_)[0], &cva(fs.ikeep_)[0]
187  ,xfull.raw_ptr(), w.raw_ptr(), mtype );
188 
189  // Scale the lhs
190  if( fs.matrix_scaling_.get() )
191  fs.matrix_scaling_->scale_rhs( M_trans, xfull.raw_ptr() );
192 
193  // Copy the solution into y
194  // Here, the solution vector is set with only those variables that
195  // are in the basis. It is important that the
196  // ordering be the same as the original ordering sent to
197  // MA28AD().
198  for( k = 1; k <= y_dim; ++k )
199  yd()(k) = xfull(col_perm(k));
200 
201 }
202 
203 // //////////////////////////////////////////////////
204 // DirectSparseSolverMA28
205 
206 // Constructors/initializers
207 
209  value_type estimated_fillin_ratio
210  ,value_type u
211  ,bool grow
212  ,value_type tol
213  ,index_type nsrch
214  ,bool lbig
215  ,bool print_ma28_outputs
216  ,const std::string& output_file_name
217  )
218  :estimated_fillin_ratio_(estimated_fillin_ratio)
219  ,u_(u)
220  ,grow_(grow)
221  ,tol_(tol)
222  ,nsrch_(nsrch)
223  ,lbig_(lbig)
224  ,print_ma28_outputs_(print_ma28_outputs)
225  ,output_file_name_(output_file_name)
226  ,file_output_num_(0)
227 {}
228 
229 // Overridden from DirectSparseSolver
230 
233 {
234  namespace mmp = MemMngPack;
236 }
237 
239  value_type estimated_fillin_ratio
240  )
241 {
242  estimated_fillin_ratio_ = estimated_fillin_ratio;
243 }
244 
245 // Overridden from DirectSparseSolverImp
246 
249 {
250  return Teuchos::rcp(new FactorizationStructureMA28);
251 }
252 
255 {
256  return Teuchos::rcp(new FactorizationNonzerosMA28);
257 }
258 
261  ,FactorizationStructure *fact_struc
262  ,FactorizationNonzeros *fact_nonzeros
263  ,DenseLinAlgPack::IVector *row_perm
264  ,DenseLinAlgPack::IVector *col_perm
265  ,size_type *rank
266  ,std::ostream *out
267  )
268 {
269  using Teuchos::dyn_cast;
270  typedef MatrixConvertToSparse MCTS;
271  using Teuchos::Workspace;
273 
274  if(out)
275  *out << "\nUsing MA28 to analyze and factor a new matrix ...\n";
276 
277  // Get the concrete factorization and nonzeros objects
278  FactorizationStructureMA28
279  &fs = dyn_cast<FactorizationStructureMA28>(*fact_struc);
280  FactorizationNonzerosMA28
281  &fn = dyn_cast<FactorizationNonzerosMA28>(*fact_nonzeros);
282 
283  // Set MA28 parameters
284  set_ma28_parameters(&fs);
285 
286  // Get the dimensions of things.
287  const index_type
288  m = A.rows(),
289  n = A.cols(),
290  nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
291 
292  // Validate input
294  n <= 0 || m <= 0 || m > n, std::invalid_argument
295  ,"DirectSparseSolverMA28::imp_analyze_and_factor(...) : Error!" );
296 
297  // Memorize the dimenstions for checks later
298  fs.m_ = m; fs.n_ = n; fs.nz_ = nz;
299  fs.max_n_ = my_max(fs.m_,fs.n_);
300 
301  // By default set licn and ircn equal to estimated_fillin_ratio * nz.
302  if( estimated_fillin_ratio_ < 1.0 ) {
303  if( out ) *out << "Warning, client set estimated_fillin_ratio = " << estimated_fillin_ratio_
304  << " < 1.0.\nSetting estimated_fillin_ratio = 10.0 ...\n";
305  estimated_fillin_ratio_ = 10.0;
306  }
307  if( fs.licn_ < fs.nz_ ) fs.licn_ = static_cast<index_type>(estimated_fillin_ratio_ * fs.nz_);
308  if( fs.lirn_ < fs.nz_ ) fs.lirn_ = static_cast<index_type>(estimated_fillin_ratio_ * fs.nz_);
309 
310  // Initialize structure storage
311  fs.ivect_.resize(fs.nz_); // perminatly stores nz row indexes
312  fs.jvect_.resize(fs.nz_); // perminatly stores nz column indexes
313 
314  index_type iflag = 0;
315  for( int num_fac = 0; num_fac < 5; ++num_fac ) {
316 
317  // Initialize matrix factorization storage and temporary storage
318  fs.icn_.resize(fs.licn_); // first nz entries stores the column indexes
319  fn.a_.resize(fs.licn_);
320  fs.ikeep_.resize( fs.ma28_.lblock() ? 5*fs.max_n_ : 4*fs.max_n_ + 1 );
321  Workspace<index_type> irn_tmp_(wss,fs.lirn_), iw(wss,8*fs.max_n_);
322  Workspace<value_type> w(wss,fs.max_n_);
323 
324  // Fill in the matrix information
326  MCTS::EXTRACT_FULL_MATRIX
327  ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
328  ,fs.nz_
329  ,&fn.a_[0]
330  ,fs.nz_
331  ,&fs.ivect_[0]
332  ,&fs.jvect_[0]
333  );
334  std::copy( &fs.ivect_[0], &fs.ivect_[0] + fs.nz_, &irn_tmp_[0] );
335  std::copy( &fs.jvect_[0], &fs.jvect_[0] + fs.nz_, &fs.icn_[0] );
336 
337  // Scale the matrix
338  if( fs.matrix_scaling_.get() )
339  fs.matrix_scaling_->scale_matrix(
340  fs.m_, fs.n_, fs.nz_, &fs.ivect_[0], &fs.jvect_[0], true
341  ,&fn.a_[0]
342  );
343 
344  // Analyze and factor the matrix
345  if(out)
346  *out << "\nCalling ma28ad(...) ...\n";
347  fs.ma28_.ma28ad(
348  fs.max_n_, fs.nz_, &fn.a_[0], fs.licn_, &irn_tmp_[0], fs.lirn_, &fs.icn_[0], fs.u_
349  ,&fs.ikeep_[0], &iw[0], &w[0], &iflag
350  );
351 
352  if(iflag != 0 && out)
353  *out << "\nMA28AD returned iflag = " << iflag << " != 0!\n";
354 
355  // Print MA28 outputs
356  print_ma28_outputs(true,iflag,fs,&w[0],out);
357 
358  if( iflag >= 0 ) break;
359 
360  switch( iflag ) {
361  case LICN_AND_LIRN_TOO_SMALL:
362  case LICN_TOO_SMALL:
363  case LICN_FAR_TOO_SMALL:
364  case LIRN_TOO_SMALL:
365  if(out)
366  *out << "\nWarning! iflag = " << iflag << ", LIRN and/or LICN are too small!\n"
367  << "Increasing lirn = " << fs.lirn_ << " amd licn = " << fs.licn_
368  << " by a factor of 10\n"
369  << "(try increasing estimated_fillin_ratio = " << estimated_fillin_ratio_
370  << " to a larger value next time)...\n";
371  fs.lirn_ = 10 * fs.lirn_;
372  fs.licn_ = 10 * fs.licn_;
373  break;
374  }
375  }
376 
377  // Check for errors and throw exception if you have to.
378  ThrowIFlagException(iflag);
379 
380  // Extract the basis matrix selection
381 
382  *rank = fs.ma28_.irank();
383 
384  row_perm->resize(fs.m_);
385  if( *rank < fs.m_ ) {
386  index_type
387  *row_perm_ikeep = &fs.ikeep_[fs.max_n_],
388  *row_perm_itr = &(*row_perm)(1),
389  *row_perm_last = row_perm_itr + fs.m_;
390  for(; row_perm_itr != row_perm_last;)
391  *row_perm_itr++ = abs(*row_perm_ikeep++);
392  // Sort partitions in assending order (required!)
393  std::sort(&(*row_perm)[0] , &(*row_perm)[0] + (*rank) );
394  std::sort(&(*row_perm)[0] + (*rank) , &(*row_perm)[0] + m );
395  }
396  else {
397  DenseLinAlgPack::identity_perm( row_perm );
398  }
399 
400  col_perm->resize(fs.n_);
401  if( *rank < fs.n_ ) {
402  index_type
403  *col_perm_ikeep = &fs.ikeep_[2*fs.max_n_],
404  *col_perm_itr = &(*col_perm)(1),
405  *col_perm_last = col_perm_itr + fs.n_;
406  for(; col_perm_itr != col_perm_last;)
407  *col_perm_itr++ = abs(*col_perm_ikeep++);
408  // Sort partitions in assending order (required!)
409  std::sort(&(*col_perm)[0] , &(*col_perm)[0] + (*rank) );
410  std::sort(&(*col_perm)[0] + (*rank) , &(*col_perm)[0] + n );
411  }
412  else {
413  DenseLinAlgPack::identity_perm( col_perm );
414  }
415 
416  // Set internal copy of basis selection
417  fs.row_perm_ = *row_perm;
418  fs.col_perm_ = *col_perm;
419  fs.rank_ = *rank;
420 
421 }
422 
425  ,const FactorizationStructure &fact_struc
426  ,FactorizationNonzeros *fact_nonzeros
427  ,std::ostream *out
428  )
429 {
430  using Teuchos::dyn_cast;
431  typedef MatrixConvertToSparse MCTS;
432  using Teuchos::Workspace;
434 
435  if(out)
436  *out << "\nUsing MA28 to factor a new matrix ...\n";
437 
438  // Get the concrete factorization and nonzeros objects
439  const FactorizationStructureMA28
440  &fs = dyn_cast<const FactorizationStructureMA28>(fact_struc);
441  FactorizationNonzerosMA28
442  &fn = dyn_cast<FactorizationNonzerosMA28>(*fact_nonzeros);
443 
444  // Get the dimensions of things.
445  const index_type
446  m = A.rows(),
447  n = A.cols(),
448  nz = A.num_nonzeros( MCTS::EXTRACT_FULL_MATRIX ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM );
449 
450  // Validate input
451 #ifdef TEUCHOS_DEBUG
453  m != fs.m_ || n != fs.n_ || fs.nz_ != nz, std::invalid_argument
454  ,"DirectSparseSolverMA28::imp_factor(...) : Error, "
455  "A is not compatible with matrix passed to imp_analyze_and_factor()!" );
456 #endif
457 
458  // Initialize matrix factorization storage and temporary storage
459  if(fn.a_.size() < fs.licn_) fn.a_.resize(fs.licn_);
460  Workspace<index_type> iw(wss,5*fs.max_n_);
461  Workspace<value_type> w(wss,fs.max_n_);
462 
463  // Fill in the matrix nonzeros (we already have the structure)
465  MCTS::EXTRACT_FULL_MATRIX
466  ,MCTS::ELEMENTS_ALLOW_DUPLICATES_SUM
467  ,fs.nz_
468  ,&fn.a_[0]
469  ,0
470  ,NULL
471  ,NULL
472  );
473 
474  // Scale the matrix
475  if( fs.matrix_scaling_.get() )
476  fs.matrix_scaling_->scale_matrix(
477  fs.m_, fs.n_, fs.nz_, &cva(fs.ivect_)[0], &cva(fs.jvect_)[0], false
478  ,&fn.a_[0]
479  );
480 
481  // Factor the matrix
482  index_type iflag = 0;
483  fs.ma28_.ma28bd(
484  fs.max_n_, fs.nz_, &fn.a_[0], fs.licn_, &cva(fs.ivect_)[0], &cva(fs.jvect_)[0], &cva(fs.icn_)[0]
485  ,&cva(fs.ikeep_)[0], &iw[0], &w[0], &iflag
486  );
487 
488  if(iflag != 0 && out)
489  *out << "\nMA28BD returned iflag = " << iflag << " != 0!\n";
490 
491  // Print MA28 outputs
492  print_ma28_outputs(false,iflag,fs,&w[0],out);
493 
494  // Check for errors and throw exception if you have to.
495  ThrowIFlagException(iflag);
496 
497 }
498 
499 // private
500 
501 void DirectSparseSolverMA28::set_ma28_parameters( FactorizationStructureMA28* fs )
502 {
503  // Set common block parameters
504  fs->ma28_.lblock( FortranTypes::F_FALSE ); // Do not permute to block triangular form (*** This is critical!)
505  fs->u_ = u_;
506  fs->ma28_.grow( grow_ ? FortranTypes::F_TRUE : FortranTypes::F_FALSE );
507  fs->ma28_.tol(tol_);
508  fs->ma28_.nsrch(nsrch_);
509  fs->ma28_.lbig( lbig_ ? FortranTypes::F_TRUE : FortranTypes::F_FALSE );
510  // Setup output file
511  if( output_file_name_.length() > 0 && fs->ma28_.lp() == 0 ) {
512  // Open a fortran file
513  index_type iout = 25; // Unique?
514  FortranTypes::f_open_file( iout, output_file_name_.c_str() );
515  fs->ma28_.mp(iout);
516  fs->ma28_.lp(iout);
517  }
518  else if( output_file_name_.length() == 0 && fs->ma28_.lp() ) {
519  fs->ma28_.mp(0);
520  fs->ma28_.lp(0);
521  }
522 }
523 
524 void DirectSparseSolverMA28::print_ma28_outputs(
525  bool ma28ad_bd
526  ,index_type iflag
527  ,const FactorizationStructureMA28 &fs
528  ,const value_type w[]
529  ,std::ostream *out
530  )
531 {
532  if( print_ma28_outputs_ && out ) {
533  *out << "\nReturn parameters from MA28 (call number = " << ++file_output_num_ << ")\n";
534  if( fs.ma28_.grow() == FortranTypes::F_TRUE )
535  *out << "w(1) = " << w[0] << std::endl;
536  *out << "rmin = " << fs.ma28_.rmin() << std::endl;
537  *out << "irncp = " << fs.ma28_.irncp() << std::endl;
538  *out << "icncp = " << fs.ma28_.icncp() << std::endl;
539  *out << "minirn = " << fs.ma28_.minirn() << std::endl;
540  *out << "minicn = " << fs.ma28_.minicn() << std::endl;
541  *out << "irank = " << fs.ma28_.irank() << std::endl;
542  *out << "themax = " << fs.ma28_.themax() << std::endl;
543  if( fs.ma28_.lbig() == FortranTypes::F_TRUE )
544  *out << "big = " << fs.ma28_.big() << std::endl;
545  *out << "ndrop = " << fs.ma28_.ndrop() << std::endl;
546  if( iflag >= 0 ) {
547  *out << "\nAnalysis:\n"
548  << "estimated_fillin_ratio can be reduced to max(minirn,minicn)/nz = "
549  << "max(" << fs.ma28_.minirn() << "," << fs.ma28_.minicn() << ")/" << fs.nz_
550  << " = " << my_max( fs.ma28_.minirn(), fs.ma28_.minicn() ) / (double)fs.nz_
551  << std::endl;
552  }
553  }
554 }
555 
556 void DirectSparseSolverMA28::ThrowIFlagException(index_type iflag)
557 {
558  E_IFlag e_iflag = static_cast<E_IFlag>(iflag);
559  const char msg_err_head[] = "DirectSparseSolverMA28::ThrowIFlagException(iflag) : Error";
560  switch(e_iflag) {
561  case SLOW_ITER_CONV :
563  true, std::runtime_error
564  ,msg_err_head << ", Convergence to slow" );
565  case MAXIT_REACHED :
567  true, std::runtime_error
568  ,msg_err_head << ", Maximum iterations exceeded");
569  case MA28BD_CALLED_WITH_DROPPED :
571  true, std::logic_error
572  ,msg_err_head << ", ma28bd called with elements dropped in ma28ad");
573  case DUPLICATE_ELEMENTS :
575  true, FactorizationFailure
576  ,msg_err_head << ", Duplicate elements have been detected");
577  case NEW_NONZERO_ELEMENT :
579  true, FactorizationFailure
580  ,msg_err_head << ", A new non-zero element has be passed to ma28bd that was not ot ma28ad");
581  case N_OUT_OF_RANGE :
583  true, FactorizationFailure
584  ,msg_err_head << ", 1 <=max(n,m) <= 32767 has been violated");
585  case NZ_LE_ZERO :
587  true, std::logic_error
588  ,msg_err_head << ", nz <= 0 has been violated");
589  case LICN_LE_NZ :
591  true, std::logic_error
592  ,msg_err_head << ", licn <= nz has been violated");
593  case LIRN_LE_NZ :
595  true, std::logic_error
596  ,msg_err_head << ", lirn <= nz has been violated");
597  case ERROR_DURRING_BLOCK_TRI :
599  true, FactorizationFailure
600  ,msg_err_head << ", An error has occured durring block triangularization");
601  case LICN_AND_LIRN_TOO_SMALL :
603  true, FactorizationFailure
604  ,msg_err_head << ", licn and lirn are to small to hold matrix factorization");
605  case LICN_TOO_SMALL :
607  true, FactorizationFailure
608  ,msg_err_head << ", licn is to small to hold matrix factorization");
609  case LICN_FAR_TOO_SMALL :
611  true, FactorizationFailure
612  ,msg_err_head << ", licn is to far small to hold matrix factorization");
613  case LIRN_TOO_SMALL :
615  true, FactorizationFailure
616  ,msg_err_head << ", lirn is to small to hold matrix factorization");
617  case NUMERICALLY_SINGULAR :
619  true, FactorizationFailure
620  ,msg_err_head << ", matrix is numerically singular, see \'abort2\'");
621  case STRUCTURALLY_SINGULAR :
623  true, FactorizationFailure
624  ,msg_err_head << ", matrix is structurally singular, see \'abort1\'");
625  default:
626  return; // We don't throw exceptions for other values of iflag.
627  }
628 }
629 
630 } // end namespace AbstractLinAlgPack
631 
632 
633 #endif // HAVE_MOOCHO_MA28
DirectSparseSolverMA28(value_type estimated_fillin_ratio=10.0, value_type u=0.1, bool grow=false, value_type tol=0.0, index_type nsrch=4, bool lbig=false, bool print_ma28_outputs=false, const std::string &output_file_name="")
Default constructor.
Teuchos::RCP< const Teuchos::AbstractFactory< BasisMatrix > > basis_matrix_factory_ptr_t
virtual void coor_extract_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness, const index_type len_Aval, value_type Aval[], const index_type len_Aij, index_type Arow[], index_type Acol[], const index_type row_offset=0, const index_type col_offset=0) const =0
Extract sparse elements in a coordinate data structure.
const Teuchos::RCP< FactorizationNonzeros > create_fact_nonzeros() const
const BasisMatrix::fact_struc_ptr_t & get_fact_struc() const
virtual index_type num_nonzeros(EExtractRegion extract_region, EElementUniqueness element_uniqueness) const =0
Returns the number of nonzeros in the matrix.
Teuchos::RCP< BasisMatrixImp > create_matrix() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T_To & dyn_cast(T_From &from)
const Teuchos::RCP< FactorizationStructure > create_fact_struc() const
void V_InvMtV(VectorMutable *v_lhs, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2) const
const basis_matrix_factory_ptr_t basis_matrix_factory() const
virtual void estimated_fillin_ratio(value_type estimated_fillin_ratio)=0
Set the estimate of the fill-in ratio of the factorization.
virtual size_type cols() const
Return the number of columns in the matrix.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void imp_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, const FactorizationStructure &fact_struc, FactorizationNonzeros *fact_nonzeros, std::ostream *out)
Mix-in interface for extracing explicit elements from a sparse matrix in one of several Fortran compa...
size_t size_type
void estimated_fillin_ratio(value_type estimated_fillin_ratio)
virtual size_type rows() const
Return the number of rows in the matrix.
Transp
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
void imp_analyze_and_factor(const AbstractLinAlgPack::MatrixConvertToSparse &A, FactorizationStructure *fact_struc, FactorizationNonzeros *fact_nonzeros, DenseLinAlgPack::IVector *row_perm, DenseLinAlgPack::IVector *col_perm, size_type *rank, std::ostream *out)