ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization Version of the Day
ConstrainedOptPack_initialize_Q_R_Q_X.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 <vector>
00043 
00044 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
00045 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00046 
00047 void ConstrainedOptPack::initialize_Q_R_Q_X(
00048   size_type            n_R
00049   ,size_type           n_X
00050   ,const size_type     i_x_free[]
00051   ,const size_type     i_x_fixed[]
00052   ,bool                test_setup
00053   ,size_type           Q_R_row_i[]
00054   ,size_type           Q_R_col_j[]
00055   ,GenPermMatrixSlice  *Q_R
00056   ,size_type           Q_X_row_i[]
00057   ,size_type           Q_X_col_j[]
00058   ,GenPermMatrixSlice  *Q_X
00059   )
00060 {
00061   namespace GPMSIP = AbstractLinAlgPack::GenPermMatrixSliceIteratorPack;
00062   const size_type
00063     n = n_R + n_X;
00064   //
00065   //                   /  i_x_free[lR-1]    ,if lR = l_x_XR[i] > 0
00066   // Setup l_x_XR[i] = |  not set yet       ,if l_x_XR[i] == 0
00067   //                   \  i_x_fixed[lX-1]   ,if lX = l_x_XR[i] < 0
00068   //
00069   typedef std::vector<long int> l_x_XR_t;   // ToDo: use temporary workspace
00070   l_x_XR_t l_x_XR(n);
00071   std::fill( l_x_XR.begin(), l_x_XR.end(), 0 );
00072   // Setup the fixed portion of l_x_XR[] first.
00073   bool i_x_fixed_is_sorted = true;
00074   if(test_setup) {
00075     size_type last_set = 0;
00076     const size_type *i_x_X = i_x_fixed;
00077     for( long int lX = 1; lX <= n_X; ++lX, ++i_x_X ) {
00078       if( *i_x_X < 1 || *i_x_X > n ) {
00079         std::ostringstream omsg;
00080         omsg
00081           << "initialize_Q_R_Q_X(...) : Error, "
00082           << "i_x_fixed[" << (lX-1) << "] = "
00083           << (*i_x_X) << " is out of bounds";
00084         throw std::invalid_argument( omsg.str() );
00085       }
00086       const long int l = l_x_XR[*i_x_X-1];
00087       if( l != 0 ) {
00088         std::ostringstream omsg;
00089         omsg
00090           << "initialize_Q_R_Q_X(...) : Error, "
00091           << "Duplicate entries for i_x_fixed[" << (lX-1) << "] = "
00092           << (*i_x_X) << " == i_x_fixed[" << (-l-1) << "]";
00093         throw std::invalid_argument( omsg.str() );
00094       }
00095       l_x_XR[*i_x_X-1] = -lX;
00096       if( *i_x_X < last_set )
00097         i_x_fixed_is_sorted = false;
00098       last_set = *i_x_X;
00099     }
00100   }
00101   else {
00102     size_type last_set = 0;
00103     const size_type *i_x_X = i_x_fixed;
00104     for( size_type lX = 1; lX <= n_X; ++lX, ++i_x_X ) {
00105       l_x_XR[*i_x_X-1] = -lX;
00106       if( *i_x_X < last_set )
00107         i_x_fixed_is_sorted = false;
00108       last_set = *i_x_X;
00109     }
00110   }
00111   // Now setup the free portion of l_x_XR[]
00112   bool i_x_free_is_sorted = true;
00113   if( i_x_free == NULL ) {
00114     long int
00115       *l_x_XR_itr = &l_x_XR[0];
00116     size_type lR = 0;
00117     for( size_type i = 1; i <= n; ++i, ++l_x_XR_itr ) {
00118       if(*l_x_XR_itr == 0) {
00119         ++lR;
00120         *l_x_XR_itr = lR;
00121       }
00122     }
00123     TEUCHOS_TEST_FOR_EXCEPT( !( lR == n_R ) );
00124   }
00125   else {
00126     if(test_setup) {
00127       size_type last_set = 0;
00128       const size_type *i_x_R = i_x_free;
00129       for( size_type lR = 1; lR <= n_R; ++lR, ++i_x_R ) {
00130         if( *i_x_R < 1 || *i_x_R > n ) {
00131           std::ostringstream omsg;
00132           omsg
00133             << "initialize_Q_R_Q_X(...) : Error, "
00134             << "i_x_free[" << (lR-1) << "] = "
00135             << (*i_x_R) << " is out of bounds";
00136           throw std::invalid_argument( omsg.str() );
00137         }
00138         const long int l = l_x_XR[*i_x_R-1];
00139         if( l != 0 ) {
00140           std::ostringstream omsg;
00141           omsg
00142             << "initialize_Q_R_Q_X(...) : Error, "
00143             << "Duplicate entries for i_x_free[" << (lR-1) << "] = "
00144             << (*i_x_R) << " == "
00145             << (l < 0 ? "i_x_fixed[" : "i_x_free[")
00146             << (l < 0 ? -l : l) << "]";
00147           throw std::invalid_argument( omsg.str() );
00148         }
00149         l_x_XR[*i_x_R-1] = lR;
00150         if( *i_x_R < last_set )
00151           i_x_free_is_sorted = false;
00152         last_set = *i_x_R;
00153       }
00154     }
00155     else {
00156       size_type last_set = 0;
00157       const size_type *i_x_R = i_x_free;
00158       for( size_type lR = 1; lR <= n_R; ++lR, ++i_x_R ) {
00159         l_x_XR[*i_x_R-1] = lR;
00160         if( *i_x_R < last_set )
00161           i_x_free_is_sorted = false;
00162         last_set = *i_x_R;
00163       }
00164     }
00165   }
00166   //
00167   // Setup Q_R and Q_X
00168   //
00169   if( i_x_free_is_sorted && Q_R_row_i == NULL && l_x_XR[n_R-1] == n_R ) {
00170     //
00171     // Here Q_R = [ I; 0 ] so lets set it
00172     //
00173     Q_R->initialize(n,n_R,GenPermMatrixSlice::IDENTITY_MATRIX);
00174     // Now we just have to set Q_X which may or may not be sorted
00175     const long int
00176       *l_x_X = &l_x_XR[n_R-1] + 1; // Just in case n_X == 0
00177     size_type
00178       *row_i_itr = Q_X_row_i,
00179       *col_j_itr = Q_X_col_j;
00180     for( size_type iX = n_R+1; iX <= n; ++iX, ++l_x_X, ++row_i_itr, ++col_j_itr ) {
00181       *row_i_itr = iX;         // This is sorted of course
00182       *col_j_itr = -(*l_x_X);  // This might be sorted
00183       TEUCHOS_TEST_FOR_EXCEPT( !(  *l_x_X < 0  ) );
00184     }         
00185     Q_X->initialize(
00186       n,n_X,n_X,0,0,i_x_fixed_is_sorted?GPMSIP::BY_ROW_AND_COL:GPMSIP::BY_ROW
00187       ,Q_X_row_i,Q_X_col_j,test_setup);
00188   }
00189   else {
00190     //
00191     // Here i_x_free[] and i_x_fixed[] may be sorted by Q_R != [ I; 0 ]
00192     //
00193     if( n_X > 0 && Q_R_row_i == NULL )
00194       throw std::invalid_argument(
00195         "initialize_Q_R_Q_X(...) : Error, "
00196         "Q_R_row_i != NULL since Q_R != [ I; 0 ]" );
00197     const long int
00198       *l_x = &l_x_XR[0];
00199     size_type
00200       *Q_R_row_i_itr = Q_R_row_i,   // Okay if NULL and n_R == 0
00201       *Q_R_col_j_itr = Q_R_col_j,
00202       *Q_X_row_i_itr = Q_X_row_i,   // Okay if NULL and n_X == 0
00203       *Q_X_col_j_itr = Q_X_col_j;
00204     for( size_type i = 1; i <= n; ++i, ++l_x ) {
00205       const long int l = *l_x;
00206       TEUCHOS_TEST_FOR_EXCEPT( !(  l != 0  ) );
00207       if( l > 0 ) {
00208         *Q_R_row_i_itr++ = i;
00209         *Q_R_col_j_itr++ = l;
00210       }
00211       else {
00212         *Q_X_row_i_itr++ = i;
00213         *Q_X_col_j_itr++ = -l;
00214       }
00215     }         
00216     TEUCHOS_TEST_FOR_EXCEPT( !(  Q_R_row_i_itr - Q_R_row_i == n_R  ) );
00217     TEUCHOS_TEST_FOR_EXCEPT( !(  Q_X_row_i_itr - Q_X_row_i == n_X  ) );
00218     Q_R->initialize(
00219       n,n_R,n_R,0,0,i_x_free_is_sorted?GPMSIP::BY_ROW_AND_COL:GPMSIP::BY_ROW
00220       ,Q_R_row_i,Q_R_col_j,test_setup);
00221     Q_X->initialize(
00222       n,n_X,n_X,0,0,i_x_fixed_is_sorted?GPMSIP::BY_ROW_AND_COL:GPMSIP::BY_ROW
00223       ,Q_X_row_i,Q_X_col_j,test_setup);
00224   }
00225 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends