MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_MatrixSymOpSerial.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 <assert.h>
00043 
00044 #include "AbstractLinAlgPack_MatrixSymOpSerial.hpp"
00045 #include "AbstractLinAlgPack_MatrixSymOpGetGMSSymMutable.hpp"
00046 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00047 #include "AbstractLinAlgPack_EtaVector.hpp"
00048 #include "DenseLinAlgPack_DMatrixOp.hpp"
00049 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00050 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00051 #include "DenseLinAlgPack_AssertOp.hpp"
00052 #include "Teuchos_dyn_cast.hpp"
00053 
00054 namespace AbstractLinAlgPack {
00055 
00056 void MatrixSymOpSerial::Mp_StPtMtP(
00057   DMatrixSliceSym* S, value_type a
00058   ,EMatRhsPlaceHolder
00059   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00060   ,value_type b
00061   ) const
00062 {
00063   using BLAS_Cpp::no_trans;
00064   using BLAS_Cpp::trans;
00065   using BLAS_Cpp::trans_not;
00066   using BLAS_Cpp::cols;
00067   //
00068   // S = b*S
00069   //
00070   // S += a*op(P')*M*op(P)
00071   //
00072   // We will perform this operation for each column in op(P) as:
00073   //
00074   // op(P)(:,j(k)) = e(i(k)) <: R^n
00075   //
00076   // S += a*op(P')*M*[ ... e(i(1)) ... e(i(k)) ... e(i(nz)) ... ]
00077   //                                     j(k)
00078   //
00079   // We will perform this by column as:
00080   //
00081   // for k = 1...nz
00082   //    S(:,j(k)) += a*y_k
00083   //
00084   //    where:
00085   //       y_k = op(P')*M*e(i(k))
00086   //
00087   // Above we only need to set the portion of S(:,j(k)) for the stored part
00088   // of the symmetric matrix (i.e. upper part for upper and lower part for lower).
00089   //
00090   DenseLinAlgPack::MtM_assert_sizes(
00091     this->rows(), this->cols(), no_trans
00092     , P.rows(), P.cols(), P_trans );
00093   DenseLinAlgPack::Mp_M_assert_sizes(
00094     S->rows(), S->cols(), no_trans
00095     , cols( P.rows(), P.cols(), P_trans )
00096     , cols( P.rows(), P.cols(), P_trans )
00097     , no_trans );
00098   //
00099   const size_type
00100     n = this->rows(),
00101     m = S->rows();
00102   // S = b*S
00103   if( b != 1.0 )
00104     DenseLinAlgPack::Mt_S( &DenseLinAlgPack::nonconst_tri_ele(S->gms(),S->uplo()), b );
00105   // Set the colums of S
00106   DVector y_k_store(m);
00107   DVectorSlice y_k = y_k_store();
00108   for( GenPermMatrixSlice::const_iterator P_itr = P.begin(); P_itr != P.end(); ++P_itr )
00109   {
00110     const size_type
00111       i_k = P_trans == no_trans ? P_itr->row_i() : P_itr->col_j(),
00112       j_k = P_trans == no_trans ? P_itr->col_j() : P_itr->row_i();
00113     // e(i(k))
00114     EtaVector
00115       e_i_k(i_k,n);
00116     // y_k = op(P')*M*e(i(k))
00117     AbstractLinAlgPack::Vp_StPtMtV( &y_k, 1.0, P, trans_not(P_trans), *this, no_trans, e_i_k(), 0.0 );
00118     // S(:,j(k)) += a*y_k
00119     if( S->uplo() == BLAS_Cpp::upper )
00120       DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(1,j_k), a, y_k(1,j_k) );
00121     else
00122       DenseLinAlgPack::Vp_StV( &S->gms().col(j_k)(j_k,m), a, y_k(j_k,m) );
00123   }
00124 }
00125 
00126 void MatrixSymOpSerial::Mp_StMtMtM(
00127   DMatrixSliceSym* sym_lhs, value_type alpha
00128   ,EMatRhsPlaceHolder dummy_place_holder
00129   ,const MatrixOpSerial& mwo_rhs, BLAS_Cpp::Transp mwo_rhs_trans
00130   ,value_type beta
00131   ) const
00132 {
00133   TEUCHOS_TEST_FOR_EXCEPT(true);  // ToDo: Give this some default implementation for 
00134         // this at some point in the future?
00135 }
00136 
00137 // Overridden from MatrixSymOp
00138 
00139 const VectorSpace& MatrixSymOpSerial::space_rows() const
00140 {
00141   return MatrixOpSerial::space_rows();
00142 }
00143 
00144 void MatrixSymOpSerial::Mp_StPtMtP(
00145   MatrixSymOp* symwo_lhs, value_type alpha
00146   ,EMatRhsPlaceHolder dummy
00147   ,const GenPermMatrixSlice& gpms_rhs, BLAS_Cpp::Transp gpms_rhs_trans
00148   ,value_type beta
00149   ) const
00150 {
00151   this->Mp_StPtMtP(
00152     &MatrixDenseSymMutableEncap(symwo_lhs)(), alpha, dummy
00153     ,gpms_rhs, gpms_rhs_trans
00154     , beta );
00155 }
00156 
00157 void MatrixSymOpSerial::Mp_StMtMtM(
00158   MatrixSymOp* symwo_lhs, value_type alpha
00159   ,EMatRhsPlaceHolder dummy
00160   ,const MatrixOp& mwo_rhs, BLAS_Cpp::Transp mwo_rhs_trans
00161   ,value_type beta
00162   ) const
00163 {
00164   using Teuchos::dyn_cast;
00165   this->Mp_StMtMtM(
00166     &MatrixDenseSymMutableEncap(symwo_lhs)(), alpha, dummy
00167     ,dyn_cast<const MatrixOpSerial>(mwo_rhs), mwo_rhs_trans
00168     ,beta );
00169 }
00170 
00171 } // end namespace AbstractLinAlgPack 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines