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 // 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 <iostream> // Debuggin only
00030 
00031 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00032 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00033 #include "AbstractLinAlgPack_VectorMutable.hpp"
00034 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00035 #include "AbstractLinAlgPack_SpVectorClass.hpp"
00036 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00037 #include "Teuchos_TestForException.hpp"
00038 
00039 namespace AbstractLinAlgPack {
00040 
00041 MatrixSymDiagStd::MatrixSymDiagStd(
00042   const VectorSpace::vec_mut_ptr_t& diag
00043   ,bool                             unique
00044   )
00045 {
00046   this->initialize(diag,unique);
00047 //  std::cerr << "MatrixSymDiagStd::rows() = " << this->rows() << std::endl; // Debugging
00048 //  std::cerr << "MatrixSymDiagStd::nz() = "   << this->nz()   << std::endl; // Debugging
00049 //  std::cerr << "MatrixSymDiagStd::cols() = " << this->cols() << std::endl; // Debugging
00050 //  std::cerr << "MatrixSymDiagStd::nz() = "   << this->nz()   << std::endl; // Debugging
00051 }
00052 
00053 void MatrixSymDiagStd::initialize(
00054   const VectorSpace::vec_mut_ptr_t& diag
00055   ,bool                             unique
00056   )
00057 {
00058   diag_   = diag;   // lazy copy!
00059   unique_ = unique;
00060 }
00061 
00062 VectorMutable& MatrixSymDiagStd::diag()
00063 {
00064   copy_unique();
00065   VectorMutable *diag = diag_.get();
00066   TEST_FOR_EXCEPTION(
00067     !diag, std::logic_error
00068     ,"MatrixSymDiagStd::diag(): Error, the diagonal vector has not been set! " );
00069   return *diag;;
00070 }
00071 
00072 const VectorSpace::vec_mut_ptr_t&
00073 MatrixSymDiagStd::diag_ptr() const
00074 {
00075   return diag_;
00076 }
00077 
00078 // Overridden from MatrixBase
00079 
00080 size_type MatrixSymDiagStd::rows() const
00081 {
00082   return diag_.get() ? diag_->dim() : 0;
00083 }
00084 
00085 size_type MatrixSymDiagStd::nz() const
00086 {
00087   return diag_.get() ? diag_->nz() : 0;
00088 }
00089 
00090 // Overridden from MatrixOp
00091 
00092 const VectorSpace& MatrixSymDiagStd::space_cols() const {
00093   return diag_->space();
00094 }
00095 
00096 const VectorSpace& MatrixSymDiagStd::space_rows() const {
00097   return diag_->space();
00098 }
00099 
00100 MatrixOp&
00101 MatrixSymDiagStd::operator=(const MatrixOp& M)
00102 {
00103   const MatrixSymDiagStd
00104     *p_M = dynamic_cast<const MatrixSymDiagStd*>(&M);
00105 
00106   TEST_FOR_EXCEPTION(
00107     p_M == NULL, std::logic_error
00108     ,"MatrixSymDiagStd::operator=(M): Error, the matrix M with concrete type "
00109     "\'" << typeName(M) << "\' does not support the MatrixSymDiagStd type! " );
00110 
00111   if( p_M == this ) return *this; // Assignment to self
00112 
00113   diag_    = p_M->diag_;  // lazy copy!
00114   unique_  = p_M->unique_;
00115 
00116   return *this;
00117 }
00118 
00119 bool MatrixSymDiagStd::Mp_StM(
00120   MatrixOp* M_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs) const
00121 {
00122   // ToDo: validate the vector spaces for the matrices!
00123   MultiVectorMutable
00124     *M_mv_lhs = dynamic_cast<MultiVectorMutable*>(M_lhs);
00125   if(!M_mv_lhs)
00126     return false;
00127   VectorSpace::vec_mut_ptr_t
00128     M_diag = M_mv_lhs->diag(0);
00129   if(!M_diag.get())
00130     return false; // Access to the diagonal is not supported!
00131   Vp_StV( M_diag.get(), alpha, *diag_ );
00132   return true;
00133 }
00134 
00135 void MatrixSymDiagStd::Vp_StMtV(
00136   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00137   , const Vector& v_rhs2, value_type beta) const
00138 {
00139   // ToDo: Validate input!
00140   if(beta == 0.0)
00141     *v_lhs = 0.0;
00142   else if(beta != 1.0)
00143     Vt_S( v_lhs, beta );
00144   ele_wise_prod( alpha, v_rhs2, *diag_, v_lhs );
00145 }
00146 
00147 void MatrixSymDiagStd::Vp_StMtV(
00148   VectorMutable* v_lhs, value_type alpha, BLAS_Cpp::Transp trans_rhs1
00149   , const SpVectorSlice& sv_rhs2, value_type beta) const
00150 {
00151   MatrixOp::Vp_StMtV(v_lhs,alpha,trans_rhs1,sv_rhs2,beta); // ToDo: Implement specialized!
00152 }
00153 
00154 // Overridden from MatrixNonsing
00155 
00156 void MatrixSymDiagStd::V_InvMtV(
00157   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00158   , const Vector& v_rhs2) const
00159 {
00160   ele_wise_divide( 1.0, v_rhs2, *diag_, v_lhs );
00161 }
00162 
00163 void MatrixSymDiagStd::V_InvMtV(
00164   VectorMutable* v_lhs, BLAS_Cpp::Transp trans_rhs1
00165   , const SpVectorSlice& sv_rhs2) const
00166 {
00167   MatrixNonsing::V_InvMtV(v_lhs,trans_rhs1,sv_rhs2 ); // ToDo: Implement specialized!
00168 }
00169 
00170 bool MatrixSymDiagStd::syrk(
00171   BLAS_Cpp::Transp   A_trans
00172   ,value_type        a
00173   ,value_type        b
00174   ,MatrixSymOp   *B
00175   ) const
00176 {
00177   MatrixSymDiagStd    *B_sd = dynamic_cast<MatrixSymDiagStd*>(B);
00178   if(!B_sd) return false;
00179   VectorMutable     &B_diag = B_sd->diag();
00180   const Vector      &A_diag = this->diag();
00181   // B = b*B + a*A*A
00182   Vt_S( &B_diag, b );
00183   ele_wise_prod( 1.0, A_diag, A_diag, &B_diag );   // B.diag(i) += a * (A.diag)(i) * (A.diag)(i)
00184   return true;
00185 }
00186 
00187 // Overridden from MatrixSymInitDiag
00188 
00189 void MatrixSymDiagStd::init_identity( const VectorSpace& space_diag, value_type alpha )
00190 {
00191   diag_ = space_diag.create_member();
00192   if( diag_->dim() )
00193     *diag_ = alpha;
00194 }
00195 
00196 void MatrixSymDiagStd::init_diagonal( const Vector& diag )
00197 {
00198   diag_ = diag.space().create_member();
00199   *diag_ = diag;
00200 }
00201 
00202 // Overridden from MatrixSymDiag
00203 
00204 const Vector& MatrixSymDiagStd::diag() const
00205 {
00206   return const_cast<MatrixSymDiagStd*>(this)->diag();
00207 }
00208 
00209 // private
00210 
00211 void MatrixSymDiagStd::copy_unique()
00212 {
00213   if( diag_.get() && diag_.count() > 1 && unique_ )
00214     diag_ = diag_->clone();
00215 }
00216 
00217 } // end namespace AbstractLinAlgPack

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