matrix.hpp 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185
  1. /**
  2. * \file matrix.hpp
  3. * \author Dhairya Malhotra, dhairya.malhotra@gmail.com
  4. * \date 2-11-2011
  5. * \brief This file contains definition of the class Matrix.
  6. */
  7. #include <stdint.h>
  8. #include <cstdlib>
  9. #include <pvfmm_common.hpp>
  10. #include <vector.hpp>
  11. #ifndef _PVFMM_MATRIX_HPP_
  12. #define _PVFMM_MATRIX_HPP_
  13. #ifdef __INTEL_OFFLOAD
  14. #pragma offload_attribute(push,target(mic))
  15. #endif
  16. namespace pvfmm{
  17. template <class T>
  18. class Permutation;
  19. template <class T>
  20. class Matrix{
  21. template <class Y>
  22. friend std::ostream& operator<<(std::ostream& output, const Matrix<Y>& M);
  23. public:
  24. struct
  25. Device{
  26. Device& operator=(Matrix& M){
  27. dim[0]=M.Dim(0);
  28. dim[1]=M.Dim(1);
  29. dev_ptr=(uintptr_t)M.data_ptr;
  30. return *this;
  31. }
  32. inline T* operator[](size_t j) const{
  33. assert(j<dim[0]);
  34. return &((T*)dev_ptr)[j*dim[1]];
  35. }
  36. size_t dim[2];
  37. uintptr_t dev_ptr;
  38. int lock_idx;
  39. };
  40. Matrix();
  41. Matrix(size_t dim1, size_t dim2, T* data_=NULL, bool own_data_=true);
  42. Matrix(const Matrix<T>& M);
  43. ~Matrix();
  44. void Swap(Matrix<T>& M);
  45. void ReInit(size_t dim1, size_t dim2, T* data_=NULL, bool own_data_=true);
  46. Device& AllocDevice(bool copy);
  47. void Device2Host(T* host_ptr=NULL);
  48. void Device2HostWait();
  49. void FreeDevice(bool copy);
  50. void Write(const char* fname);
  51. void Read(const char* fname);
  52. size_t Dim(size_t i) const;
  53. void Resize(size_t i, size_t j);
  54. void SetZero();
  55. Matrix<T>& operator=(const Matrix<T>& M);
  56. Matrix<T>& operator+=(const Matrix<T>& M);
  57. Matrix<T>& operator-=(const Matrix<T>& M);
  58. Matrix<T> operator+(const Matrix<T>& M2);
  59. Matrix<T> operator-(const Matrix<T>& M2);
  60. T& operator()(size_t i,size_t j) const;
  61. T* operator[](size_t i) const;
  62. Matrix<T> operator*(const Matrix<T>& M);
  63. static void GEMM(Matrix<T>& M_r, const Matrix<T>& A, const Matrix<T>& B, T beta=0.0);
  64. // cublasgemm wrapper
  65. static void CUBLASGEMM(Matrix<T>& M_r, const Matrix<T>& A, const Matrix<T>& B, T beta=0.0);
  66. void RowPerm(const Permutation<T>& P);
  67. void ColPerm(const Permutation<T>& P);
  68. Matrix<T> Transpose();
  69. static void Transpose(Matrix<T>& M_r, const Matrix<T>& M);
  70. // Original matrix is destroyed.
  71. void SVD(Matrix<T>& tU, Matrix<T>& tS, Matrix<T>& tVT);
  72. // Original matrix is destroyed.
  73. Matrix<T> pinv(T eps=-1);
  74. private:
  75. size_t dim[2];
  76. T* data_ptr;
  77. bool own_data;
  78. Device dev;
  79. Vector<char> dev_sig;
  80. #if defined(PVFMM_HAVE_CUDA)
  81. cudaEvent_t lock;
  82. #endif
  83. };
  84. template <class Y>
  85. std::ostream& operator<<(std::ostream& output, const Matrix<Y>& M);
  86. /**
  87. * /brief P=[e(p1)*s1 e(p2)*s2 ... e(pn)*sn],
  88. * where e(k) is the kth unit vector,
  89. * perm := [p1 p2 ... pn] is the permutation vector,
  90. * scal := [s1 s2 ... sn] is the scaling vector.
  91. */
  92. #define PERM_INT_T size_t
  93. template <class T>
  94. class Permutation{
  95. template <class Y>
  96. friend std::ostream& operator<<(std::ostream& output, const Permutation<Y>& P);
  97. public:
  98. Permutation(){}
  99. Permutation(size_t size);
  100. static Permutation<T> RandPerm(size_t size);
  101. Matrix<T> GetMatrix() const;
  102. size_t Dim() const;
  103. Permutation<T> Transpose();
  104. Permutation<T> operator*(const Permutation<T>& P);
  105. Matrix<T> operator*(const Matrix<T>& M);
  106. template <class Y>
  107. friend Matrix<Y> operator*(const Matrix<Y>& M, const Permutation<Y>& P);
  108. Vector<PERM_INT_T> perm;
  109. Vector<T> scal;
  110. };
  111. template <class T>
  112. Matrix<T> operator*(const Matrix<T>& M, const Permutation<T>& P);
  113. template <class Y>
  114. std::ostream& operator<<(std::ostream& output, const Permutation<Y>& P);
  115. }//end namespace
  116. #ifdef __INTEL_OFFLOAD
  117. #pragma offload_attribute(pop)
  118. #endif
  119. #include <matrix.txx>
  120. #endif //_PVFMM_MATRIX_HPP_