AbstractLinAlgPack_COOMatrixPartitionedViewClassDef.hpp

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

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