MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_MatrixHessianSuperBasic.cpp
Go to the documentation of this file.
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include "ConstrainedOptPack_MatrixHessianSuperBasic.hpp"
00045 #include "ConstrainedOptPack_initialize_Q_R_Q_X.hpp"
00046 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00047 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00048 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00049 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOpOut.hpp"
00050 #include "DenseLinAlgPack_DVectorClass.hpp"
00051 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00052 #include "DenseLinAlgPack_AssertOp.hpp"
00053 
00054 namespace LinAlgOpPack {
00055   using AbstractLinAlgPack::Vp_StMtV;
00056 }
00057 
00058 namespace ConstrainedOptPack {
00059 
00060 MatrixHessianSuperBasic::MatrixHessianSuperBasic()
00061   : n_(0)
00062 {}
00063 
00064 void MatrixHessianSuperBasic::initialize(
00065   size_type            n
00066   ,size_type           n_R
00067   ,const size_type     i_x_free[]
00068   ,const size_type     i_x_fixed[]
00069   ,const EBounds       bnd_fixed[]
00070   ,const B_RR_ptr_t&   B_RR_ptr
00071   ,const B_RX_ptr_t&   B_RX_ptr
00072   ,BLAS_Cpp::Transp    B_RX_trans
00073   ,const B_XX_ptr_t&   B_XX_ptr
00074   )
00075 {
00076   using DenseLinAlgPack::Mp_M_assert_sizes;
00077   using BLAS_Cpp::no_trans;
00078 
00079   const size_type
00080     n_X = n - n_R;
00081 
00082     // Validate input arguments
00083 
00084   // i_x_free
00085   if( 0 < n_R && n_R < n && i_x_free == NULL ) {
00086     throw std::invalid_argument(
00087       "MatrixHessianSuperBasic::initialize(...) : Error, "
00088       "i_x_free can not be NULL when 0 < n_R < n" );
00089   }
00090   // i_x_fixed
00091   if( 0 < n_X && n_X < n && i_x_fixed == NULL ) {
00092     throw std::invalid_argument(
00093       "MatrixHessianSuperBasic::initialize(...) : Error, "
00094       "i_x_fixed can not be NULL when 0 < n-n_R < n" );
00095   }
00096   // bnd_fixed
00097   if( 0 < n_X && bnd_fixed == NULL ) {
00098     throw std::invalid_argument(
00099       "MatrixHessianSuperBasic::initialize(...) : Error, "
00100       "bnd_fixed can not be NULL when 0 < n-n_R" );
00101   }
00102   // B_RR
00103   if(n_R > 0 ) {
00104     if( !B_RR_ptr.get() )
00105       throw std::invalid_argument(
00106         "MatrixHessianSuperBasic::initialize(...) : Error, "
00107         "B_RR_ptr.get() can not be NULL when n_R > 0" );
00108     Mp_M_assert_sizes( n_R, n_R, no_trans, B_RR_ptr->rows(), B_RR_ptr->cols(), no_trans );
00109   }
00110   // op(B_RX)
00111   if( n_R < n ) {
00112     if( B_RX_ptr.get() ) {
00113       Mp_M_assert_sizes( n_R, n_X, no_trans, B_RX_ptr->rows(), B_RX_ptr->cols(), B_RX_trans );
00114     }
00115   }
00116   // B_XX
00117   if( n_R < n ) {
00118     if( !B_XX_ptr.get() )
00119       throw std::invalid_argument(
00120         "MatrixHessianSuperBasic::initialize(...) : Error, "
00121         "B_XX_ptr.get() can not be NULL if n_R < n" );
00122     Mp_M_assert_sizes( n_X, n_X, no_trans, B_XX_ptr->rows(), B_XX_ptr->cols(), no_trans );
00123   }
00124 
00125   // Setup Q_R and Q_X and validate i_x_free[] and i_x_fixed[]
00126   const bool Q_R_is_idenity = (n_R == n && i_x_fixed == NULL );
00127   if( Q_R_is_idenity ) {
00128     Q_R_row_i_.resize(0);
00129     Q_R_col_j_.resize(0);
00130   }
00131   else {
00132     Q_R_row_i_.resize(n_R);
00133     Q_R_col_j_.resize(n_R);
00134   }
00135   Q_X_row_i_.resize(n_X);
00136   Q_X_col_j_.resize(n_X);
00137   bool test_setup = true;  // ToDo: Make this an input parameter!
00138   initialize_Q_R_Q_X(
00139     n_R,n_X,i_x_free,i_x_fixed,test_setup
00140     ,!Q_R_is_idenity ? &Q_R_row_i_[0] : NULL
00141     ,!Q_R_is_idenity ? &Q_R_col_j_[0] : NULL
00142     ,&Q_R_
00143     ,n_X ? &Q_X_row_i_[0] : NULL
00144     ,n_X ? &Q_X_col_j_[0] : NULL
00145     ,&Q_X_
00146   );
00147 
00148   // Setup bnd_fixed
00149   bnd_fixed_.resize(n_X);
00150   {for(size_type i = 0; i < n_X; ++i) bnd_fixed_[i] = bnd_fixed[i]; }
00151 
00152   // Set the rest of the arguments
00153   n_           = n;
00154   n_R_         = n_R;
00155   B_RR_ptr_    = B_RR_ptr;
00156   B_RX_ptr_    = B_RX_ptr;
00157   B_RX_trans_  = B_RX_trans;
00158   B_XX_ptr_    = B_XX_ptr;
00159 
00160 }
00161 
00162 // Overridden from Matrix
00163 
00164 size_type MatrixHessianSuperBasic::rows() const
00165 {
00166   return n_;
00167 }
00168 
00169 // Overridden from MatrixOp
00170 
00171 void MatrixHessianSuperBasic::Vp_StMtV(
00172   DVectorSlice* y, value_type a, BLAS_Cpp::Transp B_trans
00173   , const DVectorSlice& x, value_type b
00174   ) const
00175 {
00176   using BLAS_Cpp::no_trans;
00177   using BLAS_Cpp::trans;
00178   using BLAS_Cpp::trans_not;
00179   using AbstractLinAlgPack::V_MtV;
00180   using LinAlgOpPack::V_MtV;
00181   assert_initialized();
00182   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, B_trans, x.size() );
00183   if( n_ == n_R_ ) {
00184     //
00185     // B = Q_R*B_RR*Q_R'
00186     //
00187     // y = b*y + a*Q_R*B_RR*Q_R'*x
00188     //
00189     if( Q_R().is_identity() ) {
00190       AbstractLinAlgPack::Vp_StMtV(y,a,*this->B_RR_ptr(),no_trans,x,b);
00191     }
00192     else {
00193       DVector Q_R_x;
00194       V_MtV( &Q_R_x, Q_R(), trans, x );
00195       AbstractLinAlgPack::Vp_StPtMtV(y,a,Q_R(),no_trans,*this->B_RR_ptr(),no_trans,Q_R_x(),b);
00196     }
00197   }
00198   else if( n_R_ == 0 ) {
00199     //
00200     // B = Q_X *B_XX * Q_X'
00201     //
00202     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00203   }
00204   else {
00205     //
00206     // B = [ Q_R  Q_X  ] * [   B_RR      op(B_RX) ] * [ Q_R' ]
00207     //                     [ op(B_RX')      B_XX  ]   [ Q_X' ]
00208     //
00209     // y = b*y + a*op(B)*x
00210     //
00211     // y = b*y + a * [ Q_R  Q_X ] * [   B_RR      op(B_RX) ] * [ Q_R' ] * x
00212     //                              [ op(B_RX')      B_XX  ]   [ Q_X' ]
00213     //
00214     // y = b*y + a*Q_R*B_RR*x_R      + a*Q_R*op(B_RX)*x_X
00215     //         + a*Q_X*op(B_RX')*x_R + a*Q_X*B_XX*x_X
00216     // where:
00217     //     x_R = Q_R'*x
00218     //     x_X = Q_X'*x
00219     //
00220     SpVector
00221       x_R,
00222       x_X;
00223     // x_R = Q_R'*x
00224     V_MtV( &x_R, Q_R(), trans, x );
00225     // x_X = Q_X'*x
00226     V_MtV( &x_X, Q_X(), trans, x );
00227     // y = b*y + a*Q_R*B_RR*x_R
00228     AbstractLinAlgPack::Vp_StPtMtV(
00229       y, a, Q_R(), no_trans, *B_RR_ptr(), no_trans, x_R(), b );
00230     // y += a*Q_R*op(B_RX)*x_X + a*Q_X*op(B_RX')*x_R
00231     if( B_RX_ptr().get() ) {
00232       AbstractLinAlgPack::Vp_StPtMtV(
00233         y, a, Q_R(), no_trans, *B_RX_ptr(), B_RX_trans(), x_X() );
00234       AbstractLinAlgPack::Vp_StPtMtV(
00235         y, a, Q_X(), no_trans, *B_RX_ptr(), trans_not(B_RX_trans()), x_R() );
00236     }
00237     // y += a*Q_X*B_XX*x_X
00238     AbstractLinAlgPack::Vp_StPtMtV(
00239       y, a, Q_X(), no_trans, *B_XX_ptr(), no_trans, x_X() );
00240   }
00241 }
00242 
00243 void MatrixHessianSuperBasic::Vp_StMtV(
00244   DVectorSlice* y, value_type a, BLAS_Cpp::Transp B_trans
00245   , const SpVectorSlice& x, value_type b
00246   ) const
00247 {
00248   using BLAS_Cpp::no_trans;
00249   using BLAS_Cpp::trans;
00250   using BLAS_Cpp::trans_not;
00251   using AbstractLinAlgPack::V_MtV;
00252   using LinAlgOpPack::V_MtV;
00253   assert_initialized();
00254   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, B_trans, x.size() );
00255   if( n_ == n_R_ ) {
00256     //
00257     // B = Q_R*B_RR*Q_R'
00258     //
00259     // y = b*y + a*Q_R*B_RR*Q_R'*x
00260     //
00261     if( Q_R().is_identity() ) {
00262       AbstractLinAlgPack::Vp_StMtV(y,a,*this->B_RR_ptr(),no_trans,x,b);
00263     }
00264     else {
00265       SpVector Q_R_x;
00266       AbstractLinAlgPack::V_MtV( &Q_R_x, Q_R(), trans, x );
00267       AbstractLinAlgPack::Vp_StPtMtV(y,a,Q_R(),no_trans,*this->B_RR_ptr(),no_trans,Q_R_x(),b);
00268     }
00269   }
00270   else if( n_R_ == 0 ) {
00271     //
00272     // B = Q_X *B_XX * Q_X'
00273     //
00274     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00275   }
00276   else {
00277     //
00278     // B = [ Q_R  Q_X  ] * [   B_RR      op(B_RX) ] * [ Q_R' ]
00279     //                     [ op(B_RX')      B_XX  ]   [ Q_X' ]
00280     //
00281     // y = b*y + a*op(B)*x
00282     //
00283     // y = b*y + a * [ Q_R  Q_X ] * [   B_RR      op(B_RX) ] * [ Q_R' ] * x
00284     //                              [ op(B_RX')      B_XX  ]   [ Q_X' ]
00285     //
00286     // y = b*y + a*Q_R*B_RR*x_R      + a*Q_R*op(B_RX)*x_X
00287     //         + a*Q_X*op(B_RX')*x_R + a*Q_X*B_XX*x_X
00288     // where:
00289     //     x_R = Q_R'*x
00290     //     x_X = Q_X'*x
00291     //
00292     SpVector
00293       x_R,
00294       x_X;
00295     // x_R = Q_R'*x
00296     V_MtV( &x_R, Q_R(), trans, x );
00297     // x_X = Q_X'*x
00298     V_MtV( &x_X, Q_X(), trans, x );
00299     // y = b*y + a*Q_R*B_RR*x_R
00300     AbstractLinAlgPack::Vp_StPtMtV(
00301       y, a, Q_R(), no_trans, *B_RR_ptr(), no_trans, x_R(), b );
00302     // y += a*Q_R*op(B_RX)*x_X + a*Q_X*op(B_RX')*x_R
00303     if( B_RX_ptr().get() ) {
00304       AbstractLinAlgPack::Vp_StPtMtV(
00305         y, a, Q_R(), no_trans, *B_RX_ptr(), B_RX_trans(), x_X() );
00306       AbstractLinAlgPack::Vp_StPtMtV(
00307         y, a, Q_X(), no_trans, *B_RX_ptr(), trans_not(B_RX_trans()), x_R() );
00308     }
00309     // y += a*Q_X*B_XX*x_X
00310     AbstractLinAlgPack::Vp_StPtMtV(
00311       y, a, Q_X(), no_trans, *B_XX_ptr(), no_trans, x_X() );
00312   }
00313 }
00314 
00315 void MatrixHessianSuperBasic::Vp_StPtMtV(
00316   DVectorSlice* y, value_type a
00317   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00318   , BLAS_Cpp::Transp M_trans
00319   , const DVectorSlice& x, value_type b ) const
00320 {
00321   using BLAS_Cpp::no_trans;
00322   using BLAS_Cpp::trans;
00323   using BLAS_Cpp::trans_not;
00324   using AbstractLinAlgPack::V_MtV;
00325   using DenseLinAlgPack::Vt_S;
00326   using LinAlgOpPack::V_MtV;
00327   namespace slap = AbstractLinAlgPack;
00328 
00329   assert_initialized();
00330 
00331   //
00332   // y = b*y + a * op(P) * B * x
00333   //
00334   // =>
00335   //
00336   // y = b*y + a * op(P)*(Q_R*B_RR*Q_R' + Q_R*op(B_RX)*Q_X' + Q_X*op(B_RX')*Q_R + Q_X*B_XX*Q_X')*x
00337   //   
00338   //   = b*y + a*op(P)*Q_R*B_RR*Q_R'*x     + a*op(P)*Q_R*op(B_RX)*Q_X'*x
00339   //         + a*op(P)*Q_X*op(B_RX')*Q_R*x + a*op(P)*Q_X*B_XX*Q_X'*x
00340   //
00341   // In order to implement the above as efficiently as possible we need to minimize the
00342   // computations with the constituent matrices.  First off we will compute
00343   // Q_RT_x = Q_R'*x (O(n_R)) and Q_XT_x = Q_X'*x (O(n_R)) neglect any terms where
00344   // Q_RT_x.nz() == 0 or Q_XT_x.nz() == 0.  We will also determine if op(P)*Q_R == 0 (O(n_R))
00345   // or op(P)*Q_X == 0 (O(n_X)) and neglect these terms if the are zero.
00346   // Hopefully this work will allow us to skip as many computations as possible.
00347   //
00348   LinAlgOpPack::Vp_MtV_assert_sizes(y->size(),P.rows(),P.cols(),P_trans
00349     , BLAS_Cpp::rows( rows(), cols(), M_trans) );
00350   LinAlgOpPack::Vp_MtV_assert_sizes( BLAS_Cpp::cols( P.rows(), P.cols(), P_trans)
00351     ,rows(),cols(),M_trans,x.size());
00352   // Q_R'*x
00353   SpVector Q_RT_x;
00354   if(n_R_) {
00355     slap::V_MtV( &Q_RT_x, Q_R(), trans, x );
00356   }
00357   // Q_X'*x
00358   SpVector Q_XT_x;
00359   if(n_ > n_R_) {
00360     slap::V_MtV( &Q_XT_x, Q_X(), trans, x );
00361   }
00362   // op(P)*Q_R overlap
00363   size_type P_Q_R_nz = 0;
00364   AbstractLinAlgPack::intersection( P, P_trans, Q_R(), no_trans, &P_Q_R_nz );
00365   // op(P)*Q_X overlap
00366   size_type P_Q_X_nz = 0;
00367   AbstractLinAlgPack::intersection( P, P_trans, Q_X(), no_trans, &P_Q_X_nz );
00368   // y = b*y
00369   if(b==0.0)      *y = 0.0;
00370   else if(b!=1.0) Vt_S(y,b);
00371   // 
00372   DVector t; // ToDo: use workspace
00373   // y += a*op(P)*Q_R*B_RR*Q_R'*x
00374   if( P_Q_R_nz && Q_RT_x.nz() ) {
00375     t.resize(n_);
00376     slap::Vp_StPtMtV( &t(), 1.0, Q_R(), no_trans, *B_RR_ptr(), no_trans, Q_RT_x() );
00377     slap::Vp_StMtV( y, a, P, P_trans, t() );
00378   }
00379   // y += a*op(P)*Q_R*op(B_RX)*Q_X'*x
00380   if( P_Q_R_nz && B_RX_ptr().get() && Q_XT_x.nz() ) {
00381     t.resize(n_);
00382     slap::Vp_StPtMtV( &t(), 1.0, Q_R(), no_trans, *B_RX_ptr(), B_RX_trans(), Q_XT_x() );
00383     slap::Vp_StMtV( y, a, P, P_trans, t() );
00384   }
00385   // y += a*op(P)*Q_X*op(B_RX')*Q_R*x
00386   if( P_Q_X_nz && B_RX_ptr().get() && Q_RT_x.nz() ) {
00387     t.resize(n_);
00388     slap::Vp_StPtMtV( &t(), 1.0, Q_X(), no_trans, *B_RX_ptr(), trans_not(B_RX_trans()), Q_RT_x() );
00389     slap::Vp_StMtV( y, a, P, P_trans, t() );
00390   }
00391   // y += a*op(P)*Q_X*B_XX*Q_X'*x
00392   if( P_Q_X_nz && Q_XT_x.nz() ) {
00393     t.resize(n_);
00394     slap::Vp_StPtMtV( &t(), 1.0, Q_X(), no_trans, *B_XX_ptr(), no_trans, Q_XT_x() );
00395     slap::Vp_StMtV( y, a, P, P_trans, t() );
00396   }
00397 }
00398 
00399 value_type MatrixHessianSuperBasic::transVtMtV(
00400   const SpVectorSlice& x1, BLAS_Cpp::Transp B_trans
00401   , const SpVectorSlice& x2 ) const
00402 {
00403   using BLAS_Cpp::no_trans;
00404   using BLAS_Cpp::trans;
00405   assert_initialized();
00406   DenseLinAlgPack::Vp_MtV_assert_sizes( x1.size(), rows(), cols(), B_trans, x1.size() );
00407   if( n_ == n_R_ ) {
00408     //
00409     // B = Q_R*B_RR*Q_R'
00410     //
00411     // a = x1'*Q_R*B_RR*Q_R'*x2
00412     //
00413     if( Q_R().is_identity() ) {
00414       return AbstractLinAlgPack::transVtMtV( x1, *B_RR_ptr(), no_trans, x2 );
00415     }
00416     else {
00417       if( x1.overlap(x2) == DenseLinAlgPack::SAME_MEM ) {
00418         SpVector Q_RT_x2;
00419         AbstractLinAlgPack::V_MtV( &Q_RT_x2, Q_R(), trans, x2 );
00420         SpVectorSlice Q_RT_x2_slc = Q_RT_x2();
00421         return AbstractLinAlgPack::transVtMtV(
00422           Q_RT_x2_slc, *B_RR_ptr(), no_trans, Q_RT_x2_slc );
00423        }
00424       else {
00425         SpVector Q_RT_x2;
00426         AbstractLinAlgPack::V_MtV( &Q_RT_x2, Q_R(), trans, x2 );
00427         SpVector Q_RT_x1;
00428         AbstractLinAlgPack::V_MtV( &Q_RT_x1, Q_R(), trans, x1 );
00429         return AbstractLinAlgPack::transVtMtV(
00430           Q_RT_x1(), *B_RR_ptr(), no_trans, Q_RT_x2() );
00431       }
00432     }
00433   }
00434   else if( n_R_ == 0 ) {
00435     //
00436     // B = Q_X *B_XX * Q_X'
00437     //
00438     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00439   }
00440   else {
00441     //
00442     // B = [ Q_R  Q_X  ] * [   B_RR      op(B_RX) ] * [ Q_R' ]
00443     //                     [ op(B_RX')      B_XX  ]   [ Q_X' ]
00444     //
00445     //
00446     // a = x1'*B*x2
00447     // =>
00448     // a = x1' * [ Q_R  Q_X  ] * [   B_RR      op(B_RX) ] * [ Q_R' ] * x2
00449     //                           [ op(B_RX')      B_XX  ]   [ Q_X' ]
00450     //
00451     // a = x1'*Q_R*B_RR*Q_R'*x2 + 2*x1'*Q_R*op(B_RX)*Q_X'*x2 + x1'*Q_X*B_XX*Q_X'*x2
00452     //
00453     if( x1.overlap(x2) == DenseLinAlgPack::SAME_MEM ) {
00454       // a = x1'*Q_R*B_RR*Q_R'*x1 + 2*x1'*Q_R*op(B_RX)*Q_X'*x1 + x1'*Q_X*B_XX*Q_X'*x1
00455       SpVector Q_RT_x1;
00456       if( Q_R().nz() )
00457         AbstractLinAlgPack::V_MtV( &Q_RT_x1, Q_R(), trans, x1 );
00458       SpVector Q_XT_x1;
00459       if( Q_X().nz() )
00460         AbstractLinAlgPack::V_MtV( &Q_XT_x1, Q_X(), trans, x1 );
00461       SpVectorSlice Q_RT_x1_slc = Q_RT_x1();
00462       SpVectorSlice Q_XT_x1_slc = Q_XT_x1();
00463       return
00464         ( Q_R().nz()
00465           ? AbstractLinAlgPack::transVtMtV(
00466             Q_RT_x1_slc, *B_RR_ptr(), no_trans, Q_RT_x1_slc )
00467           : 0.0
00468           )
00469         + 2*(  B_RX_ptr().get() && Q_R().nz() && Q_X().nz()
00470              ? AbstractLinAlgPack::transVtMtV(
00471                Q_RT_x1_slc, *B_RX_ptr(), B_RX_trans(), Q_XT_x1_slc )
00472              : 0.0
00473           )
00474         + ( Q_X().nz()
00475           ? AbstractLinAlgPack::transVtMtV(
00476             Q_XT_x1_slc, *B_XX_ptr(), no_trans, Q_XT_x1_slc )
00477           : 0.0
00478           );
00479     }
00480     else {
00481       TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
00482     }
00483   }
00484   return 0.0; // Will never be executed!
00485 }
00486 
00487 // Private
00488 
00489 void MatrixHessianSuperBasic::assert_initialized() const
00490 {
00491   if( !n_ )
00492     throw std::logic_error(
00493       "MatrixHessianSuperBasic::assert_initialized() : Error, "
00494       "The matrix is not initialized yet" );
00495 }
00496 
00497 } // end namespace ConstrainedOptPack
00498 
00499 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines