AbstractLinAlgPack_MatrixSparseCOORSerial.cpp

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

Generated on Thu Sep 18 12:35:12 2008 for MOOCHO (Single Doxygen Collection) by doxygen 1.3.9.1