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