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 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00027 // 
00028 // ***********************************************************************
00029 // @HEADER
00030 
00031 #include <assert.h>
00032 
00033 #include <sstream>
00034 
00035 #include "ConstrainedOptPack_MatrixSymPosDefBandedChol.hpp"
00036 #include "DenseLinAlgPack_AssertOp.hpp"
00037 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00038 #include "DenseLinAlgPack_BLAS_Cpp.hpp"
00039 #include "MiReleaseResource_ref_count_ptr.h"
00040 #include "MiWorkspacePack.h"
00041 
00042 // LAPACK functions
00043 
00044 extern "C" {
00045 
00046 FORTRAN_FUNC_DECL_UL(void,DPBTRF,dpbtrf)(
00047   FORTRAN_CONST_CHAR_1_ARG(UPLO)
00048   ,const FortranTypes::f_int       &N
00049   ,const FortranTypes::f_int       &KD
00050   ,FortranTypes::f_dbl_prec        *AB
00051   ,const FortranTypes::f_int       &LDAB
00052   ,FortranTypes::f_int             *INFO
00053   );
00054 
00055 FORTRAN_FUNC_DECL_UL(void,DPBTRS,dpbtrs)(
00056   FORTRAN_CONST_CHAR_1_ARG(UPLO)
00057   ,const FortranTypes::f_int       &N
00058   ,const FortranTypes::f_int       &KD
00059   ,const FortranTypes::f_int       &NRHS
00060   ,const FortranTypes::f_dbl_prec  AB[]
00061   ,const FortranTypes::f_int       &LDAB
00062   ,FortranTypes::f_dbl_prec        *B
00063   ,const FortranTypes::f_int       &LDB
00064   ,FortranTypes::f_int             *INFO
00065   );
00066 
00067 } // end namespace LAPACK
00068 
00069 namespace ConstrainedOptPack {
00070 
00071 MatrixSymPosDefBandedChol::MatrixSymPosDefBandedChol(
00072   size_type                         n
00073   ,size_type                        kd
00074   ,DMatrixSlice                   *MB
00075   ,const release_resource_ptr_t&    MB_release_resource_ptr
00076   ,BLAS_Cpp::Uplo                   MB_uplo
00077   ,DMatrixSlice                   *UB
00078   ,const release_resource_ptr_t&    UB_release_resource_ptr
00079   ,BLAS_Cpp::Uplo                   UB_uplo
00080   ,bool                             update_factor
00081   ,value_type                       scale
00082   )
00083 {
00084   initialize(n,kd,MB,MB_release_resource_ptr,MB_uplo
00085          ,UB,UB_release_resource_ptr,UB_uplo,update_factor,scale);
00086 }
00087 
00088 void MatrixSymPosDefBandedChol::initialize(
00089   size_type                         n
00090   ,size_type                        kd
00091   ,DMatrixSlice                   *MB
00092   ,const release_resource_ptr_t&    MB_release_resource_ptr
00093   ,BLAS_Cpp::Uplo                   MB_uplo
00094   ,DMatrixSlice                   *UB
00095   ,const release_resource_ptr_t&    UB_release_resource_ptr
00096   ,BLAS_Cpp::Uplo                   UB_uplo
00097   ,bool                             update_factor
00098   ,value_type                       scale
00099   )
00100 {
00101   // Validate input
00102 
00103   if( n == 0 ) {
00104     if( kd != 0 )
00105       throw std::invalid_argument(
00106         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00107         "kd must be 0 if n == 0" );
00108     if( MB != NULL )
00109       throw std::invalid_argument(
00110         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00111         "MB must be NULL if n == 0" );
00112     if( MB_release_resource_ptr.get() != NULL )
00113       throw std::invalid_argument(
00114         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00115         "MB_release_resource_ptr.get() must be NULL if n == 0" );
00116     if( UB != NULL )
00117       throw std::invalid_argument(
00118         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00119         "UB must be NULL if n == 0" );
00120     if( UB_release_resource_ptr.get() != NULL )
00121       throw std::invalid_argument(
00122         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00123         "UB_release_resource_ptr.get() must be NULL if n == 0" );
00124   }
00125   else {
00126     if( kd + 1 > n )
00127       throw std::invalid_argument(
00128         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00129         "kd + 1 can not be larger than n" );
00130     if( MB == NULL && UB == NULL )
00131       throw std::invalid_argument(
00132         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00133         "MB and UB can not both be NULL" );
00134     if( MB != NULL && ( MB->rows() != kd + 1 || MB->cols() != n ) )
00135       throw std::invalid_argument(
00136         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00137         "MB is not the correct size" );
00138     if( UB != NULL && ( UB->rows() != kd + 1 || UB->cols() != n ) )
00139       throw std::invalid_argument(
00140         "MatrixSymPosDefBandedChol::initialize(...): Error, "
00141         "UB is not the correct size" );
00142   }
00143 
00144   // Set the members
00145 
00146   if( n == 0 ) {
00147     n_                        = 0;
00148     kd_                       = 0;
00149     MB_.bind(DMatrixSlice());
00150     MB_release_resource_ptr_  = NULL;
00151     MB_uplo_                  = BLAS_Cpp::lower;
00152     UB_.bind(DMatrixSlice());
00153     UB_release_resource_ptr_  = NULL;
00154     UB_uplo_                  = BLAS_Cpp::lower;
00155     scale_                    = 1.0;
00156   }
00157   else {
00158     // Set the members
00159     n_                        = n;
00160     kd_                       = kd;
00161     if(MB) MB_.bind(*MB);
00162     MB_release_resource_ptr_  = MB_release_resource_ptr;
00163     MB_uplo_                  = MB_uplo;
00164     if(UB) UB_.bind(*UB);
00165     UB_release_resource_ptr_  = UB_release_resource_ptr;
00166     UB_uplo_                  = UB_uplo;
00167     factor_updated_           = UB && !update_factor;
00168     scale_                    = scale;
00169     // Update the factorization if we have storage
00170     if( update_factor )
00171       update_factorization();
00172   }
00173 }
00174 
00175 // Overridden from MatrixOp
00176 
00177 size_type MatrixSymPosDefBandedChol::rows() const
00178 {
00179   return n_;
00180 }
00181 
00182 size_type MatrixSymPosDefBandedChol::nz() const
00183 {
00184   return (2 * kd_ + 1) * n_ - ( (kd_+1)*(kd_+1) - (kd_+1) );
00185 }
00186 
00187 std::ostream& MatrixSymPosDefBandedChol::output(std::ostream& out) const
00188 {
00189   return MatrixOp::output(out); // ToDo: Implement specialized version later!
00190 }
00191 
00192 void MatrixSymPosDefBandedChol::Vp_StMtV(
00193   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00194   , const DVectorSlice& x, value_type b) const
00195 {
00196   assert_initialized();
00197   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, BLAS_Cpp::no_trans, x.size() );
00198   if( MB_.rows() ) {
00199     BLAS_Cpp::sbmv(MB_uplo_,n_,kd_,a,MB_.col_ptr(1),MB_.max_rows(),x.raw_ptr(),x.stride()
00200              ,b,y->raw_ptr(),y->stride());
00201   }
00202   else if( UB_.rows() ) {
00203     TEST_FOR_EXCEPT(true); // ToDo: Implement when and if needed!
00204   }
00205   else {
00206     TEST_FOR_EXCEPT(true); // This should not happen!
00207   }
00208 }
00209 
00210 void MatrixSymPosDefBandedChol::Vp_StMtV(
00211   DVectorSlice* y, value_type a, BLAS_Cpp::Transp M_trans
00212   , const SpVectorSlice& x, value_type b) const
00213 {
00214   assert_initialized();
00215   MatrixOp::Vp_StMtV(y,a,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00216 }
00217 
00218 void MatrixSymPosDefBandedChol::Vp_StPtMtV(
00219   DVectorSlice* y, value_type a
00220   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00221   , BLAS_Cpp::Transp M_trans
00222   , const DVectorSlice& x, value_type b) const
00223 {
00224   assert_initialized();
00225   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00226 }
00227 
00228 void MatrixSymPosDefBandedChol::Vp_StPtMtV(
00229   DVectorSlice* y, value_type a
00230   , const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00231   , BLAS_Cpp::Transp M_trans
00232   , const SpVectorSlice& x, value_type b) const
00233 {
00234   assert_initialized();
00235   MatrixOp::Vp_StPtMtV(y,a,P,P_trans,M_trans,x,b); // ToDo: Implement spacialized operation when needed!
00236 }
00237 
00238 // Overridden from MatrixFactorized
00239 
00240 void MatrixSymPosDefBandedChol::V_InvMtV(
00241   DVectorSlice* y, BLAS_Cpp::Transp M_trans
00242   , const DVectorSlice& x) const
00243 {
00244   using Teuchos::Workspace;
00245   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00246 
00247   assert_initialized();
00248 
00249   DenseLinAlgPack::Vp_MtV_assert_sizes( y->size(), n_, n_, BLAS_Cpp::no_trans, x.size() );
00250 
00251   update_factorization();
00252 
00253   Workspace<value_type>  t_ws(wss,y->size());
00254   DVectorSlice                 t(&t_ws[0],t_ws.size());
00255 
00256   t = x;
00257   
00258   FortranTypes::f_int info = 0;
00259   FORTRAN_FUNC_CALL_UL(DPBTRS,dpbtrs)(
00260     FORTRAN_CHAR_1_ARG_CALL(UB_uplo_ == BLAS_Cpp::upper ? 'U' : 'L')
00261     , n_, kd_, 1, UB_.col_ptr(1), UB_.max_rows()
00262     ,&t[0], t.size(), &info );
00263   if( info > 0 )
00264     TEST_FOR_EXCEPT(true); // Should not happen!
00265   if( info < 0 ) {
00266     std::ostringstream omsg;
00267     omsg
00268       << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00269       << "The " << -info << " argument passed to xPBTRF(...) is invalid!";
00270     throw std::invalid_argument(omsg.str());
00271   }
00272   *y = t;
00273 }
00274 
00275 // Private member functions
00276 
00277 void MatrixSymPosDefBandedChol::assert_initialized() const
00278 {
00279   if( n_ == 0 )
00280     throw std::logic_error("MatrixSymPosDefBandedChol::assert_initialized(): Error, "
00281                  "not initialized!" );
00282 }
00283 
00284 void MatrixSymPosDefBandedChol::update_factorization() const
00285 {
00286   namespace rcp = MemMngPack;
00287   using Teuchos::RCP;
00288   namespace rmp = MemMngPack;
00289 
00290   if( !factor_updated_ ) {
00291     if(UB_.rows() == 0) {
00292       // Allocate our own storage for the banded cholesky factor!
00293       typedef Teuchos::RCP<DMatrix>                  UB_ptr_t;
00294       typedef rmp::ReleaseResource_ref_count_ptr<DMatrix>  UB_rel_ptr_t;
00295       typedef Teuchos::RCP<UB_rel_ptr_t>               UB_rel_ptr_ptr_t;
00296       UB_rel_ptr_ptr_t  UB_rel_ptr_ptr = new UB_rel_ptr_t(new DMatrix);
00297       UB_rel_ptr_ptr->ptr->resize(kd_+1,n_);
00298       UB_.bind( (*UB_rel_ptr_ptr->ptr)() );
00299       UB_release_resource_ptr_ = Teuchos::rcp_implicit_cast<rmp::ReleaseResource>(UB_rel_ptr_ptr);
00300     }
00301     // Update the cholesky factor
00302     LinAlgOpPack::M_StM( &UB_, scale_, MB_, BLAS_Cpp::no_trans );
00303     UB_uplo_ = MB_uplo_;
00304     FortranTypes::f_int info = 0;
00305     FORTRAN_FUNC_CALL_UL(DPBTRF,dpbtrf)(
00306       FORTRAN_CHAR_1_ARG_CALL(UB_uplo_ == BLAS_Cpp::upper ? 'U' : 'L')
00307       , n_, kd_, UB_.col_ptr(1), UB_.max_rows(), &info );
00308     if( info < 0 ) {
00309       std::ostringstream omsg;
00310       omsg
00311         << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00312         << "The " << -info << " argument passed to xPBTRF(...) is invalid!";
00313       throw std::invalid_argument(omsg.str());
00314     }
00315     if( info > 0 ) {
00316       std::ostringstream omsg;
00317       omsg
00318         << "MatrixSymPosDefBandedChol::update_factorization(): Error, "
00319         << "The leading minor of order " << info << " passed to xPBTRF(...) is not positive definite!";
00320       throw std::invalid_argument(omsg.str());
00321     }
00322     factor_updated_ = true;
00323   }
00324 }
00325 
00326 } // end namespace ConstrainedOptPack
00327 
00328 #endif // 0

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