ViennaCL - The Vienna Computing Library: /data/development/ViennaCL/dev/viennacl/linalg/row

Go to the documentation of this file.
00001 #ifndef VIENNACL_LINALG_ROW_SCALING_HPP_
00002 #define VIENNACL_LINALG_ROW_SCALING_HPP_
00003 
00004 /* =========================================================================
00005    Copyright (c) 2010-2011, Institute for Microelectronics,
00006                             Institute for Analysis and Scientific Computing,
00007                             TU Wien.
00008 
00009                             -----------------
00010                   ViennaCL - The Vienna Computing Library
00011                             -----------------
00012 
00013    Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
00014                
00015    (A list of authors and contributors can be found in the PDF manual)
00016 
00017    License:         MIT (X11), see file LICENSE in the base directory
00018 ============================================================================= */
00019 
00024 #include <vector>
00025 #include <cmath>
00026 #include "viennacl/forwards.h"
00027 #include "viennacl/vector.hpp"
00028 #include "viennacl/compressed_matrix.hpp"
00029 #include "viennacl/tools/tools.hpp"
00030 
00031 #include <map>
00032 
00033 namespace viennacl
00034 {
00035   namespace linalg
00036   {
00037     
00040     class row_scaling_tag
00041     {
00042       public:
00047         row_scaling_tag(unsigned int p = 2) : norm_(p) {}
00048         
00050         unsigned int norm() const { return norm_; }
00051         
00052       private:
00053         unsigned int norm_;
00054     };
00055     
00056 
00059     template <typename MatrixType>
00060     class row_scaling
00061     {
00062       typedef typename MatrixType::value_type      ScalarType;
00063       
00064       public:
00070         row_scaling(MatrixType const & mat, row_scaling_tag const & tag) : system_matrix(mat), tag_(tag)
00071         {
00072           assert(mat.size1() == mat.size2());
00073           diag_M_inv.resize(mat.size1());  //resize without preserving values
00074           
00075           for (typename MatrixType::const_iterator1 row_it = system_matrix.begin1();
00076                 row_it != system_matrix.end1();
00077                 ++row_it)
00078           {
00079             diag_M_inv[row_it.index1()];
00080             for (typename MatrixType::const_iterator2 col_it = row_it.begin();
00081                   col_it != row_it.end();
00082                   ++col_it)
00083             {
00084               if (tag_.norm() == 1)
00085                 diag_M_inv[col_it.index1()] += std::fabs(*col_it);
00086               else
00087                 diag_M_inv[col_it.index1()] += (*col_it) * (*col_it);
00088             }
00089             if (diag_M_inv[row_it.index1()] == 0)
00090               throw "ViennaCL: Zero row encountered while setting up row scaling preconditioner!";
00091             
00092             if (tag_.norm() == 1)
00093               diag_M_inv[row_it.index1()] = static_cast<ScalarType>(1.0) / diag_M_inv[row_it.index1()];
00094             else
00095               diag_M_inv[row_it.index1()] = static_cast<ScalarType>(1.0) / std::sqrt(diag_M_inv[row_it.index1()]);
00096           }
00097         }
00098         
00099         
00101         template <typename VectorType>
00102         void apply(VectorType & vec) const
00103         {
00104           assert(vec.size() == diag_M_inv.size());
00105           for (size_t i=0; i<vec.size(); ++i)
00106           {
00107             vec[i] *= diag_M_inv[i];
00108           }
00109         }
00110         
00111       private:
00112         MatrixType const & system_matrix;
00113         row_scaling_tag const & tag_;
00114         std::vector<ScalarType> diag_M_inv;
00115     };
00116 
00117     
00122     template <typename ScalarType, unsigned int MAT_ALIGNMENT>
00123     class row_scaling< compressed_matrix<ScalarType, MAT_ALIGNMENT> >
00124     {
00125       typedef compressed_matrix<ScalarType, MAT_ALIGNMENT>   MatrixType;
00126       
00127       public:
00133         row_scaling(MatrixType const & mat, row_scaling_tag const & tag) : system_matrix(mat), tag_(tag), diag_M_inv(mat.size1())
00134         {
00135           assert(system_matrix.size1() == system_matrix.size2());
00136           
00137           init_gpu();
00138         }
00139         
00140         /*
00141         void init_cpu()
00142         {
00143           std::vector< std::map<unsigned int, ScalarType> > cpu_check;
00144           std::vector<ScalarType> diag_M_inv_cpu(system_matrix.size1());
00145           
00146           copy(system_matrix, cpu_check);
00147           viennacl::tools::const_sparse_matrix_adapter<ScalarType> cpu_check_adapter(cpu_check);
00148           
00149           for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator1 row_it = cpu_check_adapter.begin1();
00150                 row_it != cpu_check_adapter.end1();
00151                 ++row_it)
00152           {
00153             diag_M_inv_cpu[row_it.index1()] = 0;
00154             for (typename viennacl::tools::const_sparse_matrix_adapter<ScalarType>::const_iterator2 col_it = row_it.begin();
00155                   col_it != row_it.end();
00156                   ++col_it)
00157             {
00158               if (tag_.norm() == 1)
00159                 diag_M_inv_cpu[col_it.index1()] += std::fabs(*col_it);
00160               else
00161                 diag_M_inv_cpu[col_it.index1()] += (*col_it) * (*col_it);
00162             }
00163             if (diag_M_inv_cpu[row_it.index1()] == 0)
00164               throw "ViennaCL: Zero row encountered while setting up row scaling preconditioner!";
00165             
00166             if (tag_.norm() == 1)
00167               diag_M_inv_cpu[row_it.index1()] = static_cast<ScalarType>(1.0) / diag_M_inv_cpu[row_it.index1()];
00168             else
00169               diag_M_inv_cpu[row_it.index1()] = static_cast<ScalarType>(1.0) / std::sqrt(diag_M_inv_cpu[row_it.index1()]);
00170           }
00171           
00172           diag_M_inv.resize(system_matrix.size1(), false);
00173           viennacl::fast_copy(diag_M_inv_cpu, diag_M_inv);
00174         } */
00175         
00176         void init_gpu()
00177         {
00178           if (tag_.norm() == 1)
00179           {
00180             viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(
00181                                                 viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
00182                                                 "row_scaling_1");
00183 
00184             viennacl::ocl::enqueue( k(system_matrix.handle1(), system_matrix.handle2(), system_matrix.handle(), 
00185                                       diag_M_inv, static_cast<cl_uint>(diag_M_inv.size())) );        
00186           }
00187           else
00188           {
00189             viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(
00190                                                 viennacl::linalg::kernels::compressed_matrix<ScalarType, MAT_ALIGNMENT>::program_name(),
00191                                                 "row_scaling_2");
00192 
00193             viennacl::ocl::enqueue( k(system_matrix.handle1(), system_matrix.handle2(), system_matrix.handle(), 
00194                                       diag_M_inv, static_cast<cl_uint>(diag_M_inv.size())) );        
00195           }
00196         }
00197         
00198         template <unsigned int ALIGNMENT>
00199         void apply(viennacl::vector<ScalarType, ALIGNMENT> & vec) const
00200         {
00201           assert(viennacl::traits::size1(system_matrix) == viennacl::traits::size(vec));
00202           
00203           //run kernel:
00204           viennacl::ocl::kernel & k = viennacl::ocl::get_kernel(viennacl::linalg::kernels::vector<ScalarType, ALIGNMENT>::program_name(),
00205                                                                 "diag_precond");
00206 
00207           viennacl::ocl::enqueue(
00208              k(viennacl::traits::handle(diag_M_inv), cl_uint(viennacl::traits::start(diag_M_inv)), cl_uint(viennacl::traits::size(diag_M_inv)),
00209                viennacl::traits::handle(vec), cl_uint(viennacl::traits::start(vec)), cl_uint(viennacl::traits::size(vec)) )
00210                                 );        
00211           
00212         }
00213         
00214       private:
00215         MatrixType const & system_matrix;
00216         row_scaling_tag const & tag_;
00217         viennacl::vector<ScalarType> diag_M_inv;
00218     };
00219 
00220   }
00221 }
00222 
00223 
00224 
00225 
00226 #endif
00227 
00228 
00229
/data/development/ViennaCL/dev/viennacl/linalg/row_scaling.hpp
