00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00031 #include "AbstractLinAlgPack_VectorMutableDense.hpp"
00032 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00033 #include "AbstractLinAlgPack_MatrixOpGetGMS.hpp"
00034 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00035 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
00036 #include "AbstractLinAlgPack_VectorMutable.hpp"
00037 #include "AbstractLinAlgPack_VectorSpace.hpp"
00038 #include "AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00039 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00040 #include "DenseLinAlgPack_DMatrixOp.hpp"
00041
00042 void LinAlgOpPack::Mp_StM(
00043 DMatrixSlice* C, value_type a
00044 ,const MatrixOp& B, BLAS_Cpp::Transp B_trans
00045 )
00046 {
00047 using AbstractLinAlgPack::VectorSpace;
00048 using AbstractLinAlgPack::VectorDenseEncap;
00049 using AbstractLinAlgPack::MatrixOpGetGMS;
00050 using AbstractLinAlgPack::MatrixDenseEncap;
00051 const MatrixOpGetGMS
00052 *B_get_gms = dynamic_cast<const MatrixOpGetGMS*>(&B);
00053 if(B_get_gms) {
00054 DenseLinAlgPack::Mp_StM( C, a, MatrixDenseEncap(*B_get_gms)(), B_trans );
00055 }
00056 else {
00057 const size_type num_cols = C->cols();
00058 VectorSpace::multi_vec_mut_ptr_t
00059 B_mv = ( B_trans == BLAS_Cpp::no_trans
00060 ? B.space_cols() : B.space_rows()
00061 ).create_members(num_cols);
00062 assign(B_mv.get(),B,B_trans);
00063 for( size_type j = 1; j <= num_cols; ++j ) {
00064 DenseLinAlgPack::Vp_StV(&C->col(j),a,VectorDenseEncap(*B_mv->col(j))());
00065 }
00066 }
00067 }
00068
00069 void LinAlgOpPack::Vp_StMtV(
00070 DVectorSlice* y, value_type a, const MatrixOp& M
00071 ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x, value_type b
00072 )
00073 {
00074 using BLAS_Cpp::no_trans;
00075 using AbstractLinAlgPack::VectorDenseMutableEncap;
00076 VectorSpace::vec_mut_ptr_t
00077 ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member(),
00078 ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
00079 (VectorDenseMutableEncap(*ay))() = *y;
00080 (VectorDenseMutableEncap(*ax))() = x;
00081 AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, *ax, b );
00082 *y = VectorDenseMutableEncap(*ay)();
00083 }
00084
00085 void LinAlgOpPack::Vp_StMtV(
00086 DVectorSlice* y, value_type a, const MatrixOp& M
00087 ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x, value_type b
00088 )
00089 {
00090 using BLAS_Cpp::no_trans;
00091 using AbstractLinAlgPack::VectorDenseMutableEncap;
00092 VectorSpace::vec_mut_ptr_t
00093 ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).create_member();
00094 (VectorDenseMutableEncap(*ay))() = *y;
00095 AbstractLinAlgPack::Vp_StMtV( ay.get(), a, M, M_trans, x, b );
00096 *y = VectorDenseMutableEncap(*ay)();
00097 }
00098
00099 void LinAlgOpPack::V_InvMtV(
00100 DVectorSlice* y, const MatrixOpNonsing& M
00101 ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00102 )
00103 {
00104 using BLAS_Cpp::trans;
00105 using AbstractLinAlgPack::VectorDenseMutableEncap;
00106 VectorSpace::vec_mut_ptr_t
00107 ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
00108 ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
00109 (VectorDenseMutableEncap(*ay))() = *y;
00110 (VectorDenseMutableEncap(*ax))() = x;
00111 AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
00112 *y = VectorDenseMutableEncap(*ay)();
00113 }
00114
00115 void LinAlgOpPack::V_InvMtV(
00116 DVector* y, const MatrixOpNonsing& M
00117 ,BLAS_Cpp::Transp M_trans, const DVectorSlice& x
00118 )
00119 {
00120 using BLAS_Cpp::trans;
00121 using AbstractLinAlgPack::VectorDenseMutableEncap;
00122 VectorSpace::vec_mut_ptr_t
00123 ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member(),
00124 ax = ( M_trans == trans ? M.space_rows() : M.space_cols() ).create_member();
00125 (VectorDenseMutableEncap(*ax))() = x;
00126 AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, *ax );
00127 *y = VectorDenseMutableEncap(*ay)();
00128 }
00129
00130 void LinAlgOpPack::V_InvMtV(
00131 DVectorSlice* y, const MatrixOpNonsing& M
00132 ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00133 )
00134 {
00135 using BLAS_Cpp::trans;
00136 using AbstractLinAlgPack::VectorDenseMutableEncap;
00137 VectorSpace::vec_mut_ptr_t
00138 ay = ( M_trans == trans ? M.space_cols() : M.space_rows() ).create_member();
00139 AbstractLinAlgPack::V_InvMtV( ay.get(), M, M_trans, x );
00140 *y = VectorDenseMutableEncap(*ay)();
00141 }
00142
00143 void LinAlgOpPack::V_InvMtV(
00144 DVector* y, const MatrixOpNonsing& M
00145 ,BLAS_Cpp::Transp M_trans, const SpVectorSlice& x
00146 )
00147 {
00148 y->resize(M.rows());
00149 LinAlgOpPack::V_InvMtV( &(*y)(), M, M_trans, x );
00150 }
00151
00152
00153
00154
00155
00156
00157
00158 void LinAlgOpPack::Vp_StPtMtV(
00159 DVectorSlice* y, value_type a
00160 ,const GenPermMatrixSlice& P, BLAS_Cpp::Transp P_trans
00161 ,const MatrixOp& M, BLAS_Cpp::Transp M_trans
00162 ,const DVectorSlice& x, value_type b
00163 )
00164 {
00165 namespace mmp = MemMngPack;
00166 using BLAS_Cpp::no_trans;
00167 using AbstractLinAlgPack::VectorMutableDense;
00168 using AbstractLinAlgPack::VectorDenseMutableEncap;
00169 using AbstractLinAlgPack::Vp_StPtMtV;
00170 VectorSpace::space_ptr_t
00171 ay_space = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans);
00172 VectorSpace::vec_mut_ptr_t
00173 ay = ( ay_space.get()
00174 ? ay_space->create_member()
00175 : Teuchos::rcp_implicit_cast<VectorMutable>(
00176 Teuchos::rcp(new VectorMutableDense(BLAS_Cpp::rows(P.rows(),P.cols(),P_trans)))
00177 ) ),
00178 ax = ( M_trans == no_trans ? M.space_rows() : M.space_cols() ).create_member();
00179 (VectorDenseMutableEncap(*ay))() = *y;
00180 (VectorDenseMutableEncap(*ax))() = x;
00181 Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, *ax, b );
00182 *y = VectorDenseMutableEncap(*ay)();
00183 }
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 SpVectorSlice& x, value_type b
00190 )
00191 {
00192 using BLAS_Cpp::no_trans;
00193 using AbstractLinAlgPack::VectorMutableDense;
00194 using AbstractLinAlgPack::VectorDenseMutableEncap;
00195 using AbstractLinAlgPack::Vp_StPtMtV;
00196 VectorSpace::vec_mut_ptr_t
00197 ay = ( M_trans == no_trans ? M.space_cols() : M.space_rows() ).space(P,P_trans)->create_member();
00198 (VectorDenseMutableEncap(*ay))() = *y;
00199 Vp_StPtMtV( ay.get(), a, P, P_trans, M, M_trans, x, b );
00200 *y = VectorDenseMutableEncap(*ay)();
00201 }