Line data    Source code

       1             : #ifndef SET_MATRIX4_DIAGONAL_H
       2             : #define SET_MATRIX4_DIAGONAL_H
       3             : 
       4             : #include "Util/Util.H"
       5             : #include "Base.H"
       6             : 
       7             : namespace Set
       8             : {
       9             : template<>
      10             : class Matrix4<AMREX_SPACEDIM, Sym::Diagonal>
      11             : {
      12             :     Set::Matrix A = Set::Matrix::Ones();
      13             : public:
      14       85825 :     AMREX_GPU_HOST_DEVICE Matrix4() {};
      15       26574 :     AMREX_GPU_HOST_DEVICE Matrix4(Set::Matrix a_A): A(a_A) {};
      16             :     AMREX_FORCE_INLINE
      17             :         Scalar operator () (const int i, const int j, const int k, const int l) const
      18             :     {
      19        1940 :         if (i == k && j == l) return A(i, j);
      20        1680 :         else return 0.0;
      21             :     }
      22             :     void Randomize()
      23             :     {
      24             :         A = Set::Matrix::Random();
      25             :     }
      26             :     void Print(std::ostream& os)
      27             :     {
      28             :         os << "A = \n" << A;
      29             :     }
      30          27 :     static Matrix4<AMREX_SPACEDIM, Sym::Diagonal> Zero()
      31             :     {
      32          27 :         Matrix4<AMREX_SPACEDIM, Sym::Diagonal> zero;
      33          27 :         zero.A = Set::Matrix::Zero();
      34          27 :         return zero;
      35             :     }
      36             :     friend Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix& b);
      37             :     friend Set::Vector operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix3& b);
      38             :     friend bool operator == (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
      39             :     friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator + (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
      40             :     friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator - (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b);
      41             :     friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Scalar& b);
      42             :     friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Set::Scalar& b, const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a);
      43             :     friend Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator / (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Scalar& b);
      44             : 
      45             : 
      46             :     //AMREX_GPU_HOST_DEVICE void operator =  (Matrix4<AMREX_SPACEDIM,Sym::Isotropic> &a) {lambda =  a.lambda; mu =  a.mu;}
      47        8463 :     AMREX_GPU_HOST_DEVICE void operator += (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a) { A += a.A; }
      48             :     AMREX_GPU_HOST_DEVICE void operator -= (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a) { A -= a.A; }
      49             :     //AMREX_GPU_HOST_DEVICE void operator *= (const Matrix4<AMREX_SPACEDIM,Sym::Diagonal> &a) {A *= a.A;}
      50             :     //AMREX_GPU_HOST_DEVICE void operator /= (const Matrix4<AMREX_SPACEDIM,Sym::Diagonal> &a) {A /= a.A;}
      51             :     AMREX_GPU_HOST_DEVICE void operator *= (const Set::Scalar& alpha) { A *= alpha; }
      52             :     AMREX_GPU_HOST_DEVICE void operator /= (const Set::Scalar& alpha) { A /= alpha; }
      53             : 
      54             :     Set::Scalar Norm()
      55             :     {
      56             :         return A.lpNorm<2>();
      57             :     }
      58             :     bool contains_nan() const
      59             :     {
      60             :         if (std::isnan(A.lpNorm<2>())) return true;
      61             :         return false;
      62             :     }
      63             : 
      64             : };
      65             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      66             : Set::Matrix operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix& b)
      67             : {
      68     2443536 :     Set::Matrix ret;
      69             : 
      70     7332428 :     for (int i = 0; i < AMREX_SPACEDIM; i++)
      71             :     {
      72    14672136 :         for (int j = 0; j < AMREX_SPACEDIM; j++)
      73             :         {
      74     9783204 :             ret(i, j) = a.A(i, j) * b(i, j);
      75             :         }
      76             :     }
      77             : 
      78     2443536 :     return ret;
      79             : }
      80             : 
      81             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      82             : Set::Vector operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a, const Set::Matrix3& b)
      83             : {
      84      787937 :     Set::Vector ret = Set::Vector::Zero();
      85             : 
      86     2363810 :     for (int i = 0; i < AMREX_SPACEDIM; i++)
      87     4727620 :         for (int j = 0; j < AMREX_SPACEDIM; j++)
      88     6303500 :             ret(i) += a.A(i, j) * b(i, j, j);
      89             : 
      90      787937 :     return ret;
      91             : }
      92             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      93             : bool operator == (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
      94             :     const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
      95             : {
      96          10 :     return a.A == b.A;
      97             : }
      98             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
      99             : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator + (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
     100             :     const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
     101             : {
     102        9920 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     103        9920 :     ret.A += b.A;
     104        9920 :     return ret;
     105             : }
     106             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
     107             : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator - (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
     108             :     const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& b)
     109             : {
     110     1541250 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     111     1541250 :     ret.A -= b.A;
     112     1541250 :     return ret;
     113             : }
     114             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
     115             : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
     116             :     const Set::Scalar& b)
     117             : {
     118     1583221 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     119     1583221 :     ret.A *= b;
     120     1583221 :     return ret;
     121             : }
     122             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
     123             : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator * ( const Set::Scalar& b,
     124             :                                                     const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a)
     125             : {
     126       26574 :     return a*b;
     127             : }
     128             : 
     129             : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
     130             : Matrix4<AMREX_SPACEDIM, Sym::Diagonal> operator / (const Matrix4<AMREX_SPACEDIM, Sym::Diagonal>& a,
     131             :     const Set::Scalar& b)
     132             : {
     133     1542390 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     134     1542390 :     ret.A /= b;
     135     1542390 :     return ret;
     136             : }
     137             : 
     138             : 
     139             : }
     140             : #endif