AbstractLinAlgPack_MatrixConvertToSparseEncap.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 <assert.h>
00030 
00031 #include "AbstractLinAlgPack_MatrixConvertToSparseEncap.hpp"
00032 #include "AbstractLinAlgPack_MatrixExtractSparseElements.hpp"
00033 #include "AbstractLinAlgPack_VectorSpace.hpp"
00034 #include "DenseLinAlgPack_IVector.hpp"
00035 #include "Teuchos_TestForException.hpp"
00036 
00037 namespace AbstractLinAlgPack {
00038 
00039 // Constructors/initializers
00040 
00041 MatrixConvertToSparseEncap::MatrixConvertToSparseEncap()
00042   :row_rng_(Range1D::Invalid)
00043   ,col_rng_(Range1D::Invalid)
00044   ,mese_trans_(BLAS_Cpp::no_trans)
00045   ,alpha_(0.0)
00046   ,nz_full_(0)
00047 {}
00048 
00049 MatrixConvertToSparseEncap::MatrixConvertToSparseEncap(
00050   const mese_ptr_t           &mese
00051   ,const i_vector_ptr_t      &inv_row_perm
00052   ,const Range1D             &row_rng
00053   ,const i_vector_ptr_t      &inv_col_perm
00054   ,const Range1D             &col_rng
00055   ,const BLAS_Cpp::Transp    mese_trans
00056   ,const value_type          alpha
00057   )
00058 {
00059   this->initialize(mese,inv_row_perm,row_rng,inv_col_perm,col_rng,mese_trans,alpha);
00060 }
00061 
00062 void MatrixConvertToSparseEncap::initialize(
00063   const mese_ptr_t           &mese
00064   ,const i_vector_ptr_t      &inv_row_perm
00065   ,const Range1D             &row_rng_in
00066   ,const i_vector_ptr_t      &inv_col_perm
00067   ,const Range1D             &col_rng_in
00068   ,const BLAS_Cpp::Transp    mese_trans
00069   ,const value_type          alpha
00070   )
00071 {
00072   const size_type mese_rows = mese->rows(), mese_cols = mese->cols();
00073   const Range1D row_rng = RangePack::full_range(row_rng_in,1,mese_rows);
00074   const Range1D col_rng = RangePack::full_range(col_rng_in,1,mese_cols);
00075 #ifdef TEUCHOS_DEBUG
00076   const char msg_head[] = "MatrixConvertToSparseEncap::initialize(...): Error!";
00077   TEST_FOR_EXCEPTION( mese.get() == NULL, std::logic_error, msg_head );
00078   TEST_FOR_EXCEPTION( inv_row_perm.get() != NULL && inv_row_perm->size() != mese_rows, std::logic_error, msg_head );
00079   TEST_FOR_EXCEPTION( row_rng.ubound() > mese_rows, std::logic_error, msg_head );
00080   TEST_FOR_EXCEPTION( inv_col_perm.get() != NULL && inv_col_perm->size() != mese_cols, std::logic_error, msg_head );
00081   TEST_FOR_EXCEPTION( col_rng.ubound() > mese->cols(), std::logic_error, msg_head );
00082 #endif
00083   mese_           = mese;
00084   mese_trans_     = mese_trans;
00085   alpha_          = alpha;
00086   inv_row_perm_   = inv_row_perm;
00087   inv_col_perm_   = inv_col_perm;
00088   row_rng_        = row_rng;
00089   col_rng_        = col_rng;
00090   nz_full_        = this->num_nonzeros(EXTRACT_FULL_MATRIX,ELEMENTS_ALLOW_DUPLICATES_SUM);
00091   space_cols_     = ( mese_trans_ == BLAS_Cpp::no_trans
00092             ? mese_->space_cols().sub_space(row_rng_)
00093             : mese_->space_rows().sub_space(col_rng_) );
00094   space_rows_     = ( mese_trans_ == BLAS_Cpp::no_trans
00095             ? mese_->space_rows().sub_space(col_rng_)
00096             : mese_->space_cols().sub_space(row_rng_) );
00097 }
00098 
00099 void MatrixConvertToSparseEncap::set_uninitialized()
00100 {
00101   namespace mmp   = MemMngPack;
00102   mese_           = Teuchos::null;
00103   inv_row_perm_   = Teuchos::null;
00104   row_rng_        = Range1D::Invalid;
00105   inv_col_perm_   = Teuchos::null;
00106   col_rng_        = Range1D::Invalid;
00107   mese_trans_     = BLAS_Cpp::no_trans;
00108   alpha_          = 0.0;
00109   nz_full_        = 0;
00110 }
00111 
00112 // Overridden from MatrixBase
00113 
00114 const VectorSpace& MatrixConvertToSparseEncap::space_cols() const
00115 {
00116   return *space_cols_;
00117 }
00118 
00119 const VectorSpace& MatrixConvertToSparseEncap::space_rows() const
00120 {
00121   return *space_rows_;
00122 }
00123 
00124 size_type MatrixConvertToSparseEncap::rows() const
00125 {
00126   return mese_trans_ == BLAS_Cpp::no_trans ? row_rng_.size() : col_rng_.size();
00127 }
00128 
00129 size_type MatrixConvertToSparseEncap::cols() const
00130 {
00131   return mese_trans_ == BLAS_Cpp::no_trans ? col_rng_.size() : row_rng_.size();
00132 }
00133 
00134 size_type MatrixConvertToSparseEncap::nz() const
00135 {
00136   return nz_full_;
00137 }
00138 
00139 // Overridden from MatrixConvertToSparse
00140 
00141 index_type MatrixConvertToSparseEncap::num_nonzeros(
00142   EExtractRegion        extract_region_in
00143   ,EElementUniqueness   element_uniqueness
00144   ) const
00145 {
00146   index_type dl = 0, du = 0;
00147   EExtractRegion extract_region = extract_region_in;
00148   if( mese_trans_ == BLAS_Cpp::trans )
00149     extract_region
00150       = ( extract_region_in == EXTRACT_FULL_MATRIX
00151         ? EXTRACT_FULL_MATRIX
00152         : ( extract_region_in == EXTRACT_UPPER_TRIANGULAR
00153           ? EXTRACT_LOWER_TRIANGULAR
00154           : EXTRACT_UPPER_TRIANGULAR ) );
00155   switch(extract_region) {
00156     case EXTRACT_FULL_MATRIX:
00157       dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
00158       du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
00159       break;
00160     case EXTRACT_UPPER_TRIANGULAR:
00161       dl = -(index_type)row_rng_.lbound() + (index_type)col_rng_.lbound();
00162       du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
00163       break;
00164     case EXTRACT_LOWER_TRIANGULAR:
00165       dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
00166       du = +(index_type)col_rng_.lbound() - (index_type)row_rng_.lbound();
00167       break;
00168     default:
00169       TEST_FOR_EXCEPT(true);
00170   }
00171   const index_type
00172     *inv_row_perm = inv_row_perm_.get() ? &(*inv_row_perm_)(1) : NULL,
00173     *inv_col_perm = inv_col_perm_.get() ? &(*inv_col_perm_)(1) : NULL;
00174   return mese_->count_nonzeros(
00175     element_uniqueness
00176     ,inv_row_perm
00177     ,inv_col_perm
00178     ,row_rng_
00179     ,col_rng_
00180     ,dl
00181     ,du
00182     );
00183 }
00184 
00185 void MatrixConvertToSparseEncap::coor_extract_nonzeros(
00186   EExtractRegion                extract_region
00187   ,EElementUniqueness           element_uniqueness
00188   ,const index_type             len_Aval
00189   ,value_type                   Aval[]
00190   ,const index_type             len_Aij
00191   ,index_type                   Arow[]
00192   ,index_type                   Acol[]
00193   ,const index_type             row_offset
00194   ,const index_type             col_offset
00195   ) const
00196 {
00197   index_type dl = 0, du = 0;  // This may not be correct!
00198   switch(extract_region) {
00199     case EXTRACT_FULL_MATRIX:
00200       dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
00201       du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
00202       break;
00203     case EXTRACT_UPPER_TRIANGULAR:
00204       dl = -(index_type)row_rng_.lbound() + (index_type)col_rng_.lbound();
00205       du = +(index_type)col_rng_.ubound() - (index_type)row_rng_.lbound();
00206       break;
00207     case EXTRACT_LOWER_TRIANGULAR:
00208       dl = -(index_type)row_rng_.ubound() + (index_type)col_rng_.lbound();
00209       du = +(index_type)col_rng_.lbound() - (index_type)row_rng_.lbound();
00210       break;
00211     default:
00212       TEST_FOR_EXCEPT(true);
00213   }
00214   const index_type
00215     *inv_row_perm = inv_row_perm_.get() ? &(*inv_row_perm_)(1) : NULL,
00216     *inv_col_perm = inv_col_perm_.get() ? &(*inv_col_perm_)(1) : NULL;
00217   mese_->coor_extract_nonzeros(
00218     element_uniqueness
00219     ,inv_row_perm
00220     ,inv_col_perm
00221     ,row_rng_
00222     ,col_rng_
00223     ,dl
00224     ,du
00225     ,alpha_
00226     ,len_Aval
00227     ,Aval
00228     ,len_Aij
00229     ,mese_trans_ == BLAS_Cpp::no_trans ? Arow : Acol
00230     ,mese_trans_ == BLAS_Cpp::no_trans ? Acol : Arow
00231     ,mese_trans_ == BLAS_Cpp::no_trans ? row_offset : col_offset
00232     ,mese_trans_ == BLAS_Cpp::no_trans ? col_offset : row_offset
00233     );
00234 }
00235 
00236 } // end namespace AbstractLinAlgPack 

Generated on Tue Jul 13 09:28:51 2010 for AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects by  doxygen 1.4.7