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

Generated on Thu Sep 18 12:35:14 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1