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      2443539 :     Set::Matrix ret;
      69              : 
      70      7332427 :     for (int i = 0; i < AMREX_SPACEDIM; i++)
      71              :     {
      72     14672094 :         for (int j = 0; j < AMREX_SPACEDIM; j++)
      73              :         {
      74      9783206 :             ret(i, j) = a.A(i, j) * b(i, j);
      75              :         }
      76              :     }
      77              : 
      78      2443539 :     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      2363811 :     for (int i = 0; i < AMREX_SPACEDIM; i++)
      87      4727622 :         for (int j = 0; j < AMREX_SPACEDIM; j++)
      88      6303496 :             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      1541254 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     111      1541254 :     ret.A -= b.A;
     112      1541254 :     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      1583223 :     Matrix4<AMREX_SPACEDIM, Sym::Diagonal> ret = a;
     119      1583223 :     ret.A *= b;
     120      1583223 :     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