AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects Version of the Day
AbstractLinAlgPack_COOMatrixPartitionedViewClassDef.hpp
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 #ifndef COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
00043 #define COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
00044 
00045 #include <sstream>
00046 #include <algorithm>
00047 #include <functional>
00048 
00049 #include "AbstractLinAlgPack_COOMatrixPartitionedViewClassDecl.hpp"
00050 
00051 namespace AbstractLinAlgPack {
00052 
00053 // /////////////////////////////////////////////////////////////////////////////
00054 // Template function definitions for COOMatrixPartitionedView
00055 
00056 template <class T_Indice, class T_Value>
00057 void COOMatrixPartitionedView<T_Indice,T_Value>::create_view(
00058         size_type       rows
00059       , size_type       cols
00060       , size_type       nz    
00061       , value_type      val[]
00062       , const indice_type   ivect[]
00063       , const indice_type   jvect[]
00064       , const size_type   inv_row_perm[]
00065       , const size_type   inv_col_perm[]
00066       , const size_type   num_row_part
00067       , const size_type   row_part[]
00068       , const size_type   num_col_part
00069       , const size_type   col_part[]
00070       , const EPartitionOrder partition_order )
00071 {
00072   try {
00073 
00074   // 1) Check some preconditins before we start
00075 
00076   // Check the sparsisty density of COO matrix
00077   if(nz > rows * cols)
00078     throw std::out_of_range(
00079       "COOMatrixPartitionedView<...>::create_view() : Input error, "
00080       "nz can not be greater than rows * cols");
00081   
00082   // Check row and column partition information
00083 
00084   if(num_row_part > rows || num_col_part > rows)
00085     throw std::out_of_range(
00086       "COOMatrixPartitionedView<...>::create_view() : Input error, "
00087       "num_rows_part and num_col_part can not be greater than"
00088       " rows and cols respectivly");
00089 
00090   if(row_part[num_row_part] > rows + 1)
00091     throw std::out_of_range(
00092       "COOMatrixPartitionedView<...>::create_view() : Input error, "
00093       "row_part[num_row_part] can not be greater than rows");
00094 
00095   if(col_part[num_col_part] > cols + 1)
00096     throw std::out_of_range(
00097       "COOMatrixPartitionedView<...>::create_view() : Input error, "
00098       "col_part[num_col_part] can not be greater than cols");
00099 
00100   {for(size_type i = 1; i < num_row_part + 1; ++i)
00101     if(row_part[i-1] >= row_part[i])
00102       throw std::domain_error(
00103         "COOMatrixPartitionedView<...>::create_view() : Input error, "
00104         "row_part[i-1] < row_part[i] does not hold");}
00105 
00106   {for(size_type i = 1; i < num_col_part + 1; ++i)
00107     if(col_part[i-1] >= col_part[i])
00108       throw std::domain_error(
00109         "COOMatrixPartitionedView<...>::create_view() : Input error, "
00110         "col_part[i-1] < col_part[i] does not hold");}
00111                 
00112   // Get references to referenced quantities
00113   std::vector<size_type>
00114     &_row_part    = ref_row_part_.obj(),
00115     &_col_part    = ref_col_part_.obj(),
00116     &_part_start  = ref_part_start_.obj();
00117   ele_type
00118     &_ele     = ref_ele_.obj();
00119 
00120   // 2) Initialize storage for data members and copy data
00121   num_row_part_   = num_row_part;
00122   num_col_part_   = num_col_part;
00123   _row_part.resize(num_row_part_ + 1);
00124   std::copy(row_part, row_part+ num_row_part_+1, _row_part.begin());
00125   _col_part.resize(num_col_part_ + 1);
00126   std::copy(col_part, col_part+ num_col_part_+1, _col_part.begin());
00127   partition_order_  = partition_order;
00128   _ele.resize(nz,element_type()); // hack to get around compiler error.
00129   _part_start.resize(num_row_part_ * num_col_part_+1);
00130   _part_start.assign(_part_start.size(),0); // set to 0
00131 
00132   // 3) Count the number of nonzero elements in each overall partition
00133   {
00134     // Use the storage locations _part_start[1] ... _part_start[num_part]
00135     // to count the number of nonzero elements in each overall partition.
00136     // In particular num_in_part[overall_p - 1] will give the number
00137     // of nonzero elements in the overall partition number overall_p.
00138     //
00139     // _part_start = { 0, nz1, nz2,...,nzn } 
00140     //
00141     size_type *num_in_part = &_part_start[0] + 1; 
00142 
00143     // Loop (time = O(nz), space = O(1))
00144 
00145     // read iterators
00146     const indice_type *ivect_itr    = ivect,
00147               *ivect_itr_end  = ivect + nz,
00148               *jvect_itr    = jvect;
00149 
00150     for(;ivect_itr != ivect_itr_end; ++ivect_itr, ++jvect_itr) {
00151       // Get row and column indices in the non-permuted matrix
00152       indice_type   i_org = *ivect_itr,
00153               j_org = *jvect_itr;
00154       // assert that they are in range
00155       if(i_org < 1 || i_org > rows || j_org < 1 || j_org > cols) {
00156         std::ostringstream omsg;
00157         omsg  <<  "COOMatrixPartitionedView<...>::create_view() : "
00158               " Error, element k = " << ivect_itr - ivect
00159             <<  " in the non-permuted matrix"
00160               " is out of range with rows = " << rows
00161             <<  ", cols = " << cols << ", i = " << i_org
00162             <<  ", and j = " << j_org;
00163         throw std::out_of_range(omsg.str());
00164       }
00165       // Get row and column indices in the permuted matrix
00166       indice_type   i = inv_row_perm[i_org - 1],
00167               j = inv_col_perm[j_org - 1];
00168       // assert that they are in range
00169       if(i < 1 || i > rows || j < 1 || j > cols) {
00170         std::ostringstream omsg;
00171         omsg  <<  "COOMatrixPartitionedView<...>::create_view() : "
00172               " Error, element k = " << ivect_itr - ivect
00173             <<  " in the permuted matrix"
00174               " is out of range with rows = " << rows
00175             <<  ", cols = " << cols << ", i = " << i
00176             <<  ", and j = " << j;
00177         throw std::out_of_range(omsg.str());
00178       }
00179       // get the overall partition number
00180       size_type overall_p = overall_p_from_ij(_row_part,_col_part,i,j);
00181       // Increment the number of nonzero elements.
00182       num_in_part[overall_p - 1]++;
00183     }
00184   }
00185 
00186   // 4) Set part_start[ovarall_p] equal to the start in ele
00187   // for the nonzero elements in that partition.
00188   //
00189   // _part_start = { start_1 = 0, start_2,..., start_n, ###}
00190   //
00191   {for(size_type i = 2; i < num_row_part_ * num_col_part_; ++i)
00192     _part_start[i] += _part_start[i-1];}
00193 
00194   // 5) Shift the elements over
00195   //
00196   // _part_start = { 0, start_1 = 0, start_2,..., start_n}
00197   //
00198   {for(size_type i = num_row_part_ * num_col_part_; i > 0; --i)
00199     _part_start[i] = _part_start[i-1];}
00200 
00201   // 6) Add the nonzero elements to each partition.  When we
00202   // are done we should have _part_start initialized properly.
00203   //
00204   // part_start = { start_1 = 0, start_2,..., start_n, total_nz }
00205   // 
00206   {
00207     // next_ele_insert[overall_p - 1] is the possition in ele
00208     // for the next element to ensert
00209     size_type *next_ele_insert = &_part_start[0] + 1;
00210 
00211     // Loop (time = O(nz), space = O(1))
00212 
00213     // read iterators
00214     value_type      *val_itr    = val,
00215               *val_itr_end  = val + nz;
00216     const indice_type *ivect_itr    = ivect,
00217               *jvect_itr    = jvect;
00218 
00219     for(;val_itr != val_itr_end; ++val_itr, ++ivect_itr, ++jvect_itr) {
00220       // Get row and column indices in the permuted matrix
00221       indice_type   i = inv_row_perm[*ivect_itr - 1],
00222               j = inv_col_perm[*jvect_itr - 1];
00223       // get the overall partition number
00224       size_type overall_p = overall_p_from_ij(_row_part,_col_part,i,j);
00225       // Add the element to the partition
00226       _ele[next_ele_insert[overall_p - 1]++].initialize(val_itr,i,j);
00227     }
00228   }
00229 
00230   } // end try
00231   catch(...) {
00232     free();
00233     throw;  // rethrow the exception out of here
00234   }
00235 }
00236 
00237 template <class T_Indice, class T_Value>
00238 void COOMatrixPartitionedView<T_Indice,T_Value>::bind(const COOMatrixPartitionedView& coo_view)
00239 {
00240   num_row_part_   = coo_view.num_row_part_;
00241   num_col_part_   = coo_view.num_col_part_;
00242   ref_row_part_   = coo_view.ref_row_part_;
00243   ref_col_part_   = coo_view.ref_col_part_;
00244   partition_order_  = coo_view.partition_order_;
00245   ref_ele_      = coo_view.ref_ele_;
00246   ref_part_start_   = coo_view.ref_part_start_;
00247 }
00248 
00249 template <class T_Indice, class T_Value>
00250 void COOMatrixPartitionedView<T_Indice,T_Value>::free()
00251 {
00252   // Reinitialize to uninitizlied
00253   num_row_part_ = num_col_part_ = 0;
00254   if(ref_row_part_.has_ref_set())   ref_row_part_.obj().resize(0);
00255   if(ref_col_part_.has_ref_set())   ref_col_part_.obj().resize(0);
00256   if(ref_ele_.has_ref_set())      ref_ele_.obj().resize(0,element_type());
00257   if(ref_part_start_.has_ref_set()) ref_part_start_.obj().resize(0);
00258 }
00259 
00260 template <class T_Indice, class T_Value>
00261 void COOMatrixPartitionedView<T_Indice,T_Value>::get_row_part(indice_type row_part[]) const
00262 {
00263   assert_initialized();
00264   const std::vector<size_type> &_row_part = ref_row_part_.const_obj();
00265   std::copy(_row_part.begin(), _row_part.end(), row_part);
00266 }
00267 
00268 template <class T_Indice, class T_Value>
00269 void COOMatrixPartitionedView<T_Indice,T_Value>::get_col_part(indice_type col_part[]) const
00270 {
00271   assert_initialized();
00272   const std::vector<size_type> &_col_part = ref_col_part_.const_obj();
00273   std::copy(_col_part.begin(), _col_part.end(), col_part);
00274 }
00275 
00276 template <class T_Indice, class T_Value>
00277 COOMatrixPartitionedView<T_Indice,T_Value>::size_type
00278 COOMatrixPartitionedView<T_Indice,T_Value>::part_num(const vector_size_type& part
00279   , size_type indice)
00280 {
00281   return std::upper_bound(part.begin(),part.end(),indice) - part.begin();
00282 }
00283 
00284 template <class T_Indice, class T_Value>
00285 COOMatrixPartitionedView<T_Indice,T_Value>::partition_type
00286 COOMatrixPartitionedView<T_Indice,T_Value>::create_partition(Range1D rng_overall_p) const
00287 {
00288   assert_initialized();
00289   // get reference to data structures
00290   const std::vector<size_type>
00291     &row_part = ref_row_part_.const_obj(),
00292     &col_part = ref_col_part_.const_obj(),
00293     &part_start = ref_part_start_.const_obj();
00294   ele_type
00295     &_ele   = const_cast<ref_ele_type&>(ref_ele_).obj(); // This is ok
00296   // Get upper and lower overall, row and column partition numbers
00297   rng_overall_p = DenseLinAlgPack::full_range(rng_overall_p,1,num_row_part_*num_col_part_);
00298   size_type l_p   = rng_overall_p.lbound(),
00299         u_p   = rng_overall_p.ubound(),
00300         l_r_p = imp_row_part_num(l_p),
00301         l_c_p = imp_col_part_num(l_p),
00302         u_r_p = imp_row_part_num(u_p),
00303         u_c_p = imp_col_part_num(u_p);
00304   // build argument list for creation of the partition.
00305   size_type
00306     rows  = row_part[u_r_p] - row_part[l_r_p - 1],
00307     cols  = col_part[u_c_p] - col_part[l_c_p - 1],
00308     nz    = part_start[u_p] - part_start[l_p - 1];
00309   element_type
00310     *ele    = &_ele[0] + part_start[l_p - 1];
00311   difference_type
00312     row_offset  = - (row_part[l_r_p - 1] - 1),
00313     col_offset  = - (col_part[l_c_p - 1] - 1);
00314 
00315   return partition_type(rows,cols,nz,ele,row_offset,col_offset);
00316 }
00317 
00318 } // end namespace AbstractLinAlgPack 
00319 
00320 #endif // COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends