MOOCHO (Single Doxygen Collection) Version of the Day
ConstrainedOptPack_MatrixSymPosDefBandedChol.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 <assert.h>
00045 
00046 #include <sstream>
00047 
00048 #include "ConstrainedOptPack_MatrixSymPosDefBandedChol.hpp"
00049 #include "DenseLinAlgPack_AssertOp.hpp"
00050 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00051 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
00052 #include "MiReleaseResource_ref_count_ptr.h"
00053 #include "MiWorkspacePack.h"
00054 
00055 // LAPACK functions
00056 
00057 extern "C" {
00058 
00059 FORTRAN_FUNC_DECL_UL(void,DPBTRF,dpbtrf)(
00060   FORTRAN_CONST_CHAR_1_ARG(UPLO)
00061   ,const FortranTypes::f_int       &N
00062   ,const FortranTypes::f_int       &KD
00063   ,FortranTypes::f_dbl_prec        *AB
00064   ,const FortranTypes::f_int       &LDAB
00065   ,FortranTypes::f_int             *INFO
00066   );
00067 
00068 FORTRAN_FUNC_DECL_UL(void,DPBTRS,dpbtrs)(
00069   FORTRAN_CONST_CHAR_1_ARG(UPLO)
00070   ,const FortranTypes::f_int       &N
00071   ,const FortranTypes::f_int       &KD
00072   ,const FortranTypes::f_int       &NRHS
00073   ,const FortranTypes::f_dbl_prec  AB[]
00074   ,const FortranTypes::f_int       &LDAB
00075   ,FortranTypes::f_dbl_prec        *B
00076   ,const FortranTypes::f_int       &LDB
00077   ,FortranTypes::f_int             *INFO
00078   );
00079 
00080 } // end namespace LAPACK
00081 
00082 namespace ConstrainedOptPack {
00083 
00084 MatrixSymPosDefBandedChol::MatrixSymPosDefBandedChol(
00085   size_type                         n
00086   ,size_type                        kd
00087   ,DMatrixSlice                   *MB
00088   ,const release_resource_ptr_t&    MB_release_resource_ptr
00089   ,BLAS_Cpp::Uplo                   MB_uplo
00090   ,DMatrixSlice                   *UB
00091   ,const release_resource_ptr_t&    UB_release_resource_ptr
00092   ,BLAS_Cpp::Uplo                   UB_uplo
00093   ,bool                             update_factor
00094   ,value_type                       scale
00095   )
00096 {
00097   initialize(n,kd,MB,MB_release_resource_ptr,MB_uplo
00098          ,UB,UB_release_resource_ptr,UB_uplo,update_factor,scale);
00099 }
00100 
00101 void MatrixSymPosDefBandedChol::initialize(
00102   size_type                         n
00103   ,size_type                        kd
00104   ,DMatrixSlice                   *MB
00105   ,const release_resource_ptr_t&    MB_release_resource_ptr
00106   ,BLAS_Cpp::Uplo                   MB_uplo
00107   ,DMatrixSlice                   *UB
00108   ,const release_resource_ptr_t&    UB_release_resource_ptr
00109   ,BLAS_Cpp::Uplo                   UB_uplo
00110   ,bool                             update_factor
00111   ,value_type                       scale
00112   )
00113 {
00114   // Validate input
00115 
00116   if( n == 0 ) {
00117     if( kd != 0 )
00118       throw std::invalid_argument(
00119         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00120         "kd must be 0 if n == 0" );
00121     if( MB != NULL )
00122       throw std::invalid_argument(
00123         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00124         "MB must be NULL if n == 0" );
00125     if( MB_release_resource_ptr.get() != NULL )
00126       throw std::invalid_argument(
00127         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00128         "MB_release_resource_ptr.get() must be NULL if n == 0" );
00129     if( UB != NULL )
00130       throw std::invalid_argument(
00131         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00132         "UB must be NULL if n == 0" );
00133     if( UB_release_resource_ptr.get() != NULL )
00134       throw std::invalid_argument(
00135         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00136         "UB_release_resource_ptr.get() must be NULL if n == 0" );
00137   }
00138   else {
00139     if( kd + 1 > n )
00140       throw std::invalid_argument(
00141         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00142         "kd + 1 can not be larger than n" );
00143     if( MB == NULL && UB == NULL )
00144       throw std::invalid_argument(
00145         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00146         "MB and UB can not both be NULL" );
00147     if( MB != NULL && ( MB->rows() != kd + 1 || MB->cols() != n ) )
00148       throw std::invalid_argument(
00149         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00150         "MB is not the correct size" );
00151     if( UB != NULL && ( UB->rows() != kd + 1 || UB->cols() != n ) )
00152       throw std::invalid_argument(
00153         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00154         "UB is not the correct size" );
00155   }
00156 
00157   // Set the members
00158 
00159   if( n == 0 ) {
00160     n_                        = 0;
00161     kd_                       = 0;
00162     MB_.bind(DMatrixSlice());
00163     MB_release_resource_ptr_  = NULL;
00164     MB_uplo_                  = BLAS_Cpp::lower;
00165     UB_.bind(DMatrixSlice());
00166     UB_release_resource_ptr_  = NULL;
00167     UB_uplo_                  = BLAS_Cpp::lower;
00168     scale_                    = 1.0;
00169   }
00170   else {
00171     // Set the members
00172     n_                        = n;
00173     kd_                       = kd;
00174     if(MB) MB_.bind(*MB);
00175     MB_release_resource_ptr_  = MB_release_resource_ptr;
00176     MB_uplo_                  = MB_uplo;
00177     if(UB) UB_.bind(*UB);
00178     UB_release_resource_ptr_  = UB_release_resource_ptr;
00179     UB_uplo_                  = UB_uplo;
00180     factor_updated_           = UB && !update_factor;
00181     scale_                    = scale;
00182     // Update the factorization if we have storage
00183     if( update_factor )
00184       update_factorization();
00185   }
00186 }
00187 
00188 // Overridden from MatrixOp
00189 
00190 size_type MatrixSymPosDefBandedChol::rows() const
00191 {
00192   return n_;
00193 }
00194 
00195 size_type MatrixSymPosDefBandedChol::nz() const
00196 {
00197   return (2 * kd_ + 1) * n_ - ( (kd_+1)*(kd_+1) - (kd_+1) );
00198 }
00199 
00200 std::ostream& MatrixSymPosDefBandedChol::output(std::ostream& out) const
00201 {
00202   return MatrixOp::output(out); // ToDo: Implement specialized version later!
00203 }
00204 
00205 void MatrixSymPosDefBandedChol::Vp_StMtV(
00206   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00207   , const DVectorSlice& x, value_type b) const
00208 {
00209   assert_initialized();
00210   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, BLAS_Cpp::no_trans, x.size() );
00211   if( MB_.rows() ) {
00212     BLAS_Cpp::sbmv(MB_uplo_,n_,kd_,a,MB_.col_ptr(1),MB_.max_rows(),x.raw_ptr(),x.stride()
00213              ,b,y->raw_ptr(),y->stride());
00214   }
00215   else if( UB_.rows() ) {
00216     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement when and if needed!
00217   }
00218   else {
00219     TEUCHOS_TEST_FOR_EXCEPT(true); // This should not happen!
00220   }
00221 }
00222 
00223 void MatrixSymPosDefBandedChol::Vp_StMtV(
00224   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00225   , const SpVectorSlice& x, value_type b) const
00226 {
00227   assert_initialized();
00228   MatrixOp::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00229 }
00230 
00231 void MatrixSymPosDefBandedChol::Vp_StPtMtV(
00232   DVectorSlice* y, value_type a
00233   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00234   , BLAS_Cpp::Transp M_trans
00235   , const DVectorSlice& x, value_type b) const
00236 {
00237   assert_initialized();
00238   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00239 }
00240 
00241 void MatrixSymPosDefBandedChol::Vp_StPtMtV(
00242   DVectorSlice* y, value_type a
00243   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00244   , BLAS_Cpp::Transp M_trans
00245   , const SpVectorSlice& x, value_type b) const
00246 {
00247   assert_initialized();
00248   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00249 }
00250 
00251 // Overridden from MatrixFactorized
00252 
00253 void MatrixSymPosDefBandedChol::V_InvMtV(
00254   DVectorSlice* y, BLAS_Cpp::Transp M_trans
00255   , const DVectorSlice& x) const
00256 {
00257   using Teuchos::Workspace;
00258   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00259 
00260   assert_initialized();
00261 
00262   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, BLAS_Cpp::no_trans, x.size() );
00263 
00264   update_factorization();
00265 
00266   Workspace<value_type>  t_ws(wss,y->size());
00267   DVectorSlice                 t(&t_ws[0],t_ws.size());
00268 
00269   t = x;
00270   
00271   FortranTypes::f_int info = 0;
00272   FORTRAN_FUNC_CALL_UL(DPBTRS,dpbtrs)(
00273     FORTRAN_CHAR_1_ARG_CALL(UB_uplo_ == BLAS_Cpp::upper ? 'U' : 'L')
00274     , n_, kd_, 1, UB_.col_ptr(1), UB_.max_rows()
00275     ,&t[0], t.size(), &info );
00276   if( info > 0 )
00277     TEUCHOS_TEST_FOR_EXCEPT(true); // Should not happen!
00278   if( info < 0 ) {
00279     std::ostringstream omsg;
00280     omsg
00281       << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00282       << "The " << -info << " argument passed to xPBTRF(...) is invalid!";
00283     throw std::invalid_argument(omsg.str());
00284   }
00285   *y = t;
00286 }
00287 
00288 // Private member functions
00289 
00290 void MatrixSymPosDefBandedChol::assert_initialized() const
00291 {
00292   if( n_ == 0 )
00293     throw std::logic_error("MatrixSymPosDefBandedChol::assert_initialized(): Error, "
00294                  "not initialized!" );
00295 }
00296 
00297 void MatrixSymPosDefBandedChol::update_factorization() const
00298 {
00299   namespace rcp = MemMngPack;
00300   using Teuchos::RCP;
00301   namespace rmp = MemMngPack;
00302 
00303   if( !factor_updated_ ) {
00304     if(UB_.rows() == 0) {
00305       // Allocate our own storage for the banded cholesky factor!
00306       typedef Teuchos::RCP<DMatrix>                  UB_ptr_t;
00307       typedef rmp::ReleaseResource_ref_count_ptr<DMatrix>  UB_rel_ptr_t;
00308       typedef Teuchos::RCP<UB_rel_ptr_t>               UB_rel_ptr_ptr_t;
00309       UB_rel_ptr_ptr_t  UB_rel_ptr_ptr = new UB_rel_ptr_t(new DMatrix);
00310       UB_rel_ptr_ptr->ptr->resize(kd_+1,n_);
00311       UB_.bind( (*UB_rel_ptr_ptr->ptr)() );
00312       UB_release_resource_ptr_ = Teuchos::rcp_implicit_cast<rmp::ReleaseResource>(UB_rel_ptr_ptr);
00313     }
00314     // Update the cholesky factor
00315     LinAlgOpPack::M_StM( &UB_, scale_, MB_, BLAS_Cpp::no_trans );
00316     UB_uplo_ = MB_uplo_;
00317     FortranTypes::f_int info = 0;
00318     FORTRAN_FUNC_CALL_UL(DPBTRF,dpbtrf)(
00319       FORTRAN_CHAR_1_ARG_CALL(UB_uplo_ == BLAS_Cpp::upper ? 'U' : 'L')
00320       , n_, kd_, UB_.col_ptr(1), UB_.max_rows(), &info );
00321     if( info < 0 ) {
00322       std::ostringstream omsg;
00323       omsg
00324         << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00325         << "The " << -info << " argument passed to xPBTRF(...) is invalid!";
00326       throw std::invalid_argument(omsg.str());
00327     }
00328     if( info > 0 ) {
00329       std::ostringstream omsg;
00330       omsg
00331         << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00332         << "The leading minor of order " << info << " passed to xPBTRF(...) is not positive definite!";
00333       throw std::invalid_argument(omsg.str());
00334     }
00335     factor_updated_ = true;
00336   }
00337 }
00338 
00339 } // end namespace ConstrainedOptPack
00340 
00341 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines