AbstractLinAlgPack_GenPermMatrixSlice.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 <sstream>
00032 #include <functional>
00033 #include <algorithm>
00034 
00035 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00036 #include "Teuchos_TestForException.hpp"
00037 
00038 #ifdef _WINDOWS
00039 
00040 namespace std {
00041 
00042 // Help some compilers lookup the function.
00043 inline
00044 void swap(
00045     AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<
00046       DenseLinAlgPack::size_type>& v1
00047   , AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<
00048     DenseLinAlgPack::size_type>& v2
00049   )
00050 {
00051   AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::swap(v1,v2);
00052 }
00053 
00054 } // end namespace std
00055 
00056 #endif // _WINDOWS
00057 
00058 namespace {
00059 
00060 //
00061 template< class T >
00062 inline
00063 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00064 
00065 // Return a string with the same name as the enumeration
00066 const char* ordered_by_str(
00067   AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::EOrderedBy ordered_by )
00068 {
00069   switch( ordered_by ) {
00070     case AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::BY_ROW:
00071       return "BY_ROW";
00072     case AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::BY_COL:
00073       return "BY_COL";
00074     case AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::BY_ROW_AND_COL:
00075       return "BY_ROW_AND_COL";
00076     case AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::UNORDERED:
00077       return "UNORDERED";
00078   }
00079   TEST_FOR_EXCEPT(true);  // should never be executed
00080   return 0;
00081 }
00082 
00083 // Define function objects for comparing by row and by column
00084 
00085 template<class T>
00086 struct imp_row_less
00087    : public std::binary_function<
00088      AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<T>
00089      ,AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<T>
00090      ,bool
00091      >
00092 {
00093   bool operator()(
00094       const AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::external_row_col_value_type<T>& v1
00095     , const AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::external_row_col_value_type<T>& v2
00096     )
00097   {
00098     return v1.row_i_ < v2.row_i_;
00099   }
00100 };
00101 
00102 template<class T>
00103 struct imp_col_less
00104    : public std::binary_function<
00105      AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<T>
00106      ,AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::row_col_value_type<T>
00107      ,bool
00108      >
00109 {
00110   bool operator()(
00111       const AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::external_row_col_value_type<T>& v1
00112     , const AbstractLinAlgPack::GenPermMatrixSliceIteratorPack::external_row_col_value_type<T>& v2
00113     )
00114   {
00115     return v1.col_j_ < v2.col_j_;
00116   }
00117 };
00118 
00119 } // end namespace
00120 
00121 //#ifdef _WINDOWS
00122 //namespace std {
00123 //  using imp_row_less;
00124 //}
00125 //#endif // _WINDOWS
00126 
00127 namespace AbstractLinAlgPack {
00128 
00129 GenPermMatrixSlice::GenPermMatrixSlice()
00130   : rows_(0), cols_(0), nz_(0)
00131 {}
00132 
00133 void GenPermMatrixSlice::initialize( size_type rows, size_type cols, EIdentityOrZero type )
00134 {
00135   rows_       = rows;
00136   cols_       = cols;
00137   nz_         = type == IDENTITY_MATRIX ? my_min(rows,cols) : 0;
00138   ordered_by_ = GenPermMatrixSliceIteratorPack::BY_ROW_AND_COL;
00139   row_i_      = NULL;
00140   col_j_      = NULL;   // Don't need to be set but might as well for safely
00141   row_off_    = 0;      // ...
00142   col_off_    = 0;      // ...
00143 }
00144 
00145 void GenPermMatrixSlice::initialize(
00146   size_type     rows
00147   ,size_type      cols
00148   ,size_type      nz
00149   ,difference_type  row_off
00150   ,difference_type  col_off
00151   ,EOrderedBy     ordered_by
00152   ,const size_type  row_i[]
00153   ,const size_type  col_j[]
00154   ,bool       test_setup
00155   )
00156 {
00157   namespace GPMSIP = GenPermMatrixSliceIteratorPack;
00158 
00159   if( test_setup ) {
00160     std::ostringstream omsg;
00161     omsg << "\nGenPermMatrixSlice::initialize(...) : Error: ";
00162     // Validate the input data (at least partially)
00163     validate_input_data(rows,cols,nz,row_off,col_off,ordered_by,row_i,col_j,omsg);
00164     // Validate the ordering and uniquness
00165     const size_type *ordered_sequence = NULL;
00166     if( ordered_by == GPMSIP::BY_ROW || ordered_by == GPMSIP::BY_ROW_AND_COL ) {
00167       for( size_type k = 1; k < nz; ++k ) {
00168         TEST_FOR_EXCEPTION(
00169           row_i[k-1] >= row_i[k], std::invalid_argument
00170           ,"GenPermMatrixSlice::initialize(...) : Error: "
00171           "row_i[" << k-1 << "] = " << row_i[k-1]
00172           << " >= row_i[" << k << "] = " << row_i[k]
00173           << "\nThis is not sorted by row!" );
00174       }
00175     }
00176     if( ordered_by == GPMSIP::BY_COL || ordered_by == GPMSIP::BY_ROW_AND_COL ) {
00177       for( size_type k = 1; k < nz; ++k ) {
00178         TEST_FOR_EXCEPTION(
00179           col_j[k-1] >= col_j[k], std::invalid_argument
00180           ,"GenPermMatrixSlice::initialize(...) : Error: "
00181           "col_j[" << k-1 << "] = " << col_j[k-1]
00182           << " >= col_j[" << k << "] = " << col_j[k]
00183           << "\nThis is not sorted by column!" );
00184       }
00185     }
00186   }
00187   // Set the members
00188   rows_   = rows;
00189   cols_   = cols;
00190   nz_     = nz;
00191   row_off_  = row_off;
00192   col_off_  = col_off;
00193   ordered_by_ = ordered_by;
00194   row_i_    = nz ? row_i : NULL;
00195   col_j_    = nz ? col_j : NULL;
00196 }
00197 
00198 void GenPermMatrixSlice::initialize_and_sort(
00199   size_type     rows
00200   ,size_type      cols
00201   ,size_type      nz
00202   ,difference_type  row_off
00203   ,difference_type  col_off
00204   ,EOrderedBy     ordered_by
00205   ,size_type      row_i[]
00206   ,size_type      col_j[]
00207   ,bool       test_setup
00208   )
00209 {
00210   namespace GPMSIP = GenPermMatrixSliceIteratorPack;
00211   TEST_FOR_EXCEPTION(
00212     ordered_by == GPMSIP::BY_ROW_AND_COL, std::invalid_argument
00213     ,"GenPermMatrixSlice::initialize_and_sort(...) : Error, "
00214     "ordered_by == GPMSIP::BY_ROW_AND_COL, we can not sort by row and column!" );
00215   if( test_setup ) {
00216     std::ostringstream omsg;
00217     omsg << "\nGenPermMatrixSlice::initialize_and_sort(...) : Error:\n";
00218     // Validate the input data (at least partially)
00219     validate_input_data(rows,cols,nz,row_off,col_off,ordered_by,row_i,col_j,omsg);
00220   }
00221 
00222   // Sort by row or column
00223   typedef GPMSIP::row_col_iterator<size_type> row_col_itr_t;
00224   row_col_itr_t
00225     row_col_itr = row_col_itr_t( row_off, col_off, row_i, col_j, nz );
00226   if( ordered_by == GPMSIP::BY_ROW ) {
00227     std::stable_sort( row_col_itr, row_col_itr + nz
00228       , imp_row_less<size_type>() );
00229   }
00230   else if( ordered_by == GPMSIP::BY_COL ) {
00231     std::stable_sort( row_col_itr, row_col_itr + nz
00232       , imp_col_less<size_type>() );
00233   }
00234 
00235   initialize(rows,cols,nz,row_off,col_off,ordered_by,row_i,col_j,test_setup);
00236 }
00237   
00238 void GenPermMatrixSlice::bind( const GenPermMatrixSlice& gpms )
00239 {
00240   this->initialize( gpms.rows_, gpms.cols_, gpms.nz_, gpms.row_off_
00241     , gpms.col_off_, gpms.ordered_by_, gpms.row_i_, gpms.col_j_
00242     , false );
00243 }
00244 
00245 size_type GenPermMatrixSlice::lookup_row_i(size_type col_j) const
00246 {
00247   namespace QPMSIP = GenPermMatrixSliceIteratorPack;
00248   if( col_j < 1 || cols_ < col_j )
00249     std::out_of_range(
00250       "GenPermMatrixSlice::lookup_row_i(col_j) : Error, "
00251       "col_j is out of bounds" );
00252   if(!nz_)
00253     return 0;
00254   if(is_identity())
00255     return col_j <= nz_ ? col_j : 0;
00256   const size_type
00257     *itr = NULL;
00258   if( ordered_by() == QPMSIP::BY_COL || ordered_by() == QPMSIP::BY_ROW_AND_COL )
00259     itr = std::lower_bound( col_j_, col_j_ + nz_, col_j );
00260   else
00261     itr = std::find( col_j_, col_j_ + nz_, col_j );
00262   return (itr != col_j_ + nz_ && *itr == col_j) ? row_i_[itr - col_j_] : 0;
00263 }
00264 
00265 size_type GenPermMatrixSlice::lookup_col_j(size_type row_i) const
00266 {
00267   namespace QPMSIP = GenPermMatrixSliceIteratorPack;
00268   if( row_i < 1 || rows_ < row_i )
00269     std::out_of_range(
00270       "GenPermMatrixSlice::lookup_col_j(row_i) : Error, "
00271       "row_i is out of bounds" );
00272   if(!nz_)
00273     return 0;
00274   if(is_identity())
00275     return row_i <= nz_ ? row_i : 0;
00276   const size_type
00277     *itr = NULL;
00278   if( ordered_by() == QPMSIP::BY_ROW || ordered_by() == QPMSIP::BY_ROW_AND_COL )
00279     itr = std::lower_bound( row_i_, row_i_ + nz_, row_i );
00280   else
00281     itr = std::find( row_i_, row_i_ + nz_, row_i );
00282   return (itr != row_i_ + nz_ && *itr == row_i) ? col_j_[itr - row_i_] : 0;
00283 }
00284 
00285 GenPermMatrixSlice::const_iterator GenPermMatrixSlice::begin() const
00286 {
00287   validate_not_identity();
00288   return const_iterator(row_off_,col_off_,row_i_,col_j_,nz_);
00289 }
00290 
00291 GenPermMatrixSlice::const_iterator GenPermMatrixSlice::end() const
00292 {
00293   validate_not_identity();
00294   return const_iterator(row_off_,col_off_,row_i_+nz_,col_j_+nz_,0);
00295 }
00296 
00297 const GenPermMatrixSlice GenPermMatrixSlice::create_submatrix(
00298   const Range1D& rng, EOrderedBy ordered_by ) const
00299 {
00300   namespace GPMSIP = GenPermMatrixSliceIteratorPack;
00301 
00302   validate_not_identity();
00303 
00304   // Validate the input
00305   TEST_FOR_EXCEPTION(
00306     ordered_by == GPMSIP::BY_ROW_AND_COL, std::invalid_argument
00307     ,"GenPermMatrixSlice::initialize_and_sort(...) : Error, "
00308     "ordered_by == GPMSIP::BY_ROW_AND_COL, we can not sort by row and column!" );
00309   TEST_FOR_EXCEPTION(
00310     rng.full_range(), std::logic_error,
00311     "GenPermMatrixSlice::create_submatrix(...) : Error, "
00312     "The range argument can not be rng.full_range() == true" );
00313   TEST_FOR_EXCEPTION(
00314     ordered_by == GPMSIP::BY_ROW && rng.ubound() > rows(), std::logic_error
00315     ,"GenPermMatrixSlice::create_submatrix(...) : Error, "
00316     "rng.ubound() can not be larger than this->rows()" );
00317   TEST_FOR_EXCEPTION(
00318     ordered_by == GPMSIP::BY_COL && rng.ubound() > cols(), std::logic_error
00319     ,"GenPermMatrixSlice::create_submatrix(...) : Error, "
00320     "rng.ubound() can not be larger than this->cols()" );
00321   TEST_FOR_EXCEPTION(
00322     ordered_by == GPMSIP::UNORDERED, std::logic_error
00323     ,"GenPermMatrixSlice::create_submatrix(...) : Error, "
00324     "You can have ordered_by == GPMSIP::UNORDERED" );
00325 
00326   // Find the upper and lower k for row[k], col[k] indices
00327   size_type
00328     k_l = 0,    // zero based 
00329     k_u = nz() + 1; // zero based (== nz + 1 flag that there are no nonzeros to even search)
00330   size_type
00331     rows = 0,
00332     cols = 0;
00333   difference_type
00334     row_off = 0,
00335     col_off = 0;
00336   switch( ordered_by ) {
00337     case GPMSIP::BY_ROW:
00338     case GPMSIP::BY_COL:
00339     {
00340       TEST_FOR_EXCEPTION(
00341         this->ordered_by() != GPMSIP::BY_ROW_AND_COL
00342         && ( nz() > 1 && ordered_by != this->ordered_by() )
00343         ,std::logic_error
00344         ,"GenPermMatrixSlice::create_submatrix(...) : Error, "
00345         << "nz = " << nz() << " > 1 and "
00346         << "ordered_by = " << ordered_by_str(ordered_by)
00347         << " != this->ordered_by() = "
00348         << ordered_by_str(this->ordered_by()) );
00349       // Search the rows or the columns.
00350       const size_type
00351         *search_k = NULL;
00352       difference_type
00353         search_k_off;
00354       if( ordered_by == GPMSIP::BY_ROW ) {
00355         search_k = row_i_;  // may be null
00356         search_k_off = row_off_;
00357         rows = rng.size();
00358         cols = this->cols();
00359         row_off = row_off_ - (difference_type)(rng.lbound() - 1);
00360         col_off = col_off_;
00361       }
00362       else {  // BY_COL
00363         search_k = col_j_;  // may be null
00364         search_k_off = col_off_;
00365         rows = this->rows();
00366         cols = rng.size();
00367         row_off = row_off_;
00368         col_off = col_off_ - (difference_type)(rng.lbound() - 1);;
00369       }
00370       if( search_k ) {
00371         const size_type
00372           *l = std::lower_bound( search_k, search_k + nz()
00373               , rng.lbound() - search_k_off );
00374         k_l = l - search_k;
00375         // If we did not find the lower bound in the range, don't even bother
00376         // looking for the upper bound.
00377         if( k_l != nz() ) {
00378           const size_type
00379             *u = std::upper_bound( search_k, search_k + nz()
00380                 , rng.ubound() - search_k_off );
00381           k_u = u - search_k;
00382           // Here, if there are any nonzero elements in this range then
00383           // k_u - k_l > 0 will always be true!
00384         }
00385       }
00386       break;
00387     }
00388     case GPMSIP::UNORDERED:
00389       TEST_FOR_EXCEPT(true);
00390   }
00391   GenPermMatrixSlice gpms;
00392   if( k_u - k_l > 0 && k_u != nz() + 1 ) {
00393     gpms.initialize(
00394         rows, cols
00395       , k_u - k_l
00396       , row_off, col_off
00397       , ordered_by
00398       , row_i_ + k_l
00399       , col_j_ + k_l
00400       );
00401   }
00402   else  {
00403     gpms.initialize(
00404         rows, cols
00405       , 0
00406       , row_off, col_off
00407       , ordered_by
00408       , NULL
00409       , NULL
00410       );
00411   }
00412   return gpms;
00413 }
00414 
00415 // static members
00416 
00417 void GenPermMatrixSlice::validate_input_data(
00418   size_type     rows
00419   ,size_type      cols
00420   ,size_type      nz
00421   ,difference_type  row_off
00422   ,difference_type  col_off
00423   ,EOrderedBy     ordered_by
00424   ,const size_type  row_i[]
00425   ,const size_type  col_j[]
00426   ,std::ostringstream &omsg
00427   )
00428 {
00429   namespace GPMSIP = GenPermMatrixSliceIteratorPack;
00430 
00431   TEST_FOR_EXCEPTION(
00432     nz > rows * cols, std::invalid_argument
00433     ,omsg.str() << "nz = " << nz << " can not be greater than rows * cols = "
00434     << rows << " * " << cols << " = " << rows * cols );
00435   
00436   // First see if everything is in range.
00437   for( size_type k = 0; k < nz; ++k ) {
00438     TEST_FOR_EXCEPTION(
00439       row_i[k] + row_off < 1 || rows < row_i[k] + row_off, std::invalid_argument
00440       ,omsg.str() << "row_i[" << k << "] + row_off = " << row_i[k] << " + " << row_off
00441       << " = " << (row_i[k] + row_off)
00442       << " is out of range [1,rows] = [1," << rows << "]" );
00443     TEST_FOR_EXCEPTION(
00444       col_j[k] + col_off < 1 || cols < col_j[k] + col_off, std::invalid_argument
00445       ,omsg.str() << "col_j[" << k << "] + col_off = " << col_j[k] << " + " << col_off
00446       << " = " << (col_j[k] + col_off)
00447       << " is out of range [1,cols] = [1," << cols << "]" );
00448   }
00449   
00450   // ToDo: Technically, we need to validate that the nonzero elements row_i[k], col_j[k] are
00451   // unique but this is much harder to do!
00452   
00453 }
00454 
00455 // private
00456 
00457 void GenPermMatrixSlice::validate_not_identity() const
00458 {
00459   TEST_FOR_EXCEPTION(
00460     is_identity(), std::logic_error
00461     ,"GenPermMatrixSlice::validate_not_identity() : "
00462     "Error, this->is_identity() is true" );
00463 }
00464 
00465 } // end namespace AbstractLinAlgPack

Generated on Tue Oct 20 12:51:42 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7