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