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

Generated on Wed May 12 21:52:28 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7