MOOCHO (Single Doxygen Collection) Version of the Day
AbstractLinAlgPack_LinAlgOpPackHack.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 
00043 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00044 #include "AbstractLinAlgPack_VectorMutableDense.hpp"
00045 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00046 #include "AbstractLinAlgPack_MatrixOpGetGMS.hpp"
00047 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00048 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00049 #include "AbstractLinAlgPack_VectorMutable.hpp"
00050 #include "AbstractLinAlgPack_VectorSpace.hpp"
00051 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00052 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00053 #include "DenseLinAlgPack_DMatrixOp.hpp"
00054 
00055 
00056 void LinAlgOpPack::assign(
00057   DMatrixSlice* gms_lhs, const MatrixOp& M_rhs,
00058   BLAS_Cpp::Transp trans_rhs
00059   )
00060 {
00061   Mp_M_assert_sizes(
00062     gms_lhs->rows(), gms_lhs->cols(), BLAS_Cpp::no_trans,
00063     M_rhs.rows(), M_rhs.cols(), trans_rhs );
00064   (*gms_lhs) = 0.0;
00065   Mp_StM(gms_lhs,1.0,M_rhs,trans_rhs);
00066 }
00067 
00068 
00069 void LinAlgOpPack::Mp_StM(
00070   DMatrixSlice* C, value_type a
00071   ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
00072   )
00073 {
00074   using AbstractLinAlgPack::VectorSpace;
00075   using AbstractLinAlgPack::VectorDenseEncap;
00076   using AbstractLinAlgPack::MatrixOpGetGMS;
00077   using AbstractLinAlgPack::MatrixDenseEncap;
00078   const MatrixOpGetGMS
00079     *B_get_gms = dynamic_cast<const MatrixOpGetGMS*>(&B);
00080   if(B_get_gms) {
00081     DenseLinAlgPack::Mp_StM( C, a, MatrixDenseEncap(*B_get_gms)(), B_trans );   
00082   }
00083   else {
00084     const size_type num_cols = C->cols();
00085     VectorSpace::multi_vec_mut_ptr_t
00086       B_mv = ( B_trans == BLAS_Cpp::no_trans 
00087           ? B.space_cols() : B.space_rows()
00088           ).create_members(num_cols);
00089     assign(B_mv.get(),B,B_trans);
00090     for( size_type j = 1; j <= num_cols; ++j ) {
00091       DenseLinAlgPack::Vp_StV(&C->col(j),a,VectorDenseEncap(*B_mv->col(j))());
00092     }
00093   }
00094 }
00095 
00096 void LinAlgOpPack::Vp_StMtV(
00097   DVectorSlice* y, value_type a, const MatrixOp& M
00098   ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x, value_type b
00099   )
00100 {
00101   using BLAS_Cpp::no_trans;
00102   using AbstractLinAlgPack::VectorDenseMutableEncap;
00103   VectorSpace::vec_mut_ptr_t
00104     ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member(),
00105     ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
00106   (VectorDenseMutableEncap(*ay))() = *y;
00107   (VectorDenseMutableEncap(*ax))() = x;
00108   AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, *ax, b );
00109   *y = VectorDenseMutableEncap(*ay)();
00110 }
00111 
00112 void LinAlgOpPack::Vp_StMtV(
00113   DVectorSlice* y, value_type a, const MatrixOp& M
00114   ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x, value_type b
00115   )
00116 {
00117   using BLAS_Cpp::no_trans;
00118   using AbstractLinAlgPack::VectorDenseMutableEncap;
00119   VectorSpace::vec_mut_ptr_t
00120     ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member();
00121   (VectorDenseMutableEncap(*ay))() = *y;
00122   AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, x, b );
00123   *y = VectorDenseMutableEncap(*ay)();
00124 }
00125 
00126 void LinAlgOpPack::V_InvMtV(
00127   DVectorSlice* y, const MatrixOpNonsing& M
00128   ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00129   )
00130 {
00131   using BLAS_Cpp::trans;
00132   using AbstractLinAlgPack::VectorDenseMutableEncap;
00133   VectorSpace::vec_mut_ptr_t
00134     ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
00135     ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
00136   (VectorDenseMutableEncap(*ay))() = *y;
00137   (VectorDenseMutableEncap(*ax))() = x;
00138   AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
00139   *y = VectorDenseMutableEncap(*ay)();
00140 }
00141 
00142 void LinAlgOpPack::V_InvMtV(
00143   DVector* y, const MatrixOpNonsing& M
00144   ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00145   )
00146 {
00147   using BLAS_Cpp::trans;
00148   using AbstractLinAlgPack::VectorDenseMutableEncap;
00149   VectorSpace::vec_mut_ptr_t
00150     ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
00151     ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
00152   (VectorDenseMutableEncap(*ax))() = x;
00153   AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
00154   *y = VectorDenseMutableEncap(*ay)();
00155 }
00156 
00157 void LinAlgOpPack::V_InvMtV(
00158   DVectorSlice* y, const MatrixOpNonsing& M
00159   ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00160   )
00161 {
00162   using BLAS_Cpp::trans;
00163   using AbstractLinAlgPack::VectorDenseMutableEncap;
00164   VectorSpace::vec_mut_ptr_t
00165     ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member();
00166   AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, x );
00167   *y = VectorDenseMutableEncap(*ay)();
00168 }
00169 
00170 void LinAlgOpPack::V_InvMtV(
00171   DVector* y, const MatrixOpNonsing& M
00172   ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00173   )
00174 {
00175   y->resize(M.rows());
00176   LinAlgOpPack::V_InvMtV( &(*y)(), M, M_trans, x );
00177 }
00178 
00179 // These methods below are a real problem to implement in general.
00180 //
00181 // If the column space of op(M) could not return the above vector space
00182 // for op(M).space_cols().space(P,P_trans) then we will try, as a last
00183 // resort, to a dense serial vector and see what happens.
00184 
00185 void LinAlgOpPack::Vp_StPtMtV(
00186   DVectorSlice* y, value_type a
00187   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00188   ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
00189   ,const DVectorSlice& x, value_type b
00190   )
00191 {
00192   namespace mmp = MemMngPack;
00193   using BLAS_Cpp::no_trans;
00194   using AbstractLinAlgPack::VectorMutableDense;
00195   using AbstractLinAlgPack::VectorDenseMutableEncap;
00196   using AbstractLinAlgPack::Vp_StPtMtV;
00197   VectorSpace::space_ptr_t
00198     ay_space = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans);
00199   VectorSpace::vec_mut_ptr_t
00200     ay =  ( ay_space.get()
00201         ? ay_space->create_member()
00202         : Teuchos::rcp_implicit_cast<VectorMutable>(
00203           Teuchos::rcp(new VectorMutableDense(BLAS_Cpp::rows(P.rows(),P.cols(),P_trans)))
00204           ) ),
00205     ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
00206   (VectorDenseMutableEncap(*ay))() = *y;
00207   (VectorDenseMutableEncap(*ax))() = x;
00208   Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, *ax, b );
00209   *y = VectorDenseMutableEncap(*ay)();
00210 }
00211 
00212 void LinAlgOpPack::Vp_StPtMtV(
00213   DVectorSlice* y, value_type a
00214   ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00215   ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
00216   ,const SpVectorSlice& x, value_type b
00217   )
00218 {
00219   using BLAS_Cpp::no_trans;
00220   using AbstractLinAlgPack::VectorMutableDense;
00221   using AbstractLinAlgPack::VectorDenseMutableEncap;
00222   using AbstractLinAlgPack::Vp_StPtMtV;
00223   VectorSpace::vec_mut_ptr_t
00224     ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans)->create_member();
00225   (VectorDenseMutableEncap(*ay))() = *y;
00226   Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, x, b );
00227   *y = VectorDenseMutableEncap(*ay)();
00228 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines