daily pastebin goal
42%
SHARE
TWEET

Untitled

a guest Dec 7th, 2017 50 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # include <initializer_list>
  2. # include <iomanip>
  3. # include <iostream>
  4. # include <cassert>
  5. # include <cmath>
  6. # include <vector>
  7. # include <string>
  8. # include <set>
  9. # include <algorithm>
  10. # include <map>
  11. # include <exception>
  12.  
  13. # define __TESTING__
  14.  
  15. template<typename T>
  16. class DIAmatrix;
  17.  
  18.  
  19. template <typename U>
  20. std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m );
  21.  
  22. //- SpMV
  23. //
  24. template<typename U>
  25. std::vector<U> operator*(const DIAmatrix<U>& , const std::vector<U>& ) ;
  26.  
  27.  
  28. template<typename Type >
  29. class DIAmatrix
  30. {
  31.     template <typename U>
  32.     friend std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m );
  33.  
  34.      //- SpMV
  35.      //
  36.      template<typename U>
  37.      friend std::vector<U> operator*(const DIAmatrix<U>& , const std::vector<U>& ) ;
  38.  
  39.  
  40.   public:
  41.  
  42.      constexpr DIAmatrix(std::initializer_list<std::vector<Type>> ) ;
  43.  
  44.      Type& operator()(const std::size_t , const std::size_t ) noexcept ;
  45.  
  46.       auto constexpr printDIA() const noexcept ;
  47.  
  48.  
  49.     private:
  50.  
  51.       std::size_t denseRows ;
  52.       std::size_t denseCols ;
  53.  
  54.       mutable Type dummy    ;
  55.  
  56.       std::size_t nnz ;
  57.  
  58.       std::size_t dim ;      
  59.  
  60.       std::map<int ,  std::vector<Type>> value ;   // value
  61.       std::set<int>                        dig ;   // off-diagonal distance
  62.  
  63.       Type constexpr findValue(const std::size_t , const std::size_t ) const noexcept ;
  64. };
  65.  
  66. //-------       implemantation      ------- //
  67.  
  68. template<typename T>
  69. constexpr DIAmatrix<T>::DIAmatrix(std::initializer_list<std::vector<T>> rows )  
  70. {
  71.     denseRows = rows.size();
  72.     denseCols = (*rows.begin()).size();
  73.  
  74.     if(denseRows != denseCols )
  75.     {
  76.        throw std::runtime_error("DIA-Matrix , SIZE EXCEPTION THROWN :n>>> Matrix Must be square! <<<");  
  77.     }
  78.  
  79.     dim = denseRows ;
  80.  
  81.     value[0].resize(dim);
  82.     value[1].resize(dim-1);       // diagonal
  83.     value[-1].resize(dim-1);      // tri-bands as default
  84.  
  85.     dig.insert(-1);
  86.     dig.insert(0);                // tribands as default
  87.     dig.insert(1);
  88.  
  89.  
  90.     auto i =0 , j=0 ;  
  91.     for(auto row : rows )
  92.     {
  93.        j=0;    
  94.        for(auto val : row )
  95.        {
  96.          if(val != 0)
  97.          {
  98.                 if( i==j )   value[0].at(i) = val  ;
  99.            else if(j==i-1)   value[-1].at(j) = val  ;
  100.            else if(j==i+1)   value[1].at(i) = val ;
  101.            else if(j < i )
  102.            {
  103.               auto offset =  j-i ;
  104.               if( !(dig.find(offset) != dig.end()) )
  105.               {
  106.                   dig.insert(offset);
  107.                   value[offset].resize(dim-(fabs(j-i)));
  108.                   value[offset].at(j) = val ;
  109.               }
  110.               else
  111.                   value[offset].at(j) = val ;
  112.            }
  113.            else if(j > i)
  114.            {
  115.              auto offset = j-i ;
  116.              if( !(dig.find(offset) != dig.end()))
  117.              {    
  118.                  dig.insert(offset);
  119.                  value[offset].resize(dim-(fabs(j-i)));
  120.                  value[offset].at(i) = val ;
  121.              }
  122.              else
  123.                value[offset].at(i) = val;
  124.            }
  125.          }
  126.          j++;
  127.        }
  128.        i++;
  129.     }
  130. #  ifdef __TESTING__
  131.       printDIA();
  132. #  endif
  133. }
  134.  
  135.  
  136. template<typename T>
  137. T constexpr DIAmatrix<T>::findValue(const std::size_t r, const std::size_t c) const noexcept
  138. {
  139.     T val = 0.0;  
  140.     if( value.find(c-r) != value.end() )  
  141.        val  = r < c ? value.at(c-r).at(r) : value.at(c-r).at(c) ;      
  142.     return val ;    
  143. }
  144.  
  145.  
  146. template<typename T>
  147. auto constexpr DIAmatrix<T>::printDIA() const noexcept
  148. {
  149.  
  150.    auto it = dig.begin() ;
  151.  
  152.    for( auto& n : dig)
  153.     {
  154.       std::cout << std::setw(4) <<*it << " | :  " ;
  155.       for(auto& x : value.at(n) )
  156.       {
  157.          std::cout << x << " " ;  
  158.       }
  159.       ++it ;
  160.       std::cout << std::endl;
  161.     }
  162. }
  163.  
  164. template <typename T>
  165. T& DIAmatrix<T>::operator()(const std::size_t r, const std::size_t c) noexcept
  166. {
  167.       assert( r >0 && r <= dim && c > 0 && c <= dim ) ;
  168.       dummy = findValue(r-1,c-1);
  169.       return dummy ;
  170. }
  171.  
  172.  
  173. template <typename U>
  174. std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m )
  175. {
  176.       for(auto i=1; i <= m.denseRows ; i++ )
  177.       {
  178.           for(auto j=1 ; j <= m.denseCols ; j++)
  179.           {
  180.                os << std::setw(6) << m(i,j) << " ";  
  181.           }
  182.           os << std::endl;  
  183.       }
  184. }
  185.  
  186.  
  187.  
  188. template<typename T>
  189. std::vector<T> operator*(const DIAmatrix<T>& m, const std::vector<T>& x )
  190. {
  191.     std::vector<T> y(x.size());
  192.  
  193.     if(m.denseCols != x.size())
  194.     {
  195.        throw std::runtime_error("Erorr size of matrix and vector doesn't match!");
  196.     }
  197.     else
  198.     {
  199.       // for(auto i=0; i< x.size() ; i++ )
  200.           for(auto offset : m.dig )
  201.           {
  202.             for(auto j = 0 ; j < m.value.at(offset).size() ; j++ )
  203.             {
  204.                  y.at(j) += m.value.at(offset).at(j) * x.at(j+fabs(offset))  ;
  205.             }
  206.          }
  207.       // }  
  208.        //}
  209.     }
  210.     return y;
  211. }
  212.    
  213. # include "DiaMatrix.H"
  214.  
  215. using namespace std;
  216.  
  217. int main() {
  218.  
  219.    DIAmatrix<int> dm = {{11,12,0,0,0,0,0,0} ,{0,22,0,0,0,0,0,0} ,{31,32,33,0,0,0,0,0},
  220.                               {41,42,43,44,0,0,0,0}, {0,0,0,0,55,56,0,0},{0,0,0,0,0,66,67,0},{0,0,0,0,0,0,77,78},{0,0,0,0,0,0,87,88}};
  221.  
  222.      cout << dm ;
  223.  
  224.     std::vector<int> v1 = {3,4,0,1,6,8,1,19};
  225.  
  226.  
  227.     vector<int> v2 = dm * v1 ;
  228.  
  229.    for(auto& x : v2)
  230.      cout << x << ' ' ;
  231.    cout << endl;  
  232.  
  233.  
  234.   return 0;
  235. }
  236.    
  237. 122 130 43 44 778 595 3212 1672
  238.    
  239. 81 88 221 335 778 595 1559 1759
RAW Paste Data
Top