AbstractLinAlgPack_MatrixPermAggr.cpp

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 "AbstractLinAlgPack_MatrixPermAggr.hpp"
00030 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00031 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00032 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00033 #include "AbstractLinAlgPack_VectorSpace.hpp"
00034 #include "AbstractLinAlgPack_Permutation.hpp"
00035 #include "AbstractLinAlgPack_PermutationOut.hpp"
00036 #include "Teuchos_TestForException.hpp"
00037 #include "Teuchos_dyn_cast.hpp"
00038 
00039 namespace AbstractLinAlgPack {
00040 
00041 // Constructors / initializers
00042 
00043 MatrixPermAggr::MatrixPermAggr()
00044 {} // Nothing to explicitly initialize
00045 
00046 MatrixPermAggr::MatrixPermAggr(
00047   const mat_ptr_t      &mat_orig
00048   ,const perm_ptr_t    &row_perm
00049   ,const perm_ptr_t    &col_perm
00050   ,const mat_ptr_t     &mat_perm
00051   )
00052 {
00053   this->initialize(mat_orig,row_perm,col_perm,mat_perm);
00054 }
00055 
00056 void MatrixPermAggr::initialize(
00057   const mat_ptr_t      &mat_orig
00058   ,const perm_ptr_t    &row_perm
00059   ,const perm_ptr_t    &col_perm
00060   ,const mat_ptr_t     &mat_perm
00061   )
00062 {
00063 #ifdef TEUCHOS_DEBUG
00064   TEST_FOR_EXCEPTION(
00065     mat_orig.get() == NULL, std::invalid_argument
00066     ,"MatrixPermAggr::initialize(...): Error!" );
00067 #endif
00068 #ifdef ABSTRACTLINALGPACK_ASSERT_COMPATIBILITY
00069   bool is_compatible = false;
00070   if(row_perm.get()) {
00071     is_compatible = mat_orig->space_cols().is_compatible(row_perm->space());
00072     TEST_FOR_EXCEPTION(
00073       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00074       ,"MatrixPermAggr::initialize(...): Error, " 
00075       "mat_orig->space_cols().is_compatible(row_perm->space()) == false" );
00076   }
00077   if(col_perm.get()) {
00078     is_compatible = mat_orig->space_rows().is_compatible(col_perm->space());
00079     TEST_FOR_EXCEPTION(
00080       !is_compatible, VectorSpace::IncompatibleVectorSpaces
00081       ,"MatrixPermAggr::initialize(...): Error, " 
00082       "mat_orig->space_rows().is_compatible(col_perm->space()) == false" );
00083   }
00084 #endif
00085   mat_orig_   = mat_orig;
00086   row_perm_   = row_perm;
00087   col_perm_   = col_perm;
00088   mat_perm_   = mat_perm;
00089 }
00090 
00091 void MatrixPermAggr::set_uninitialized()
00092 {
00093   namespace rcp = MemMngPack;
00094   mat_orig_   = Teuchos::null;
00095   row_perm_   = Teuchos::null;
00096   col_perm_   = Teuchos::null;
00097   mat_perm_   = Teuchos::null;
00098 }
00099 
00100 // Overridden from MatrixBase
00101 
00102 size_type MatrixPermAggr::rows() const
00103 {
00104   return mat_orig_.get() ? mat_orig_->rows() : 0;
00105 }
00106 
00107 size_type MatrixPermAggr::cols() const
00108 {
00109   return mat_orig_.get() ? mat_orig_->cols() : 0;
00110 }
00111 
00112 size_type MatrixPermAggr::nz() const
00113 {
00114   return mat_orig_.get() ? mat_orig_->nz() : 0;
00115 }
00116 
00117 // Overridden from MatrixOp
00118 
00119 const VectorSpace& MatrixPermAggr::space_cols() const
00120 {
00121   return mat_orig_->space_cols();
00122 }
00123 
00124 const VectorSpace& MatrixPermAggr::space_rows() const
00125 {
00126   return mat_orig_->space_rows();
00127 }
00128 
00129 MatrixOp::mat_ptr_t
00130 MatrixPermAggr::sub_view(const Range1D& row_rng, const Range1D& col_rng) const
00131 {
00132   if(mat_perm_.get())
00133     return mat_perm_->sub_view(row_rng,col_rng);
00134   if(!row_perm_.get() && !col_perm_.get())
00135     return mat_orig_->sub_view(row_rng,col_rng);
00136   return MatrixOp::sub_view(row_rng,col_rng); // ToDo: Speicalized!
00137 }
00138 
00139 MatrixOp& MatrixPermAggr::operator=(const MatrixOp& M)
00140 {
00141   using Teuchos::dyn_cast;
00142   const MatrixPermAggr
00143     Mp = dyn_cast<const MatrixPermAggr>(M);
00144   if( this == &Mp )
00145     return *this; // Assignment to self
00146   // Shallow copy is okay as long as client is careful!
00147   mat_orig_  = Mp.mat_orig_;
00148   row_perm_  = Mp.row_perm_;
00149   col_perm_  = Mp.col_perm_;
00150   mat_perm_  = Mp.mat_perm_;
00151   return *this;
00152 }
00153 
00154 std::ostream& MatrixPermAggr::output(std::ostream& out) const
00155 {
00156   out << "Matrix with permuted view:\n";
00157   out << "mat_orig =\n" << *mat_orig_;
00158   out << "row_perm =";
00159   if( row_perm_.get() )
00160     out << "\n" << *row_perm_;
00161   else
00162     out << " NULL\n";
00163   out << "col_perm =";
00164   if( col_perm_.get() )
00165     out << "\n" << *col_perm_;
00166   else
00167     out << " NULL\n";
00168   out << "mat_perm =";
00169   if( mat_perm_.get() )
00170     out << "\n" << *mat_perm_;
00171   else
00172     out << " NULL\n";
00173   return out;
00174 }
00175 
00176 bool MatrixPermAggr::Mp_StM(
00177   MatrixOp* mwo_lhs, value_type alpha
00178   ,BLAS_Cpp::Transp trans_rhs
00179   ) const
00180 {
00181   if(!row_perm_.get() && !col_perm_.get()) {
00182     AbstractLinAlgPack::Mp_StM(mwo_lhs,alpha,*mat_orig_,trans_rhs);
00183     return true;
00184   }
00185   AbstractLinAlgPack::Mp_StM(mwo_lhs,alpha,*mat_orig_,trans_rhs); // ToDo: Specialize!
00186   return true;
00187 }
00188 
00189 bool MatrixPermAggr::Mp_StMtP(
00190   MatrixOp* mwo_lhs, value_type alpha
00191   ,BLAS_Cpp::Transp M_trans
00192   ,const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00193   ) const
00194 {
00195   if(!row_perm_.get() && !col_perm_.get()) {
00196     AbstractLinAlgPack::Mp_StMtP(mwo_lhs,alpha,*mat_orig_,M_trans,P_rhs,P_rhs_trans);
00197     return true;
00198   }
00199   AbstractLinAlgPack::Mp_StMtP(mwo_lhs,alpha,*mat_orig_,M_trans,P_rhs,P_rhs_trans); // ToDo: Specialize!
00200   return true;
00201 }
00202 
00203 bool MatrixPermAggr::Mp_StPtM(
00204   MatrixOp* mwo_lhs, value_type alpha
00205   , const GenPermMatrixSlice& P_rhs, BLAS_Cpp::Transp P_rhs_trans
00206   , BLAS_Cpp::Transp M_trans
00207   ) const
00208 {
00209   if(!row_perm_.get() && !col_perm_.get()) {
00210     AbstractLinAlgPack::Mp_StPtM(mwo_lhs,alpha,P_rhs,P_rhs_trans,*mat_orig_,M_trans);
00211     return true;
00212   }
00213   AbstractLinAlgPack::Mp_StPtM(mwo_lhs,alpha,P_rhs,P_rhs_trans,*mat_orig_,M_trans); // ToDo: Specialize!
00214   return true;
00215 }
00216 
00217 bool MatrixPermAggr::Mp_StPtMtP(
00218   MatrixOp* mwo_lhs, value_type alpha
00219   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00220   ,BLAS_Cpp::Transp M_trans
00221   ,const GenPermMatrixSlice& P_rhs2, BLAS_Cpp::Transp P_rhs2_trans
00222   ) const
00223 {
00224   if(!row_perm_.get() && !col_perm_.get()) {
00225     AbstractLinAlgPack::Mp_StPtMtP(mwo_lhs,alpha,P_rhs1,P_rhs1_trans,*mat_orig_,M_trans,P_rhs2,P_rhs2_trans);
00226     return true;
00227   }
00228   AbstractLinAlgPack::Mp_StPtMtP(mwo_lhs,alpha,P_rhs1,P_rhs1_trans,*mat_orig_,M_trans,P_rhs2,P_rhs2_trans); // ToDo: Specialize!
00229   return true;
00230 }
00231 
00232 void MatrixPermAggr::Vp_StMtV(
00233   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00234   ,const Vector& x, value_type b
00235   ) const
00236 {
00237   using BLAS_Cpp::no_trans;
00238   using BLAS_Cpp::trans;
00239   using LinAlgOpPack::V_MtV;
00240 
00241   if(mat_perm_.get()) {
00242     AbstractLinAlgPack::Vp_StMtV(y,a,*mat_perm_,M_trans,x,b);
00243     return;
00244   }
00245   if(!row_perm_.get() && !col_perm_.get()) {
00246     AbstractLinAlgPack::Vp_StMtV(y,a,*mat_orig_,M_trans,x,b);
00247     return;
00248   }
00249 
00250   //
00251   // y = a*op(Pr*M*Pc')*x + b*y
00252   //
00253   // =>
00254   //
00255   // y = a*P1*op(M)*P2'*x + b*y
00256   //
00257   // =>
00258   //
00259   // ta = P2'*x
00260   // tb = op(M)*ta
00261   // tc = P1*tb
00262   // y = a*tc + b*y
00263   //
00264   const Permutation
00265     *P1 = ( M_trans == no_trans ? row_perm_.get() : col_perm_.get() ),
00266     *P2 = ( M_trans == no_trans ? col_perm_.get() : row_perm_.get() );
00267   VectorSpace::vec_mut_ptr_t  ta, tb, tc;
00268   // ta = P2'*x
00269   if( P2 && !P2->is_identity() )
00270     P2->permute( trans, x, (ta = P2->space().create_member()).get() );
00271   else
00272     *(tb = ( M_trans == no_trans ? mat_orig_->space_rows() : mat_orig_->space_cols() ).create_member() ) = x;
00273   // tb = op(M)*ta
00274   V_MtV(
00275     (tb = ( M_trans == no_trans ? mat_orig_->space_cols() : mat_orig_->space_rows() ).create_member() ).get()
00276     ,*mat_orig_, M_trans, *ta
00277     );
00278   // tc = P1*tb
00279   if( P1 && !P1->is_identity() )
00280     P1->permute( no_trans, *tb, (tc = P1->space().create_member()).get() );
00281   else
00282     tc = tb->clone();
00283   // y = b*y + a*tc
00284   AbstractLinAlgPack::Vt_S( y, b );
00285   AbstractLinAlgPack::Vp_StV( y, a, *tc );
00286 }
00287 
00288 void MatrixPermAggr::Vp_StMtV(
00289   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00290   , const SpVectorSlice& x, value_type b) const
00291 {
00292   if(mat_perm_.get()) {
00293     AbstractLinAlgPack::Vp_StMtV(y,a,*mat_perm_,M_trans,x,b);
00294     return;
00295   }
00296   if(!row_perm_.get() && !col_perm_.get()) {
00297     AbstractLinAlgPack::Vp_StMtV(y,a,*mat_orig_,M_trans,x,b);
00298     return;
00299   }
00300   MatrixOp::Vp_StMtV(y,a,M_trans,x,b);
00301 }
00302 
00303 void MatrixPermAggr::Vp_StPtMtV(
00304   VectorMutable* vs_lhs, value_type alpha
00305   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00306   ,BLAS_Cpp::Transp M_rhs2_trans
00307   ,const Vector& v_rhs3, value_type beta
00308   ) const
00309 {
00310   if(!row_perm_.get() && !col_perm_.get()) {
00311     AbstractLinAlgPack::Vp_StPtMtV(vs_lhs,alpha,P_rhs1,P_rhs1_trans,*mat_orig_,M_rhs2_trans,v_rhs3,beta);
00312     return;
00313   }
00314   MatrixOp::Vp_StPtMtV(vs_lhs,alpha,P_rhs1,P_rhs1_trans,M_rhs2_trans,v_rhs3,beta);
00315 }
00316 
00317 void MatrixPermAggr::Vp_StPtMtV(
00318   VectorMutable* vs_lhs, value_type alpha
00319   ,const GenPermMatrixSlice& P_rhs1, BLAS_Cpp::Transp P_rhs1_trans
00320   ,BLAS_Cpp::Transp M_rhs2_trans
00321   ,const SpVectorSlice& sv_rhs3, value_type beta
00322   ) const
00323 {
00324   if(!row_perm_.get() && !col_perm_.get()) {
00325     AbstractLinAlgPack::Vp_StPtMtV(vs_lhs,alpha,P_rhs1,P_rhs1_trans,*mat_orig_,M_rhs2_trans,sv_rhs3,beta);
00326     return;
00327   }
00328   MatrixOp::Vp_StPtMtV(vs_lhs,alpha,P_rhs1,P_rhs1_trans,M_rhs2_trans,sv_rhs3,beta);
00329 }
00330 
00331 value_type MatrixPermAggr::transVtMtV(
00332   const Vector& v_rhs1, BLAS_Cpp::Transp trans_rhs2
00333   ,const Vector& v_rhs3
00334   ) const
00335 {
00336   if(!row_perm_.get() && !col_perm_.get())
00337     return AbstractLinAlgPack::transVtMtV(v_rhs1,*mat_orig_,trans_rhs2,v_rhs3);
00338   return MatrixOp::transVtMtV(v_rhs1,trans_rhs2,v_rhs3);
00339 }
00340 
00341 value_type MatrixPermAggr::transVtMtV(
00342   const SpVectorSlice& sv_rhs1, BLAS_Cpp::Transp trans_rhs2
00343   ,const SpVectorSlice& sv_rhs3
00344   ) const
00345 {
00346   if(!row_perm_.get() && !col_perm_.get())
00347     return AbstractLinAlgPack::transVtMtV(sv_rhs1,*mat_orig_,trans_rhs2,sv_rhs3);
00348   return MatrixOp::transVtMtV(sv_rhs1,trans_rhs2,sv_rhs3);
00349 }
00350 
00351 void MatrixPermAggr::syr2k(
00352   BLAS_Cpp::Transp M_trans, value_type alpha
00353   ,const GenPermMatrixSlice& P1, BLAS_Cpp::Transp P1_trans
00354   ,const GenPermMatrixSlice& P2, BLAS_Cpp::Transp P2_trans
00355   ,value_type beta, MatrixSymOp* symwo_lhs
00356   ) const
00357 {
00358   if(!row_perm_.get() && !col_perm_.get()) {
00359     AbstractLinAlgPack::syr2k(*mat_orig_,M_trans,alpha,P1,P1_trans,P2,P2_trans,beta,symwo_lhs);
00360     return;
00361   }
00362   MatrixOp::syr2k(M_trans,alpha,P1,P1_trans,P2,P2_trans,beta,symwo_lhs);
00363 }
00364 
00365 bool MatrixPermAggr::Mp_StMtM(
00366   MatrixOp* mwo_lhs, value_type alpha
00367   ,BLAS_Cpp::Transp trans_rhs1, const MatrixOp& mwo_rhs2
00368   ,BLAS_Cpp::Transp trans_rhs2, value_type beta
00369   ) const
00370 {
00371   if(!row_perm_.get() && !col_perm_.get()) {
00372     AbstractLinAlgPack::Mp_StMtM(mwo_lhs,alpha,*mat_orig_,trans_rhs1,mwo_rhs2,trans_rhs2,beta);
00373     return true;
00374   }
00375   return MatrixOp::Mp_StMtM(mwo_lhs,alpha,trans_rhs1,mwo_rhs2,trans_rhs2,beta);
00376 }
00377 
00378 bool MatrixPermAggr::Mp_StMtM(
00379   MatrixOp* mwo_lhs, value_type alpha
00380   ,const MatrixOp& mwo_rhs1, BLAS_Cpp::Transp trans_rhs1
00381   ,BLAS_Cpp::Transp trans_rhs2, value_type beta
00382   ) const
00383 {
00384   if(!row_perm_.get() && !col_perm_.get()) {
00385     AbstractLinAlgPack::Mp_StMtM(mwo_lhs,alpha,mwo_rhs1,trans_rhs1,*mat_orig_,trans_rhs2,beta);
00386     return true;
00387   }
00388   return MatrixOp::Mp_StMtM(mwo_lhs,alpha,mwo_rhs1,trans_rhs1,trans_rhs2,beta);
00389 }
00390 
00391 bool MatrixPermAggr::syrk(
00392   BLAS_Cpp::Transp M_trans, value_type alpha
00393   ,value_type beta, MatrixSymOp* sym_lhs
00394   ) const
00395 {
00396   if(!row_perm_.get() && !col_perm_.get()) {
00397     AbstractLinAlgPack::syrk(*mat_orig_,M_trans,alpha,beta,sym_lhs);
00398     return true;
00399   }
00400   return MatrixOp::syrk(M_trans,alpha,beta,sym_lhs);
00401 }
00402 
00403 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:09:15 2011 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.6.3