AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_MatrixSparseCOORSerial.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_MatrixSparseCOORSerial.hpp"
00045 #include "AbstractLinAlgPack_MatrixCOORTmplItfc.hpp"
00046 #include "AbstractLinAlgPack_COOMatrixTmplOp.hpp"
00047 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00048 #include "AbstractLinAlgPack_AssertOp.hpp"
00049 #include "AbstractLinAlgPack_LinAlgOpPackHack.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 #include "Teuchos_dyn_cast.hpp"
00052 
00053 namespace AbstractLinAlgPack {
00054 
00055 MatrixSparseCOORSerial::ReleaseValRowColArrays::~ReleaseValRowColArrays()
00056 {
00057   if(owns_mem_) {
00058     if(val_)   delete [] val_;
00059     if(row_i_) delete [] row_i_;
00060     if(col_j_) delete [] col_j_;
00061   }
00062 }
00063 
00064 bool MatrixSparseCOORSerial::ReleaseValRowColArrays::resource_is_bound() const
00065 {
00066   return val_ != NULL;
00067 }
00068 
00069 // static members
00070 
00071 MatrixSparseCOORSerial::release_resource_ptr_t
00072 MatrixSparseCOORSerial::release_resource_null_ = Teuchos::null;
00073 
00074 // Constructors / initializers
00075 
00076 MatrixSparseCOORSerial::MatrixSparseCOORSerial()
00077   :rows_(0)
00078   ,cols_(0)
00079   ,max_nz_(0)
00080   ,nz_(0)
00081   ,val_(NULL)
00082   ,row_i_(NULL)
00083   ,col_j_(NULL)
00084   ,self_allocate_(true)
00085 {}
00086 
00087 void MatrixSparseCOORSerial::set_buffers(
00088   size_type                      max_nz
00089   ,value_type                    *val
00090   ,index_type                    *row_i
00091   ,index_type                    *col_j
00092   ,const release_resource_ptr_t  &release_resource
00093   ,size_type                     rows
00094   ,size_type                     cols
00095   ,size_type                     nz
00096   ,bool                          check_input
00097   )
00098 {
00099 #ifdef TEUCHOS_DEBUG
00100   const char msg_err[] = "MatrixSparseCOORSerial::set_buffer(...) : Error,!";
00101   TEUCHOS_TEST_FOR_EXCEPTION( max_nz <= 0, std::invalid_argument, msg_err );
00102   TEUCHOS_TEST_FOR_EXCEPTION( val == NULL || row_i == NULL || col_j == NULL, std::invalid_argument, msg_err );
00103   TEUCHOS_TEST_FOR_EXCEPTION( rows > 0 && cols <= 0 , std::invalid_argument, msg_err );
00104   TEUCHOS_TEST_FOR_EXCEPTION( rows > 0 && (nz < 0 || nz > max_nz), std::invalid_argument, msg_err );
00105 #endif
00106   max_nz_           = max_nz;
00107   val_              = val;
00108   row_i_            = row_i;
00109   col_j_            = col_j;
00110   release_resource_ = release_resource;
00111   self_allocate_    = false;
00112   if(rows) {
00113     rows_ = rows;
00114     cols_ = cols;
00115     nz_   = nz;
00116     space_cols_.initialize(rows);
00117     space_rows_.initialize(cols);
00118     if( nz && check_input ) {
00119       TEUCHOS_TEST_FOR_EXCEPT(true); // Todo: Check that row_i[] and col_j[] are in bounds
00120     }
00121   }
00122   else {
00123     rows_ = 0;
00124     cols_ = 0;
00125     nz_   = 0;
00126     space_cols_.initialize(0);
00127     space_rows_.initialize(0);
00128   }
00129 }
00130 
00131 void MatrixSparseCOORSerial::set_uninitialized()
00132 {
00133   max_nz_           = 0;
00134   val_              = NULL;
00135   row_i_            = NULL;
00136   col_j_            = NULL;
00137   release_resource_ = Teuchos::null;
00138   self_allocate_    = true;
00139   rows_             = 0;
00140   cols_             = 0;
00141   nz_               = 0;
00142   space_cols_.initialize(0);
00143   space_rows_.initialize(0);
00144 }
00145 
00146 // Overridden from MatrixBase
00147 
00148 size_type MatrixSparseCOORSerial::rows() const
00149 {
00150   return rows_;
00151 }
00152 
00153 size_type MatrixSparseCOORSerial::cols() const
00154 {
00155   return cols_;
00156 }
00157 
00158 size_type MatrixSparseCOORSerial::nz() const
00159 {
00160   return nz_;
00161 }
00162 
00163 // Overridden from MatrixOp
00164 
00165 const VectorSpace& MatrixSparseCOORSerial::space_cols() const
00166 {
00167   return space_cols_;
00168 }
00169 
00170 const VectorSpace& MatrixSparseCOORSerial::space_rows() const
00171 {
00172   return space_rows_;
00173 }
00174 
00175 MatrixOp& MatrixSparseCOORSerial::operator=(const MatrixOp& M)
00176 {
00177   using Teuchos::dyn_cast;
00178   const MatrixSparseCOORSerial
00179     &Mc = dyn_cast<const MatrixSparseCOORSerial>(M);
00180   if( this == &Mc )
00181     return *this; // assignment to self
00182   // A shallow copy is fine as long as we are carefull.
00183   max_nz_           = Mc.max_nz_;
00184   val_              = Mc.val_;
00185   row_i_            = Mc.row_i_;
00186   col_j_            = Mc.col_j_;
00187   release_resource_ = Mc.release_resource_;
00188   self_allocate_    = Mc.self_allocate_;
00189   rows_             = Mc.rows_;
00190   cols_             = Mc.cols_;
00191   nz_               = Mc.nz_;
00192   space_cols_.initialize(rows_);
00193   space_rows_.initialize(cols_);
00194   return *this;
00195 }
00196 
00197 std::ostream& MatrixSparseCOORSerial::output(std::ostream& out) const
00198 {
00199   const size_type
00200     rows = this->rows(),
00201     cols = this->cols(),
00202     nz   = this->nz();
00203   out
00204     << "Sparse " << rows << " x " << cols << " matrix with "
00205     << nz << " nonzero entries:\n";
00206   const value_type
00207     *itr_val      = &val_[0],
00208     *itr_val_end  = itr_val + nz_;
00209   const index_type
00210     *itr_ivect    = &row_i_[0],
00211     *itr_jvect    = &col_j_[0];
00212   for(; itr_val != itr_val_end; ++itr_val, ++itr_ivect, ++itr_jvect)
00213     out << *itr_val << ":" << *itr_ivect << ":" << *itr_jvect << " ";
00214   out << "\n";
00215   
00216   if( rows * cols <= 400 ) {
00217     out << "Converted to dense =\n";
00218     MatrixOp::output(out);
00219   }
00220 
00221   return out;
00222 }
00223 
00224 void MatrixSparseCOORSerial::Vp_StMtV(
00225   VectorMutable* y, value_type a, BLAS_Cpp::Transp M_trans
00226   , const Vector& x, value_type b
00227   ) const
00228 {
00229   AbstractLinAlgPack::Vp_MtV_assert_compatibility( y, *this, M_trans, x );
00230   VectorDenseMutableEncap   y_d(*y);
00231   VectorDenseEncap          x_d(x);
00232   DenseLinAlgPack::Vt_S(&y_d(),b);
00233   Vp_StCOOMtV(
00234     &y_d(),a,MatrixCOORTmplItfc<value_type,index_type>(
00235       rows_,cols_,nz_,0,0,val_,row_i_,col_j_ )
00236     ,M_trans, x_d()
00237     );
00238 }
00239 
00240 // Overridden from MatrixLoadSparseElements
00241 
00242 void MatrixSparseCOORSerial::reinitialize(
00243   size_type                  rows
00244   ,size_type                 cols
00245   ,size_type                 max_nz
00246   ,EAssumeElementUniqueness  element_uniqueness
00247   )
00248 {
00249   namespace rcp = MemMngPack;
00250 #ifdef TEUCHOS_DEBUG
00251   const char msg_err_head[] = "MatrixSparseCOORSerial::reinitialize(...) : Error";
00252   TEUCHOS_TEST_FOR_EXCEPTION( max_nz <= 0, std::invalid_argument, msg_err_head<<"!" );
00253   TEUCHOS_TEST_FOR_EXCEPTION( rows <= 0 || cols <= 0 , std::invalid_argument, msg_err_head<<"!" );
00254 #endif
00255   rows_               = rows;
00256   cols_               = cols; 
00257   element_uniqueness_ = element_uniqueness;
00258   if( self_allocate_ ) {
00259     if(max_nz_ < max_nz) {
00260       release_resource_ = Teuchos::rcp(
00261         new ReleaseValRowColArrays(
00262           val_    = new value_type[max_nz]
00263           ,row_i_ = new index_type[max_nz]
00264           ,col_j_ = new index_type[max_nz]
00265           )
00266         );
00267       max_nz_ = max_nz;
00268     }
00269 
00270   }
00271   else {
00272 #ifdef TEUCHOS_DEBUG
00273     TEUCHOS_TEST_FOR_EXCEPTION(
00274       max_nz <= max_nz_, std::invalid_argument
00275       ,msg_err_head << "Buffers set up by client in set_buffers() only allows storage for "
00276       "max_nz_ = " << max_nz_ << " nonzero entries while client requests storage for "
00277       "max_nz = " << max_nz << " nonzero entries!" );
00278 #endif
00279   }
00280   reload_val_only_         = false;
00281   reload_val_only_nz_last_ = 0;
00282   nz_                      = 0;
00283   max_nz_load_             = 0;
00284 }
00285 
00286 void MatrixSparseCOORSerial::reset_to_load_values()
00287 {
00288 #ifdef TEUCHOS_DEBUG
00289   TEUCHOS_TEST_FOR_EXCEPTION(
00290     rows_ == 0 || cols_ == 0, std::invalid_argument
00291     ,"MatrixSparseCOORSerial::reset_to_load_values(...) : Error, "
00292     "this matrix is not initialized so it can't be rest to load "
00293     "new values for nonzero entries." );
00294 #endif
00295   reload_val_only_         = true;
00296   reload_val_only_nz_last_ = nz_;
00297   nz_                      = 0;
00298   max_nz_load_             = 0;
00299 }
00300 
00301 void MatrixSparseCOORSerial::get_load_nonzeros_buffers(
00302   size_type      max_nz_load
00303   ,value_type    **val
00304   ,index_type    **row_i
00305   ,index_type    **col_j
00306   )
00307 {
00308 #ifdef TEUCHOS_DEBUG
00309   TEUCHOS_TEST_FOR_EXCEPTION(
00310     max_nz_load_ != 0 , std::logic_error
00311     ,"MatrixSparseCOORSerial::get_load_nonzeros_buffers(...) : Error, "
00312     "You must call commit_load_nonzeros_buffers() between calls to this method!" );
00313   TEUCHOS_TEST_FOR_EXCEPTION(
00314     max_nz_load <= 0 || max_nz_load > max_nz_ - nz_, std::invalid_argument
00315     ,"MatrixSparseCOORSerial::get_load_nonzeros_buffers(...) : Error, "
00316     "The number of nonzeros to load max_nz_load = " << max_nz_load << " can not "
00317     "be greater than max_nz - nz = " << max_nz_ << " - " << nz_ << " = " << (max_nz_-nz_) <<
00318     " entries!" );
00319   TEUCHOS_TEST_FOR_EXCEPTION(
00320     reload_val_only_ && (row_i != NULL || col_j != NULL), std::invalid_argument
00321     ,"MatrixSparseCOORSerial::get_load_nonzeros_buffers(...) : Error, "
00322     "reset_to_load_values() was called and therefore the structure of the matrix "
00323     "can not be set!" );
00324   TEUCHOS_TEST_FOR_EXCEPTION(
00325     !reload_val_only_  && (row_i == NULL || col_j == NULL), std::invalid_argument
00326     ,"MatrixSparseCOORSerial::get_load_nonzeros_buffers(...) : Error, "
00327     "both *row_i and *col_j must be non-NULL since reinitialize() was called" );
00328 #endif
00329   max_nz_load_ = max_nz_load;
00330   *val   = val_   + nz_;
00331   if(!reload_val_only_)
00332     *row_i = row_i_ + nz_;
00333   if(!reload_val_only_)
00334     *col_j = col_j_ + nz_;
00335 }
00336 
00337 void MatrixSparseCOORSerial::commit_load_nonzeros_buffers(
00338   size_type      nz_commit
00339   ,value_type    **val
00340   ,index_type    **row_i
00341   ,index_type    **col_j
00342   )
00343 {
00344 #ifdef TEUCHOS_DEBUG
00345   TEUCHOS_TEST_FOR_EXCEPTION(
00346     max_nz_load_ == 0 , std::logic_error
00347     ,"MatrixSparseCOORSerial::commit_load_nonzeros_buffers(...) : Error, "
00348     "You must call get_load_nonzeros_buffers() before calling this method!" );
00349   TEUCHOS_TEST_FOR_EXCEPTION(
00350     nz_commit > max_nz_load_ , std::logic_error
00351     ,"MatrixSparseCOORSerial::commit_load_nonzeros_buffers(...) : Error, "
00352     "You can not commit more nonzero entries than you requested buffer space for in "
00353     "get_load_nonzeros_buffers(...)!" );
00354   TEUCHOS_TEST_FOR_EXCEPTION(
00355     *val != val_ + nz_
00356     , std::logic_error
00357     ,"MatrixSparseCOORSerial::commit_load_nonzeros_buffers(...) : Error, "
00358     "This is not the buffer I give you in get_load_nonzeros_buffers(...)!" );
00359   TEUCHOS_TEST_FOR_EXCEPTION(
00360     reload_val_only_ && (row_i != NULL || col_j != NULL), std::invalid_argument
00361     ,"MatrixSparseCOORSerial::commit_load_nonzeros_buffers(...) : Error, "
00362     "reset_to_load_values() was called and therefore the structure of the matrix "
00363     "can not be set!" );
00364 #endif
00365   nz_            += nz_commit;
00366   max_nz_load_    = 0;
00367 }
00368 
00369 void MatrixSparseCOORSerial::finish_construction( bool test_setup )
00370 {
00371   TEUCHOS_TEST_FOR_EXCEPTION(
00372     reload_val_only_ == true && reload_val_only_nz_last_ != nz_, std::logic_error
00373     ,"MatrixSparseCOORSerial::finish_construction() : Error, the number of nonzeros on"
00374     " the initial load with row and column indexes was = " << reload_val_only_nz_last_ <<
00375     " and does not agree with the number of nonzero values = " << nz_ << " loaded this time!" );
00376   if( test_setup ) {
00377     for( size_type k = 0; k < nz_; ++k ) {
00378       const index_type
00379         i = row_i_[k],
00380         j = col_j_[k];
00381       TEUCHOS_TEST_FOR_EXCEPTION(
00382         i < 1 || rows_ < i, std::logic_error
00383         ,"MatrixSparseCOORSerial::finish_construction(true) : Error, "
00384         "row_i[" << k << "] = " << i << " is not in the range [1,rows] = [1,"<<rows_<<"]!" );
00385       TEUCHOS_TEST_FOR_EXCEPTION(
00386         j < 1 || cols_ < j, std::logic_error
00387         ,"MatrixSparseCOORSerial::finish_construction(true) : Error, "
00388         "col_j[" << k << "] = " << j << " is not in the range [1,cols] = [1,"<<cols_<<"]!" );
00389     }
00390   }
00391   space_cols_.initialize(rows_);
00392   space_rows_.initialize(cols_);
00393 }
00394 
00395 // Overridden from MatrixExtractSparseElements
00396 
00397 #ifdef TEUCHOS_DEBUG
00398 #define VALIDATE_ROW_COL_IN_RANGE() \
00399 TEUCHOS_TEST_FOR_EXCEPTION( \
00400   i < 1 || rows_ < i, std::invalid_argument \
00401   ,err_msg_head<<", i = inv_row_perm[(row_i["<<k<<"]=="<<*row_i<<")-1] = "<<i<<" > rows = "<<rows_ ); \
00402 TEUCHOS_TEST_FOR_EXCEPTION( \
00403   j < 1 || cols_ < j, std::invalid_argument \
00404   ,err_msg_head<<", j = inv_col_perm[(col_j["<<k<<"]=="<<*col_j<<")-1] = "<<j<<" > rows = "<<cols_ );
00405 #else
00406 #define VALIDATE_ROW_COL_IN_RANGE()
00407 #endif
00408 
00409 index_type MatrixSparseCOORSerial::count_nonzeros(
00410   EElementUniqueness    element_uniqueness
00411   ,const index_type     inv_row_perm[]
00412   ,const index_type     inv_col_perm[]
00413   ,const Range1D        &row_rng_in
00414   ,const Range1D        &col_rng_in
00415   ,index_type           dl
00416   ,index_type           du
00417   ) const
00418 {
00419 #ifdef TEUCHOS_DEBUG
00420   const char err_msg_head[] = "MatrixSparseCOORSerial::count_nonzeros(...): Error";
00421   TEUCHOS_TEST_FOR_EXCEPTION(
00422     element_uniqueness_ == ELEMENTS_ASSUME_DUPLICATES_SUM && element_uniqueness == ELEMENTS_FORCE_UNIQUE
00423     ,std::logic_error
00424     ,err_msg_head << ", the client requests a count for unique "
00425     "elements but this sparse matrix object is not allowed to assume this!" );
00426 #endif
00427   const Range1D
00428     row_rng = RangePack::full_range(row_rng_in,1,rows_),
00429     col_rng = RangePack::full_range(col_rng_in,1,rows_),
00430     row_rng_full(1,rows_),
00431     col_rng_full(1,cols_);
00432   const index_type
00433     *row_i    = row_i_,
00434     *col_j    = col_j_;
00435   index_type
00436     cnt_nz    = 0,
00437     k         = 0;
00438   if( dl == -row_rng.ubound() + col_rng.lbound() && du == +col_rng.ubound() - row_rng.lbound() ) {
00439     // The diagonals are not limiting so we can ignore them
00440     if( row_rng == row_rng_full && col_rng == col_rng_full ) {
00441       // The row and column ranges are not limiting either
00442       cnt_nz = nz_; // Just return the count of all the elements!
00443     }
00444     else {
00445       // The row or column range is limiting
00446       if( inv_row_perm == NULL && inv_col_perm == NULL ) {
00447         // Neither the rows nor columns are permuted
00448         for( k = 0; k < nz_; ++row_i, ++col_j, ++k ) {
00449           const index_type
00450             i = (*row_i),
00451             j = (*col_j);
00452           VALIDATE_ROW_COL_IN_RANGE();
00453           cnt_nz += row_rng.in_range(i) && col_rng.in_range(j) ? 1 : 0;
00454         }
00455       }
00456       else if ( inv_row_perm != NULL && inv_col_perm == NULL ) {
00457         // Only the rows are permuted 
00458         for( k = 0; k < nz_; ++row_i, ++col_j, ++k ) {
00459           const index_type
00460             i = inv_row_perm[(*row_i)-1],
00461             j = (*col_j);
00462           VALIDATE_ROW_COL_IN_RANGE();
00463           cnt_nz += row_rng.in_range(i) && col_rng.in_range(j) ? 1 : 0;
00464         }
00465       }
00466       else if ( inv_row_perm == NULL && inv_col_perm != NULL ) {
00467         // Only the columns are permuted 
00468         for( k = 0; k < nz_; ++row_i, ++col_j, ++k ) {
00469           const index_type
00470             i = (*row_i),
00471             j = inv_col_perm[(*col_j)-1];
00472           VALIDATE_ROW_COL_IN_RANGE();
00473           cnt_nz += row_rng.in_range(i) && col_rng.in_range(j) ? 1 : 0;
00474         }
00475       }
00476       else {
00477         // Both the rows and columns are permuted!
00478         for( k = 0; k < nz_; ++row_i, ++col_j, ++k ) {
00479           const index_type
00480             i = inv_row_perm[(*row_i)-1],
00481             j = inv_col_perm[(*col_j)-1];
00482           VALIDATE_ROW_COL_IN_RANGE();
00483           cnt_nz += row_rng.in_range(i) && col_rng.in_range(j) ? 1 : 0;
00484         }
00485       }
00486     }
00487   }
00488   else {
00489     // We have to consider the diagonals dl and du
00490     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00491   }
00492   return cnt_nz;
00493 }
00494 
00495 void MatrixSparseCOORSerial::coor_extract_nonzeros(
00496   EElementUniqueness    element_uniqueness
00497   ,const index_type     inv_row_perm[]
00498   ,const index_type     inv_col_perm[]
00499   ,const Range1D        &row_rng_in
00500   ,const Range1D        &col_rng_in
00501   ,index_type           dl
00502   ,index_type           du
00503   ,value_type           alpha
00504   ,const index_type     len_Aval
00505   ,value_type           Aval[]
00506   ,const index_type     len_Aij
00507   ,index_type           Arow[]
00508   ,index_type           Acol[]
00509   ,const index_type     row_offset
00510   ,const index_type     col_offset
00511   ) const
00512 {
00513 #ifdef TEUCHOS_DEBUG
00514   const char err_msg_head[] = "MatrixSparseCOORSerial::count_nonzeros(...): Error";
00515   TEUCHOS_TEST_FOR_EXCEPTION(
00516     element_uniqueness_ == ELEMENTS_ASSUME_DUPLICATES_SUM && element_uniqueness == ELEMENTS_FORCE_UNIQUE
00517     ,std::logic_error
00518     ,err_msg_head << ", the client requests extraction of unique "
00519     "elements but this sparse matrix object can not guarantee this!" );
00520 #endif
00521   const Range1D
00522     row_rng = RangePack::full_range(row_rng_in,1,rows_),
00523     col_rng = RangePack::full_range(col_rng_in,1,rows_),
00524     row_rng_full(1,rows_),
00525     col_rng_full(1,cols_);
00526   value_type
00527     *val      = val_;
00528   const index_type
00529     *row_i    = row_i_,
00530     *col_j    = col_j_;
00531   index_type
00532     cnt_nz    = 0,
00533     k         = 0;
00534   if( dl == -row_rng.ubound() + col_rng.lbound() && du == +col_rng.ubound() - row_rng.lbound() ) {
00535     // The diagonals are not limiting so we can ignore them
00536     if( row_rng == row_rng_full && col_rng == col_rng_full ) {
00537       // The row and column ranges are not limiting either
00538       if( inv_row_perm == NULL && inv_col_perm == NULL ) {
00539         // We are just extracting the whole, unpermuted matrix
00540         for( k = 0; k < nz_; ++val, ++row_i, ++col_j, ++k ) {
00541           ++cnt_nz;
00542           if( len_Aval )
00543             *Aval++ = *val;           // ToDo: Split into different loops (no inner if())
00544           if( len_Aij ) {
00545             *Arow++ = *row_i + row_offset;
00546             *Acol++ = *col_j + col_offset;
00547           }
00548         }
00549       }
00550       else {
00551         TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00552       }
00553     }
00554     else {
00555       // The row or column range is limiting
00556       if( inv_row_perm == NULL && inv_col_perm == NULL ) {
00557         // There are no permutations to consider
00558         for( k = 0; k < nz_; ++val, ++row_i, ++col_j, ++k ) {
00559           const index_type
00560             i = (*row_i),
00561             j = (*col_j);
00562           VALIDATE_ROW_COL_IN_RANGE();
00563           if( row_rng.in_range(i) && col_rng.in_range(j) ) {
00564             ++cnt_nz;
00565             if( len_Aval )
00566               *Aval++ = *val;           // ToDo: Split into different loops (no inner if())
00567             if( len_Aij ) {
00568               *Arow++ = i + row_offset;
00569               *Acol++ = j + col_offset;
00570             }
00571           }
00572         }
00573       }
00574       else if( inv_row_perm != NULL && inv_col_perm == NULL ) {
00575         // We must consider only row permutations
00576         for( k = 0; k < nz_; ++val, ++row_i, ++col_j, ++k ) {
00577           const index_type
00578             i = inv_row_perm[(*row_i)-1],
00579             j = (*col_j);
00580           VALIDATE_ROW_COL_IN_RANGE();
00581           if( row_rng.in_range(i) && col_rng.in_range(j) ) {
00582             ++cnt_nz;
00583             if( len_Aval )
00584               *Aval++ = *val;           // ToDo: Split into different loops (no inner if())
00585             if( len_Aij ) {
00586               *Arow++ = i + row_offset;
00587               *Acol++ = j + col_offset;
00588             }
00589           }
00590         }
00591       }
00592       else if( inv_row_perm == NULL && inv_col_perm != NULL ) {
00593         // We must consider only column permutations
00594         for( k = 0; k < nz_; ++val, ++row_i, ++col_j, ++k ) {
00595           const index_type
00596             i = (*row_i),
00597             j = inv_col_perm[(*col_j)-1];
00598           VALIDATE_ROW_COL_IN_RANGE();
00599           if( row_rng.in_range(i) && col_rng.in_range(j) ) {
00600             ++cnt_nz;
00601             if( len_Aval )
00602               *Aval++ = *val;           // ToDo: Split into different loops (no inner if())
00603             if( len_Aij ) {
00604               *Arow++ = i + row_offset;
00605               *Acol++ = j + col_offset;
00606             }
00607           }
00608         }
00609       }
00610       else {
00611         // We must consider row and column permutations
00612         for( k = 0; k < nz_; ++val, ++row_i, ++col_j, ++k ) {
00613           const index_type
00614             i = inv_row_perm[(*row_i)-1],
00615             j = inv_col_perm[(*col_j)-1];
00616           VALIDATE_ROW_COL_IN_RANGE();
00617           if( row_rng.in_range(i) && col_rng.in_range(j) ) {
00618             ++cnt_nz;
00619             if( len_Aval )
00620               *Aval++ = *val;           // ToDo: Split into different loops (no inner if())
00621             if( len_Aij ) {
00622               *Arow++ = i + row_offset;
00623               *Acol++ = j + col_offset;
00624             }
00625           }
00626         }
00627       }
00628     }
00629   }
00630   else {
00631     // We have to consider the diagonals dl and du
00632     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement!
00633   }
00634   TEUCHOS_TEST_FOR_EXCEPT( !(  len_Aval == 0 || len_Aval == cnt_nz  ) );
00635   TEUCHOS_TEST_FOR_EXCEPT( !(  len_Aij == 0  || len_Aij  == cnt_nz  ) );
00636 }
00637 
00638 // private
00639 
00640 void MatrixSparseCOORSerial::make_storage_unique()
00641 {
00642   if( release_resource_.count() > 1 ) {
00643     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Allocate new storage and copy this memory.
00644     self_allocate_ = true;
00645   }
00646 }
00647 
00648 } // end namespace AbstractLinAlgPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends