AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_dense_Vp_StPtMtV.hpp
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 #ifndef DENSE_V_P_S_T_P_T_M_T_V_H
00043 #define DENSE_V_P_S_T_P_T_M_T_V_H
00044 
00045 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00046 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00047 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_EtaVector.hpp"
00048 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00049 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00050 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00051 #include "DenseLinAlgPack_DMatrixClass.hpp"
00052 #include "DenseLinAlgPack_DMatrixOut.hpp"
00053 #include "DenseLinAlgPack_AssertOp.hpp"
00054 #include "MiWorkspacePack.h"
00055 
00056 namespace AbstractLinAlgPack {
00057 
00064 template<class M_t, class V_t>
00065 void dense_Vp_StPtMtV(
00066   DVectorSlice                 *y
00067   ,value_type                 a
00068   ,const GenPermMatrixSlice   &P
00069   ,BLAS_Cpp::Transp           P_trans
00070   ,const M_t                  &M
00071   ,BLAS_Cpp::Transp           M_trans
00072   ,const V_t                  &x
00073   ,value_type                 b
00074   )
00075 {
00076   using BLAS_Cpp::no_trans;
00077   using BLAS_Cpp::trans;
00078   using BLAS_Cpp::trans_not;
00079   using BLAS_Cpp::rows;
00080   using BLAS_Cpp::cols;
00081   using DenseLinAlgPack::dot;
00082   using DenseLinAlgPack::DVector;
00083   using DenseLinAlgPack::Vt_S;
00084   using AbstractLinAlgPack::dot;
00085   using AbstractLinAlgPack::Vp_StMtV;
00086   using AbstractLinAlgPack::GenPermMatrixSlice;
00087   typedef AbstractLinAlgPack::EtaVector eta;
00088   using Teuchos::Workspace;
00089   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00090   
00091   // Validate the sizes
00092   // 
00093   // y += op(P)*op(M)*x
00094   // 
00095   const DenseLinAlgPack::size_type
00096     ny = y->size(),
00097     nx = x.size(),
00098     opM_rows = rows( M.rows(), M.cols(), M_trans ),
00099     opM_cols = cols( M.rows(), M.cols(), M_trans );
00100   if(    ny != rows( P.rows(), P.cols(), P_trans )
00101        || nx != opM_cols
00102        || cols( P.rows(), P.cols(), P_trans ) != opM_rows )
00103     throw std::length_error( "MatrixOp::Vp_StPtMtV(...) : Error, "
00104       "sizes of arguments does not match up" );
00105   //
00106   // Compute y = b*y + a*op(P)*op(M)*x in a resonably efficient manner.  This
00107   // implementation will assume that M is stored as a dense matrix.  Either
00108   // t = op(M)*x is computed first (O(opM_rows*nx) flops) then y = b*y + a*op(P)*t
00109   // (O(ny) + O(P_nz) flops) or each row of t' = e(j)' * op(M) (O(nx) flops)
00110   // is computed one at a time and then y(i) = b*y(i) + a*t'*x (O(nx)  flops)
00111   // where op(P)(i,j) = 1.0.  In the latter case, there are P_nz rows
00112   // of op(M) that have to be generated so the total cost is O(P_nz*nx).
00113   // Therefore, we will do the former if opM_rows < P_nz and the latter otherwise.
00114   //
00115   if( !P.nz() ) {
00116     // y = b*y
00117     if(b==0.0)       *y = 0.0;
00118     else if(b!=1.0)  Vt_S(y,b);
00119   }
00120   else if( opM_rows > P.nz() || P.is_identity() ) {
00121     // t = op(M)*x
00122     Workspace<value_type> t_ws(wss,opM_rows);
00123     DVectorSlice t(&t_ws[0],t_ws.size());
00124     LinAlgOpPack::V_MtV( &t, M, M_trans, x );
00125     // y = b*y + a*op(P)*t
00126     Vp_StMtV( y, a, P, P_trans, t(), b );
00127   }
00128   else {
00129     // y = b*y
00130     if(b==0.0)       *y = 0.0;
00131     else if(b!=1.0)  Vt_S(y,b);
00132     // Compute t' = e(j)' * op(M) then y(i) += a*t'*x where op(P)(i,j) = 1.0
00133     Workspace<value_type> t_ws(wss,opM_cols);
00134     DVectorSlice t(&t_ws[0],t_ws.size());
00135     if( P.is_identity() ) {
00136       for( size_type k = 1; k <= P.nz(); ++k ) {
00137         const size_type
00138           i = k,
00139           j = k;
00140         // t = op(M') * e(j)      
00141         LinAlgOpPack::V_MtV( &t, M, trans_not(M_trans), eta(j,opM_rows)() );
00142         // y(i) += a*t'*x
00143         (*y)(i) += a * dot( t(), x );
00144       }
00145     }
00146     else {
00147       for( GenPermMatrixSlice::const_iterator itr = P.begin(); itr != P.end(); ++itr ) {
00148         const DenseLinAlgPack::size_type
00149           i = P_trans == no_trans ? itr->row_i() : itr->col_j(),
00150           j = P_trans == no_trans ? itr->col_j() : itr->row_i();
00151         // t = op(M') * e(j)      
00152         LinAlgOpPack::V_MtV( &t, M, trans_not(M_trans), eta(j,opM_rows)() );
00153         // y(i) += a*t'*x
00154         (*y)(i) += a * dot( t(), x );
00155       }
00156     }
00157   }
00158 }
00159 
00160 } // end namespace AbstractLinAlgPack
00161 
00162 #endif // DENSE_V_P_S_T_P_T_M_T_V_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends