daily pastebin goal
49%
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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top