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