MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MatrixSymDiagStd.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <iostream> // Debuggin only
00043 
00044 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00045 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00046 #include "AbstractLinAlgPack_VectorMutable.hpp"
00047 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00048 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00049 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 
00052 namespace AbstractLinAlgPack {
00053 
00054 MatrixSymDiagStd::MatrixSymDiagStd(
00055   const VectorSpace::vec_mut_ptr_t& diag
00056   ,bool                             unique
00057   )
00058 {
00059   this->initialize(diag,unique);
00060 //  std::cerr << "MatrixSymDiagStd::rows() = " << this->rows() << std::endl; // Debugging
00061 //  std::cerr << "MatrixSymDiagStd::nz() = "   << this->nz()   << std::endl; // Debugging
00062 //  std::cerr << "MatrixSymDiagStd::cols() = " << this->cols() << std::endl; // Debugging
00063 //  std::cerr << "MatrixSymDiagStd::nz() = "   << this->nz()   << std::endl; // Debugging
00064 }
00065 
00066 void MatrixSymDiagStd::initialize(
00067   const VectorSpace::vec_mut_ptr_t& diag
00068   ,bool                             unique
00069   )
00070 {
00071   diag_   = diag;   // lazy copy!
00072   unique_ = unique;
00073 }
00074 
00075 VectorMutable& MatrixSymDiagStd::diag()
00076 {
00077   copy_unique();
00078   VectorMutable *diag = diag_.get();
00079   TEUCHOS_TEST_FOR_EXCEPTION(
00080     !diag, std::logic_error
00081     ,"MatrixSymDiagStd::diag(): Error, the diagonal vector has not been set! " );
00082   return *diag;;
00083 }
00084 
00085 const VectorSpace::vec_mut_ptr_t&
00086 MatrixSymDiagStd::diag_ptr() const
00087 {
00088   return diag_;
00089 }
00090 
00091 // Overridden from MatrixBase
00092 
00093 size_type MatrixSymDiagStd::rows() const
00094 {
00095   return diag_.get() ? diag_->dim() : 0;
00096 }
00097 
00098 size_type MatrixSymDiagStd::nz() const
00099 {
00100   return diag_.get() ? diag_->nz() : 0;
00101 }
00102 
00103 // Overridden from MatrixOp
00104 
00105 const VectorSpace& MatrixSymDiagStd::space_cols() const {
00106   return diag_->space();
00107 }
00108 
00109 const VectorSpace& MatrixSymDiagStd::space_rows() const {
00110   return diag_->space();
00111 }
00112 
00113 MatrixOp&
00114 MatrixSymDiagStd::operator=(const MatrixOp& M)
00115 {
00116   const MatrixSymDiagStd
00117     *p_M = dynamic_cast<const MatrixSymDiagStd*>(&M);
00118 
00119   TEUCHOS_TEST_FOR_EXCEPTION(
00120     p_M == NULL, std::logic_error
00121     ,"MatrixSymDiagStd::operator=(M): Error, the matrix M with concrete type "
00122     "\'" << typeName(M) << "\' does not support the MatrixSymDiagStd type! " );
00123 
00124   if( p_M == this ) return *this; // Assignment to self
00125 
00126   diag_    = p_M->diag_;  // lazy copy!
00127   unique_  = p_M->unique_;
00128 
00129   return *this;
00130 }
00131 
00132 bool MatrixSymDiagStd::Mp_StM(
00133   MatrixOp* M_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
00134 {
00135   // ToDo: validate the vector spaces for the matrices!
00136   MultiVectorMutable
00137     *M_mv_lhs = dynamic_cast<MultiVectorMutable*>(M_lhs);
00138   if(!M_mv_lhs)
00139     return false;
00140   VectorSpace::vec_mut_ptr_t
00141     M_diag = M_mv_lhs->diag(0);
00142   if(!M_diag.get())
00143     return false; // Access to the diagonal is not supported!
00144   Vp_StV( M_diag.get(), alpha, *diag_ );
00145   return true;
00146 }
00147 
00148 void MatrixSymDiagStd::Vp_StMtV(
00149   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00150   , const Vector& v_rhs2, value_type beta) const
00151 {
00152   // ToDo: Validate input!
00153   if(beta == 0.0)
00154     *v_lhs = 0.0;
00155   else if(beta != 1.0)
00156     Vt_S( v_lhs, beta );
00157   ele_wise_prod( alpha, v_rhs2, *diag_, v_lhs );
00158 }
00159 
00160 void MatrixSymDiagStd::Vp_StMtV(
00161   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00162   , const SpVectorSlice& sv_rhs2, value_type beta) const
00163 {
00164   MatrixOp::Vp_StMtV(v_lhs,alpha,trans_rhs1,sv_rhs2,beta); // ToDo: Implement specialized!
00165 }
00166 
00167 // Overridden from MatrixNonsing
00168 
00169 void MatrixSymDiagStd::V_InvMtV(
00170   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00171   , const Vector& v_rhs2) const
00172 {
00173   ele_wise_divide( 1.0, v_rhs2, *diag_, v_lhs );
00174 }
00175 
00176 void MatrixSymDiagStd::V_InvMtV(
00177   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00178   , const SpVectorSlice& sv_rhs2) const
00179 {
00180   MatrixNonsing::V_InvMtV(v_lhs,trans_rhs1,sv_rhs2 ); // ToDo: Implement specialized!
00181 }
00182 
00183 bool MatrixSymDiagStd::syrk(
00184   BLAS_Cpp::Transp   A_trans
00185   ,value_type        a
00186   ,value_type        b
00187   ,MatrixSymOp   *B
00188   ) const
00189 {
00190   MatrixSymDiagStd    *B_sd = dynamic_cast<MatrixSymDiagStd*>(B);
00191   if(!B_sd) return false;
00192   VectorMutable     &B_diag = B_sd->diag();
00193   const Vector      &A_diag = this->diag();
00194   // B = b*B + a*A*A
00195   Vt_S( &B_diag, b );
00196   ele_wise_prod( 1.0, A_diag, A_diag, &B_diag );   // B.diag(i) += a * (A.diag)(i) * (A.diag)(i)
00197   return true;
00198 }
00199 
00200 // Overridden from MatrixSymInitDiag
00201 
00202 void MatrixSymDiagStd::init_identity( const VectorSpace& space_diag, value_type alpha )
00203 {
00204   diag_ = space_diag.create_member();
00205   if( diag_->dim() )
00206     *diag_ = alpha;
00207 }
00208 
00209 void MatrixSymDiagStd::init_diagonal( const Vector& diag )
00210 {
00211   diag_ = diag.space().create_member();
00212   *diag_ = diag;
00213 }
00214 
00215 // Overridden from MatrixSymDiag
00216 
00217 const Vector& MatrixSymDiagStd::diag() const
00218 {
00219   return const_cast<MatrixSymDiagStd*>(this)->diag();
00220 }
00221 
00222 // private
00223 
00224 void MatrixSymDiagStd::copy_unique()
00225 {
00226   if( diag_.get() && diag_.total_count() > 1 && unique_ )
00227     diag_ = diag_->clone();
00228 }
00229 
00230 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines