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

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