#include "macrodef.h"
#include "bssn_rhs.h"
#include "share_func.h"
#include "tool.h"
#include <time.h>
// 0-based i,j,k
// #define IDX_F(i,j,k,nx,ny) ((i) + (j)*(nx) + (k)*(nx)*(ny))
// ex(1)=nx, ex(2)=ny, ex(3)=nz

// 用法：a[ IDX_F(i,j,k,nx,ny) ]

#ifndef BSSN_KERNEL_FINE_TIMING
#define BSSN_KERNEL_FINE_TIMING 0
#endif

#if BSSN_KERNEL_FINE_TIMING
namespace rhs_kernel_timing
{
    enum Bucket
    {
        KB_SETUP_DERIVS = 0,
        KB_GEOM_GAMMA,
        KB_RICCI_METRIC,
        KB_CHI_LAPSE,
        KB_AIJ_TRK_GAUGE,
        KB_KO_CONSTRAINT,
        KB_COUNT
    };

    static double local_bucket_seconds[KB_COUNT];

    static const char *bucket_labels[KB_COUNT] =
    {
        "setup_derivs",
        "geom_gamma",
        "ricci_metric",
        "chi_lapse",
        "aij_trk_gauge",
        "ko_constraint"
    };

    static inline double now_seconds()
    {
        struct timespec ts;
        clock_gettime(CLOCK_MONOTONIC, &ts);
        return double(ts.tv_sec) + 1.0e-9 * double(ts.tv_nsec);
    }
}

extern "C" void f_bssn_rhs_kernel_timing_reset()
{
    for (int i = 0; i < rhs_kernel_timing::KB_COUNT; ++i)
        rhs_kernel_timing::local_bucket_seconds[i] = 0.0;
}

extern "C" int f_bssn_rhs_kernel_timing_bucket_count()
{
    return rhs_kernel_timing::KB_COUNT;
}

extern "C" const double *f_bssn_rhs_kernel_timing_local_seconds()
{
    return rhs_kernel_timing::local_bucket_seconds;
}

extern "C" const char *f_bssn_rhs_kernel_timing_label(int bucket_index)
{
    if (bucket_index < 0 || bucket_index >= rhs_kernel_timing::KB_COUNT)
        return "unknown";
    return rhs_kernel_timing::bucket_labels[bucket_index];
}

#define RHS_KERNEL_TIMER_DECL(var_name) const double var_name = rhs_kernel_timing::now_seconds()
#define RHS_KERNEL_TIMER_ADD(bucket_name, var_name) \
    rhs_kernel_timing::local_bucket_seconds[int(rhs_kernel_timing::bucket_name)] += \
        rhs_kernel_timing::now_seconds() - (var_name)
#else
extern "C" void f_bssn_rhs_kernel_timing_reset() {}
extern "C" int f_bssn_rhs_kernel_timing_bucket_count() { return 0; }
extern "C" const double *f_bssn_rhs_kernel_timing_local_seconds() { return 0; }
extern "C" const char *f_bssn_rhs_kernel_timing_label(int) { return "disabled"; }

#define RHS_KERNEL_TIMER_DECL(var_name)
#define RHS_KERNEL_TIMER_ADD(bucket_name, var_name)
#endif

// C function that calculates the right-hand side for BSSN equations
int f_compute_rhs_bssn(int *ex, double &T, 
                       double *X, double *Y, double *Z,
                       double *chi, double *trK,
                       double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
                       double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
                       double *Gamx, double *Gamy, double *Gamz,
                       double *Lap, double *betax, double *betay, double *betaz,
                       double *dtSfx, double *dtSfy, double *dtSfz,
                       double *chi_rhs, double *trK_rhs,
                       double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
                       double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
                       double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
                       double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
                       double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
                       double *rho, double *Sx, double *Sy, double *Sz,
                       double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
                       double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
                       double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
                       double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
                       double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
                       double *ham_Res, double *movx_Res, double *movy_Res, double *movz_Res,
                       double *Gmx_Res, double *Gmy_Res, double *Gmz_Res,
                       int &Symmetry, int &Lev, double &eps, int &co
                       )  // return gont
{   
    int nx = ex[0], ny = ex[1], nz=ex[2];
    int all = nx*ny*nz;
    // printf("nx=%d ny=%d nz=%d all=%d\n", nx, ny, nz, all);

    // temp variable
    double chix[all],chiy[all],chiz[all];
    double gxxx[all],gxyx[all],gxzx[all],gyyx[all],gyzx[all],gzzx[all];
    double gxxy[all],gxyy[all],gxzy[all],gyyy[all],gyzy[all],gzzy[all];
    double gxxz[all],gxyz[all],gxzz[all],gyyz[all],gyzz[all],gzzz[all];
    double Lapx[all], Lapy[all], Lapz[all];
    double betaxx[all], betaxy[all], betaxz[all];
    double betayx[all], betayy[all], betayz[all];
    double betazx[all], betazy[all], betazz[all];
    double Gamxx[all],Gamxy[all],Gamxz[all];
    double Gamyx[all],Gamyy[all],Gamyz[all];
    double Gamzx[all],Gamzy[all],Gamzz[all];
    double Kx[all], Ky[all], Kz[all], S[all];
    double f[all], fxx[all], fxy[all], fxz[all], fyy[all], fyz[all], fzz[all];
    double alpn1[all], chin1[all];
    double gupxx[all], gupxy[all], gupxz[all];
    double gupyy[all], gupyz[all], gupzz[all];
    double SSS[3] = { 1.0,  1.0,  1.0};
    double AAS[3] = {-1.0, -1.0,  1.0};
    double ASA[3] = {-1.0,  1.0, -1.0};
    double SAA[3] = { 1.0, -1.0, -1.0};
    double ASS[3] = {-1.0,  1.0,  1.0};
    double SAS[3] = { 1.0, -1.0,  1.0};
    double SSA[3] = { 1.0,  1.0, -1.0};
    double dX, dY, dZ, PI;
    const double ZEO = 0.0, ONE = 1.0, TWO = 2.0, FOUR = 4.0;
    const double EIGHT = 8.0, HALF = 0.5, THR = 3.0;
    const double SYM = 1.0, ANTI = -1.0;
    const double FF = 0.75, eta = 2.0;
    const double F1o3 = 1.0/3.0, F2o3 = 2.0/3.0, F3o2 = 1.5, F1o6 = 1.0/6.0;
    const double F16 = 16.0, F8 = 8.0;
    #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7)
    double reta[all];
    /* 使用时：reta[idx]，其中 idx = i + nx*(j + ny*k)  (Fortran列主序) */
    #endif

    #if (GAUGE == 6 || GAUGE == 7)
    int BHN = 0;
    double Porg[9] = {0.0};
    double Mass[3] = {0.0};

    extern "C" {
        #ifdef fortran1
        void getpbh(int &, double *, double *);
        #elif defined(fortran2)
        void GETPBH(int &, double *, double *);
        #else
        void getpbh_(int &, double *, double *);
        #endif
    }

    #ifdef fortran1
    getpbh(BHN, Porg, Mass);
    #elif defined(fortran2)
    GETPBH(BHN, Porg, Mass);
    #else
    getpbh_(BHN, Porg, Mass);
    #endif
    #endif
    PI = acos(-1.0); 
    dX = X[1] - X[0];
    dY = Y[1] - Y[0];
    dZ = Z[1] - Z[0];

        RHS_KERNEL_TIMER_DECL(timer_setup_derivs);
        // 1ms //
        for(int i=0;i<all;i+=1){
            alpn1[i] = Lap[i] + 1.0;
            chin1[i] = chi[i] + 1.0;
        }
        // 9ms //
        fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev);
        fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev);
        fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev);
        fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
        fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
        fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev);
        fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev);
        fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
        fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev);
        fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
        fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
        fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);

        // 3ms //
        for(int i=0;i<all;i+=1){
            const double divb = betaxx[i] + betayy[i] + betazz[i];
            chi_rhs[i] = F2o3 * chin1[i] * (alpn1[i] * trK[i] - divb); 
            gxx_rhs[i] = -TWO * alpn1[i] * Axx[i] - F2o3 * (dxx[i] + ONE) * divb + 
                    TWO * ((dxx[i] + ONE) * betaxx[i] + gxy[i] * betayx[i] + gxz[i] * betazx[i]);
            gyy_rhs[i] = -TWO * alpn1[i] * Ayy[i] - F2o3 * (dyy[i] + ONE) * divb + 
                    TWO * (gxy[i] * betaxy[i] + (dyy[i] + ONE) * betayy[i] + gyz[i] * betazy[i]);
            gzz_rhs[i] = -TWO * alpn1[i] * Azz[i] - F2o3 * (dzz[i] + ONE) * divb + 
                        TWO * (gxz[i] * betaxz[i] + gyz[i] * betayz[i] + (dzz[i] + ONE) * betazz[i]);
            gxy_rhs[i] = -TWO * alpn1[i] * Axy[i] + F1o3 * gxy[i] * divb + 
                        (dxx[i] + ONE) * betaxy[i] + gxz[i] * betazy[i] + (dyy[i] + ONE) * betayx[i] 
                        + gyz[i] * betazx[i] - gxy[i] * betazz[i];        
            gyz_rhs[i] = -TWO * alpn1[i] * Ayz[i] + F1o3 * gyz[i] * divb + 
                        gxy[i] * betaxz[i] + (dyy[i] + ONE) * betayz[i] + gxz[i] * betaxy[i] 
                        + (dzz[i] + ONE) * betazy[i] - gyz[i] * betaxx[i];
            gxz_rhs[i] = -TWO * alpn1[i] * Axz[i] + F1o3 * gxz[i] * divb + 
                        (dxx[i] + ONE) * betaxz[i] + gxy[i] * betayz[i] + gyz[i] * betayx[i] 
                        + (dzz[i] + ONE) * betazx[i] - gxz[i] * betayy[i];                  
        }
        RHS_KERNEL_TIMER_ADD(KB_SETUP_DERIVS, timer_setup_derivs);
        RHS_KERNEL_TIMER_DECL(timer_geom_gamma);
        // Fused: inverse metric + Gamma constraint + Christoffel (3 loops -> 1)
        for(int i=0;i<all;i+=1){
            double det = (dxx[i] + ONE) * (dyy[i] + ONE) * (dzz[i] + ONE) + gxy[i] * gyz[i] * gxz[i] + gxz[i] * gxy[i] * gyz[i] -
                        gxz[i] * (dyy[i] + ONE) * gxz[i] - gxy[i] * gxy[i] * (dzz[i] + ONE) - (dxx[i] + ONE) * gyz[i] * gyz[i];
            double lg_xx = ((dyy[i] + ONE) * (dzz[i] + ONE) - gyz[i] * gyz[i]) / det;
            double lg_xy = -(gxy[i] * (dzz[i] + ONE) - gyz[i] * gxz[i]) / det;
            double lg_xz = (gxy[i] * gyz[i] - (dyy[i] + ONE) * gxz[i]) / det;
            double lg_yy = ((dxx[i] + ONE) * (dzz[i] + ONE) - gxz[i] * gxz[i]) / det;
            double lg_yz = -((dxx[i] + ONE) * gyz[i] - gxy[i] * gxz[i]) / det;
            double lg_zz = ((dxx[i] + ONE) * (dyy[i] + ONE) - gxy[i] * gxy[i]) / det;
            gupxx[i] = lg_xx; gupxy[i] = lg_xy; gupxz[i] = lg_xz;
            gupyy[i] = lg_yy; gupyz[i] = lg_yz; gupzz[i] = lg_zz;

            if(co==0){
                Gmx_Res[i] = Gamx[i] - (
                    lg_xx * (lg_xx*gxxx[i] + lg_xy*gxyx[i] + lg_xz*gxzx[i]) +
                    lg_xy * (lg_xx*gxyx[i] + lg_xy*gyyx[i] + lg_xz*gyzx[i]) +
                    lg_xz * (lg_xx*gxzx[i] + lg_xy*gyzx[i] + lg_xz*gzzx[i]) +
                    lg_xx * (lg_xy*gxxy[i] + lg_yy*gxyy[i] + lg_yz*gxzy[i]) +
                    lg_xy * (lg_xy*gxyy[i] + lg_yy*gyyy[i] + lg_yz*gyzy[i]) +
                    lg_xz * (lg_xy*gxzy[i] + lg_yy*gyzy[i] + lg_yz*gzzy[i]) +
                    lg_xx * (lg_xz*gxxz[i] + lg_yz*gxyz[i] + lg_zz*gxzz[i]) +
                    lg_xy * (lg_xz*gxyz[i] + lg_yz*gyyz[i] + lg_zz*gyzz[i]) +
                    lg_xz * (lg_xz*gxzz[i] + lg_yz*gyzz[i] + lg_zz*gzzz[i])
                );
                Gmy_Res[i] = Gamy[i] - (
                    lg_xx * (lg_xy*gxxx[i] + lg_yy*gxyx[i] + lg_yz*gxzx[i]) +
                    lg_xy * (lg_xy*gxyx[i] + lg_yy*gyyx[i] + lg_yz*gyzx[i]) +
                    lg_xz * (lg_xy*gxzx[i] + lg_yy*gyzx[i] + lg_yz*gzzx[i]) +
                    lg_xy * (lg_xy*gxxy[i] + lg_yy*gxyy[i] + lg_yz*gxzy[i]) +
                    lg_yy * (lg_xy*gxyy[i] + lg_yy*gyyy[i] + lg_yz*gyzy[i]) +
                    lg_yz * (lg_xy*gxzy[i] + lg_yy*gyzy[i] + lg_yz*gzzy[i]) +
                    lg_xy * (lg_xz*gxxz[i] + lg_yz*gxyz[i] + lg_zz*gxzz[i]) +
                    lg_yy * (lg_xz*gxyz[i] + lg_yz*gyyz[i] + lg_zz*gyzz[i]) +
                    lg_yz * (lg_xz*gxzz[i] + lg_yz*gyzz[i] + lg_zz*gzzz[i])
                );
                Gmz_Res[i] = Gamz[i] - (
                    lg_xx * (lg_xz*gxxx[i] + lg_yz*gxyx[i] + lg_zz*gxzx[i]) +
                    lg_xy * (lg_xz*gxyx[i] + lg_yz*gyyx[i] + lg_zz*gyzx[i]) +
                    lg_xz * (lg_xz*gxzx[i] + lg_yz*gyzx[i] + lg_zz*gzzx[i]) +
                    lg_xy * (lg_xz*gxxy[i] + lg_yz*gxyy[i] + lg_zz*gxzy[i]) +
                    lg_yy * (lg_xz*gxyy[i] + lg_yz*gyyy[i] + lg_zz*gyzy[i]) +
                    lg_yz * (lg_xz*gxzy[i] + lg_yz*gyzy[i] + lg_zz*gzzy[i]) +
                    lg_xz * (lg_xz*gxxz[i] + lg_yz*gxyz[i] + lg_zz*gxzz[i]) +
                    lg_yz * (lg_xz*gxyz[i] + lg_yz*gyyz[i] + lg_zz*gyzz[i]) +
                    lg_zz * (lg_xz*gxzz[i] + lg_yz*gyzz[i] + lg_zz*gzzz[i])
                );
            }

            Gamxxx[i] = HALF * ( lg_xx*gxxx[i]
                            + lg_xy*(TWO*gxyx[i] - gxxy[i])
                            + lg_xz*(TWO*gxzx[i] - gxxz[i]) );
            Gamyxx[i] = HALF * ( lg_xy*gxxx[i]
                            + lg_yy*(TWO*gxyx[i] - gxxy[i])
                            + lg_yz*(TWO*gxzx[i] - gxxz[i]) );
            Gamzxx[i] = HALF * ( lg_xz*gxxx[i]
                            + lg_yz*(TWO*gxyx[i] - gxxy[i])
                            + lg_zz*(TWO*gxzx[i] - gxxz[i]) );
            Gamxyy[i] = HALF * ( lg_xx*(TWO*gxyy[i] - gyyx[i])
                            + lg_xy*gyyy[i]
                            + lg_xz*(TWO*gyzy[i] - gyyz[i]) );
            Gamyyy[i] = HALF * ( lg_xy*(TWO*gxyy[i] - gyyx[i])
                            + lg_yy*gyyy[i]
                            + lg_yz*(TWO*gyzy[i] - gyyz[i]) );
            Gamzyy[i] = HALF * ( lg_xz*(TWO*gxyy[i] - gyyx[i])
                            + lg_yz*gyyy[i]
                            + lg_zz*(TWO*gyzy[i] - gyyz[i]) );
            Gamxzz[i] = HALF * ( lg_xx*(TWO*gxzz[i] - gzzx[i])
                            + lg_xy*(TWO*gyzz[i] - gzzy[i])
                            + lg_xz*gzzz[i] );
            Gamyzz[i] = HALF * ( lg_xy*(TWO*gxzz[i] - gzzx[i])
                            + lg_yy*(TWO*gyzz[i] - gzzy[i])
                            + lg_yz*gzzz[i] );
            Gamzzz[i] = HALF * ( lg_xz*(TWO*gxzz[i] - gzzx[i])
                            + lg_yz*(TWO*gyzz[i] - gzzy[i])
                            + lg_zz*gzzz[i] );
            Gamxxy[i] = HALF * ( lg_xx*gxxy[i]
                            + lg_xy*gyyx[i]
                            + lg_xz*(gxzy[i] + gyzx[i] - gxyz[i]) );
            Gamyxy[i] = HALF * ( lg_xy*gxxy[i]
                            + lg_yy*gyyx[i]
                            + lg_yz*(gxzy[i] + gyzx[i] - gxyz[i]) );
            Gamzxy[i] = HALF * ( lg_xz*gxxy[i]
                            + lg_yz*gyyx[i]
                            + lg_zz*(gxzy[i] + gyzx[i] - gxyz[i]) );
            Gamxxz[i] = HALF * ( lg_xx*gxxz[i]
                            + lg_xy*(gxyz[i] + gyzx[i] - gxzy[i])
                            + lg_xz*gzzx[i] );
            Gamyxz[i] = HALF * ( lg_xy*gxxz[i]
                            + lg_yy*(gxyz[i] + gyzx[i] - gxzy[i])
                            + lg_yz*gzzx[i] );
            Gamzxz[i] = HALF * ( lg_xz*gxxz[i]
                            + lg_yz*(gxyz[i] + gyzx[i] - gxzy[i])
                            + lg_zz*gzzx[i] );
            Gamxyz[i] = HALF * ( lg_xx*(gxyz[i] + gxzy[i] - gyzx[i])
                            + lg_xy*gyyz[i]
                            + lg_xz*gzzy[i] );
            Gamyyz[i] = HALF * ( lg_xy*(gxyz[i] + gxzy[i] - gyzx[i])
                            + lg_yy*gyyz[i]
                            + lg_yz*gzzy[i] );
            Gamzyz[i] = HALF * ( lg_xz*(gxyz[i] + gxzy[i] - gyzx[i])
                            + lg_yz*gyyz[i]
                            + lg_zz*gzzy[i] );
        }
        // Fused: A^{ij} raise-index + Gamma_rhs part 1 (2 loops -> 1)
        for (int i=0;i<all;i+=1) {
            double axx = gupxx[i]*gupxx[i]*Axx[i]
                + gupxy[i]*gupxy[i]*Ayy[i]
                + gupxz[i]*gupxz[i]*Azz[i]
                + TWO * ( gupxx[i]*gupxy[i]*Axy[i]
                        + gupxx[i]*gupxz[i]*Axz[i]
                        + gupxy[i]*gupxz[i]*Ayz[i] );
            double ayy = gupxy[i]*gupxy[i]*Axx[i]
                + gupyy[i]*gupyy[i]*Ayy[i]
                + gupyz[i]*gupyz[i]*Azz[i]
                + TWO * ( gupxy[i]*gupyy[i]*Axy[i]
                        + gupxy[i]*gupyz[i]*Axz[i]
                        + gupyy[i]*gupyz[i]*Ayz[i] );
            double azz = gupxz[i]*gupxz[i]*Axx[i]
                + gupyz[i]*gupyz[i]*Ayy[i]
                + gupzz[i]*gupzz[i]*Azz[i]
                + TWO * ( gupxz[i]*gupyz[i]*Axy[i]
                        + gupxz[i]*gupzz[i]*Axz[i]
                        + gupyz[i]*gupzz[i]*Ayz[i] );
            double axy = gupxx[i]*gupxy[i]*Axx[i]
                + gupxy[i]*gupyy[i]*Ayy[i]
                + gupxz[i]*gupyz[i]*Azz[i]
                + ( gupxx[i]*gupyy[i] + gupxy[i]*gupxy[i] ) * Axy[i]
                + ( gupxx[i]*gupyz[i] + gupxz[i]*gupxy[i] ) * Axz[i]
                + ( gupxy[i]*gupyz[i] + gupxz[i]*gupyy[i] ) * Ayz[i];
            double axz = gupxx[i]*gupxz[i]*Axx[i]
                + gupxy[i]*gupyz[i]*Ayy[i]
                + gupxz[i]*gupzz[i]*Azz[i]
                + ( gupxx[i]*gupyz[i] + gupxy[i]*gupxz[i] ) * Axy[i]
                + ( gupxx[i]*gupzz[i] + gupxz[i]*gupxz[i] ) * Axz[i]
                + ( gupxy[i]*gupzz[i] + gupxz[i]*gupyz[i] ) * Ayz[i];
            double ayz = gupxy[i]*gupxz[i]*Axx[i]
                + gupyy[i]*gupyz[i]*Ayy[i]
                + gupyz[i]*gupzz[i]*Azz[i]
                + ( gupxy[i]*gupyz[i] + gupyy[i]*gupxz[i] ) * Axy[i]
                + ( gupxy[i]*gupzz[i] + gupyz[i]*gupxz[i] ) * Axz[i]
                + ( gupyy[i]*gupzz[i] + gupyz[i]*gupyz[i] ) * Ayz[i];
            Gamx_rhs[i] = - TWO * ( Lapx[i]*axx + Lapy[i]*axy + Lapz[i]*axz ) +
                TWO * alpn1[i] * (
                -F3o2/chin1[i] * ( chix[i]*axx + chiy[i]*axy + chiz[i]*axz ) -
                    gupxx[i] * ( F2o3*Kx[i] + EIGHT*PI*Sx[i] ) -
                    gupxy[i] * ( F2o3*Ky[i] + EIGHT*PI*Sy[i] ) -
                    gupxz[i] * ( F2o3*Kz[i] + EIGHT*PI*Sz[i] ) +
                                Gamxxx[i]*axx + Gamxyy[i]*ayy + Gamxzz[i]*azz +
                        TWO * ( Gamxxy[i]*axy + Gamxxz[i]*axz + Gamxyz[i]*ayz ) );

            Gamy_rhs[i] = -TWO * ( Lapx[i]*axy + Lapy[i]*ayy + Lapz[i]*ayz )
                        +  TWO * alpn1[i] * (
                            -F3o2/chin1[i] * ( chix[i]*axy + chiy[i]*ayy + chiz[i]*ayz )
                            - gupxy[i] * ( F2o3*Kx[i] + EIGHT*PI*Sx[i] )
                            - gupyy[i] * ( F2o3*Ky[i] + EIGHT*PI*Sy[i] )
                            - gupyz[i] * ( F2o3*Kz[i] + EIGHT*PI*Sz[i] )
                            + Gamyxx[i]*axx + Gamyyy[i]*ayy + Gamyzz[i]*azz
                            + TWO * ( Gamyxy[i]*axy + Gamyxz[i]*axz + Gamyyz[i]*ayz )
                        );

            Gamz_rhs[i] = -TWO * ( Lapx[i]*axz + Lapy[i]*ayz + Lapz[i]*azz )
                        +  TWO * alpn1[i] * (
                            -F3o2/chin1[i] * ( chix[i]*axz + chiy[i]*ayz + chiz[i]*azz )
                            - gupxz[i] * ( F2o3*Kx[i] + EIGHT*PI*Sx[i] )
                            - gupyz[i] * ( F2o3*Ky[i] + EIGHT*PI*Sy[i] )
                            - gupzz[i] * ( F2o3*Kz[i] + EIGHT*PI*Sz[i] )
                            + Gamzxx[i]*axx + Gamzyy[i]*ayy + Gamzzz[i]*azz
                            + TWO * ( Gamzxy[i]*axy + Gamzxz[i]*axz + Gamzyz[i]*ayz )
                        );
        }
        RHS_KERNEL_TIMER_ADD(KB_GEOM_GAMMA, timer_geom_gamma);
        RHS_KERNEL_TIMER_DECL(timer_ricci_metric);
        // 22.3ms //
        fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,
        X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev);
        fdderivs(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,
        X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev);
        fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,
        X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev);
        fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev);
        fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev);
        fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev);

        // Fused: fxx/Gamxa + Gamma_rhs part 2 (2 loops -> 1)
        for(int i=0;i<all;i+=1){
            const double divb = betaxx[i] + betayy[i] + betazz[i];
            double lfxx = gxxx[i] + gxyy[i] + gxzz[i];
            double lfxy = gxyx[i] + gyyy[i] + gyzz[i];
            double lfxz = gxzx[i] + gyzy[i] + gzzz[i];

            double gxa = gupxx[i]*Gamxxx[i] + gupyy[i]*Gamxyy[i] + gupzz[i]*Gamxzz[i]
                    + TWO * ( gupxy[i]*Gamxxy[i] + gupxz[i]*Gamxxz[i] + gupyz[i]*Gamxyz[i] );
            double gya = gupxx[i]*Gamyxx[i] + gupyy[i]*Gamyyy[i] + gupzz[i]*Gamyzz[i]
                    + TWO * ( gupxy[i]*Gamyxy[i] + gupxz[i]*Gamyxz[i] + gupyz[i]*Gamyyz[i] );
            double gza = gupxx[i]*Gamzxx[i] + gupyy[i]*Gamzyy[i] + gupzz[i]*Gamzzz[i]
                    + TWO * ( gupxy[i]*Gamzxy[i] + gupxz[i]*Gamzxz[i] + gupyz[i]*Gamzyz[i] );
            Gamx_rhs[i] = Gamx_rhs[i]
                        + F2o3 * gxa * divb
                        - gxa * betaxx[i] - gya * betaxy[i] - gza * betaxz[i]
                        + F1o3 * ( gupxx[i] * lfxx + gupxy[i] * lfxy + gupxz[i] * lfxz )
                        + gupxx[i] * gxxx[i] + gupyy[i] * gyyx[i] + gupzz[i] * gzzx[i]
                        + TWO * ( gupxy[i] * gxyx[i] + gupxz[i] * gxzx[i] + gupyz[i] * gyzx[i] );

            Gamy_rhs[i] = Gamy_rhs[i]
                        + F2o3 * gya * divb
                        - gxa * betayx[i] - gya * betayy[i] - gza * betayz[i]
                        + F1o3 * ( gupxy[i] * lfxx + gupyy[i] * lfxy + gupyz[i] * lfxz )
                        + gupxx[i] * gxxy[i] + gupyy[i] * gyyy[i] + gupzz[i] * gzzy[i]
                        + TWO * ( gupxy[i] * gxyy[i] + gupxz[i] * gxzy[i] + gupyz[i] * gyzy[i] );

            Gamz_rhs[i] = Gamz_rhs[i]
                        + F2o3 * gza * divb
                        - gxa * betazx[i] - gya * betazy[i] - gza * betazz[i]
                        + F1o3 * ( gupxz[i] * lfxx + gupyz[i] * lfxy + gupzz[i] * lfxz )
                        + gupxx[i] * gxxz[i] + gupyy[i] * gyyz[i] + gupzz[i] * gzzz[i]
                        + TWO * ( gupxy[i] * gxyz[i] + gupxz[i] * gxzz[i] + gupyz[i] * gyzz[i] );
        }
        // 4.4ms //
        for (int i=0;i<all;i+=1) {
            gxxx[i] = (dxx[i] + ONE)*Gamxxx[i] + gxy[i]*Gamyxx[i] + gxz[i]*Gamzxx[i];
            gxyx[i] = (dxx[i] + ONE)*Gamxxy[i] + gxy[i]*Gamyxy[i] + gxz[i]*Gamzxy[i];
            gxzx[i] = (dxx[i] + ONE)*Gamxxz[i] + gxy[i]*Gamyxz[i] + gxz[i]*Gamzxz[i];
            gyyx[i] = (dxx[i] + ONE)*Gamxyy[i] + gxy[i]*Gamyyy[i] + gxz[i]*Gamzyy[i];
            gyzx[i] = (dxx[i] + ONE)*Gamxyz[i] + gxy[i]*Gamyyz[i] + gxz[i]*Gamzyz[i];
            gzzx[i] = (dxx[i] + ONE)*Gamxzz[i] + gxy[i]*Gamyzz[i] + gxz[i]*Gamzzz[i];

            gxxy[i] = gxy[i]*Gamxxx[i] + (dyy[i] + ONE)*Gamyxx[i] + gyz[i]*Gamzxx[i];
            gxyy[i] = gxy[i]*Gamxxy[i] + (dyy[i] + ONE)*Gamyxy[i] + gyz[i]*Gamzxy[i];
            gxzy[i] = gxy[i]*Gamxxz[i] + (dyy[i] + ONE)*Gamyxz[i] + gyz[i]*Gamzxz[i];
            gyyy[i] = gxy[i]*Gamxyy[i] + (dyy[i] + ONE)*Gamyyy[i] + gyz[i]*Gamzyy[i];
            gyzy[i] = gxy[i]*Gamxyz[i] + (dyy[i] + ONE)*Gamyyz[i] + gyz[i]*Gamzyz[i];
            gzzy[i] = gxy[i]*Gamxzz[i] + (dyy[i] + ONE)*Gamyzz[i] + gyz[i]*Gamzzz[i];

            gxxz[i] = gxz[i]*Gamxxx[i] + gyz[i]*Gamyxx[i] + (dzz[i] + ONE)*Gamzxx[i];
            gxyz[i] = gxz[i]*Gamxxy[i] + gyz[i]*Gamyxy[i] + (dzz[i] + ONE)*Gamzxy[i];
            gxzz[i] = gxz[i]*Gamxxz[i] + gyz[i]*Gamyxz[i] + (dzz[i] + ONE)*Gamzxz[i];
            gyyz[i] = gxz[i]*Gamxyy[i] + gyz[i]*Gamyyy[i] + (dzz[i] + ONE)*Gamzyy[i];
            gyzz[i] = gxz[i]*Gamxyz[i] + gyz[i]*Gamyyz[i] + (dzz[i] + ONE)*Gamzyz[i];
            gzzz[i] = gxz[i]*Gamxzz[i] + gyz[i]*Gamyzz[i] + (dzz[i] + ONE)*Gamzzz[i];
        }
        // 22.2ms //
            fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);

        // 0.7ms //
        for (int i = 0; i < all; i += 1) {
            Rxx[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }
        // 23ms //
            fdderivs(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
        for (int i = 0; i < all; i += 1) {
            Ryy[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }

        // 23ms //
            fdderivs(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev);
        for (int i = 0; i < all; i += 1) {
            Rzz[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }

        // 23ms //
            fdderivs(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,Symmetry,Lev);
        for (int i = 0; i < all; i += 1) {
            Rxy[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }

        // 23ms //
            fdderivs(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,Symmetry,Lev);
        for (int i = 0; i < all; i += 1) {
            Rxz[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }

        // 23ms //
            fdderivs(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,Symmetry,Lev);
        for (int i = 0; i < all; i += 1) {
            Ryz[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    + (gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i]) * TWO;
        }

        // 14ms //
        /* 假设 all = ex1*ex2*ex3，所有量都是 length=all 的 double 数组（已按同一扁平化规则排布） */
        for (int i = 0; i < all; i += 1) {
            const double gxa = gupxx[i]*Gamxxx[i] + gupyy[i]*Gamxyy[i] + gupzz[i]*Gamxzz[i]
                            + TWO * ( gupxy[i]*Gamxxy[i] + gupxz[i]*Gamxxz[i] + gupyz[i]*Gamxyz[i] );
            const double gya = gupxx[i]*Gamyxx[i] + gupyy[i]*Gamyyy[i] + gupzz[i]*Gamyzz[i]
                            + TWO * ( gupxy[i]*Gamyxy[i] + gupxz[i]*Gamyxz[i] + gupyz[i]*Gamyyz[i] );
            const double gza = gupxx[i]*Gamzxx[i] + gupyy[i]*Gamzyy[i] + gupzz[i]*Gamzzz[i]
                            + TWO * ( gupxy[i]*Gamzxy[i] + gupxz[i]*Gamzxz[i] + gupyz[i]*Gamzyz[i] );

            Rxx[i] =
                -HALF * Rxx[i]
                + (dxx[i] + ONE) * Gamxx[i] + gxy[i] * Gamyx[i] + gxz[i] * Gamzx[i]
                + gxa * gxxx[i] + gya * gxyx[i] + gza * gxzx[i]
                + gupxx[i] * (
                    TWO * (Gamxxx[i] * gxxx[i] + Gamyxx[i] * gxyx[i] + Gamzxx[i] * gxzx[i]) +
                        (Gamxxx[i] * gxxx[i] + Gamyxx[i] * gxxy[i] + Gamzxx[i] * gxxz[i])
                )
                + gupxy[i] * (
                    TWO * (Gamxxx[i] * gxyx[i] + Gamyxx[i] * gyyx[i] + Gamzxx[i] * gyzx[i] +
                        Gamxxy[i] * gxxx[i] + Gamyxy[i] * gxyx[i] + Gamzxy[i] * gxzx[i]) +
                        (Gamxxy[i] * gxxx[i] + Gamyxy[i] * gxxy[i] + Gamzxy[i] * gxxz[i]) +
                        (Gamxxx[i] * gxyx[i] + Gamyxx[i] * gxyy[i] + Gamzxx[i] * gxyz[i])
                )
                + gupxz[i] * (
                    TWO * (Gamxxx[i] * gxzx[i] + Gamyxx[i] * gyzx[i] + Gamzxx[i] * gzzx[i] +
                        Gamxxz[i] * gxxx[i] + Gamyxz[i] * gxyx[i] + Gamzxz[i] * gxzx[i]) +
                        (Gamxxz[i] * gxxx[i] + Gamyxz[i] * gxxy[i] + Gamzxz[i] * gxxz[i]) +
                        (Gamxxx[i] * gxzx[i] + Gamyxx[i] * gxzy[i] + Gamzxx[i] * gxzz[i])
                )
                + gupyy[i] * (
                    TWO * (Gamxxy[i] * gxyx[i] + Gamyxy[i] * gyyx[i] + Gamzxy[i] * gyzx[i]) +
                        (Gamxxy[i] * gxyx[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxyz[i])
                )
                + gupyz[i] * (
                    TWO * (Gamxxy[i] * gxzx[i] + Gamyxy[i] * gyzx[i] + Gamzxy[i] * gzzx[i] +
                        Gamxxz[i] * gxyx[i] + Gamyxz[i] * gyyx[i] + Gamzxz[i] * gyzx[i]) +
                        (Gamxxz[i] * gxyx[i] + Gamyxz[i] * gxyy[i] + Gamzxz[i] * gxyz[i]) +
                        (Gamxxy[i] * gxzx[i] + Gamyxy[i] * gxzy[i] + Gamzxy[i] * gxzz[i])
                )
                + gupzz[i] * (
                    TWO * (Gamxxz[i] * gxzx[i] + Gamyxz[i] * gyzx[i] + Gamzxz[i] * gzzx[i]) +
                        (Gamxxz[i] * gxzx[i] + Gamyxz[i] * gxzy[i] + Gamzxz[i] * gxzz[i])
                );

            Ryy[i] =
                -HALF * Ryy[i]
                + gxy[i] * Gamxy[i] + (dyy[i] + ONE) * Gamyy[i] + gyz[i] * Gamzy[i]
                + gxa * gxyy[i] + gya * gyyy[i] + gza * gyzy[i]
                + gupxx[i] * (
                    TWO * (Gamxxy[i] * gxxy[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxzy[i]) +
                        (Gamxxy[i] * gxyx[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxyz[i])
                )
                + gupxy[i] * (
                    TWO * (Gamxxy[i] * gxyy[i] + Gamyxy[i] * gyyy[i] + Gamzxy[i] * gyzy[i] +
                        Gamxyy[i] * gxxy[i] + Gamyyy[i] * gxyy[i] + Gamzyy[i] * gxzy[i]) +
                        (Gamxyy[i] * gxyx[i] + Gamyyy[i] * gxyy[i] + Gamzyy[i] * gxyz[i]) +
                        (Gamxxy[i] * gyyx[i] + Gamyxy[i] * gyyy[i] + Gamzxy[i] * gyyz[i])
                )
                + gupxz[i] * (
                    TWO * (Gamxxy[i] * gxzy[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gzzy[i] +
                        Gamxyz[i] * gxxy[i] + Gamyyz[i] * gxyy[i] + Gamzyz[i] * gxzy[i]) +
                        (Gamxyz[i] * gxyx[i] + Gamyyz[i] * gxyy[i] + Gamzyz[i] * gxyz[i]) +
                        (Gamxxy[i] * gyzx[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gyzz[i])
                )
                + gupyy[i] * (
                    TWO * (Gamxyy[i] * gxyy[i] + Gamyyy[i] * gyyy[i] + Gamzyy[i] * gyzy[i]) +
                        (Gamxyy[i] * gyyx[i] + Gamyyy[i] * gyyy[i] + Gamzyy[i] * gyyz[i])
                )
                + gupyz[i] * (
                    TWO * (Gamxyy[i] * gxzy[i] + Gamyyy[i] * gyzy[i] + Gamzyy[i] * gzzy[i] +
                        Gamxyz[i] * gxyy[i] + Gamyyz[i] * gyyy[i] + Gamzyz[i] * gyzy[i]) +
                        (Gamxyz[i] * gyyx[i] + Gamyyz[i] * gyyy[i] + Gamzyz[i] * gyyz[i]) +
                        (Gamxyy[i] * gyzx[i] + Gamyyy[i] * gyzy[i] + Gamzyy[i] * gyzz[i])
                )
                + gupzz[i] * (
                    TWO * (Gamxyz[i] * gxzy[i] + Gamyyz[i] * gyzy[i] + Gamzyz[i] * gzzy[i]) +
                        (Gamxyz[i] * gyzx[i] + Gamyyz[i] * gyzy[i] + Gamzyz[i] * gyzz[i])
                );

            Rzz[i] =
                -HALF * Rzz[i]
                + gxz[i] * Gamxz[i] + gyz[i] * Gamyz[i] + (dzz[i] + ONE) * Gamzz[i]
                + gxa * gxzz[i] + gya * gyzz[i] + gza * gzzz[i]
                + gupxx[i] * (
                    TWO * (Gamxxz[i] * gxxz[i] + Gamyxz[i] * gxyz[i] + Gamzxz[i] * gxzz[i]) +
                        (Gamxxz[i] * gxzx[i] + Gamyxz[i] * gxzy[i] + Gamzxz[i] * gxzz[i])
                )
                + gupxy[i] * (
                    TWO * (Gamxxz[i] * gxyz[i] + Gamyxz[i] * gyyz[i] + Gamzxz[i] * gyzz[i] +
                        Gamxyz[i] * gxxz[i] + Gamyyz[i] * gxyz[i] + Gamzyz[i] * gxzz[i]) +
                        (Gamxyz[i] * gxzx[i] + Gamyyz[i] * gxzy[i] + Gamzyz[i] * gxzz[i]) +
                        (Gamxxz[i] * gyzx[i] + Gamyxz[i] * gyzy[i] + Gamzxz[i] * gyzz[i])
                )
                + gupxz[i] * (
                    TWO * (Gamxxz[i] * gxzz[i] + Gamyxz[i] * gyzz[i] + Gamzxz[i] * gzzz[i] +
                        Gamxzz[i] * gxxz[i] + Gamyzz[i] * gxyz[i] + Gamzzz[i] * gxzz[i]) +
                        (Gamxzz[i] * gxzx[i] + Gamyzz[i] * gxzy[i] + Gamzzz[i] * gxzz[i]) +
                        (Gamxxz[i] * gzzx[i] + Gamyxz[i] * gzzy[i] + Gamzxz[i] * gzzz[i])
                )
                + gupyy[i] * (
                    TWO * (Gamxyz[i] * gxyz[i] + Gamyyz[i] * gyyz[i] + Gamzyz[i] * gyzz[i]) +
                        (Gamxyz[i] * gyzx[i] + Gamyyz[i] * gyzy[i] + Gamzyz[i] * gyzz[i])
                )
                + gupyz[i] * (
                    TWO * (Gamxyz[i] * gxzz[i] + Gamyyz[i] * gyzz[i] + Gamzyz[i] * gzzz[i] +
                        Gamxzz[i] * gxyz[i] + Gamyzz[i] * gyyz[i] + Gamzzz[i] * gyzz[i]) +
                        (Gamxzz[i] * gyzx[i] + Gamyzz[i] * gyzy[i] + Gamzzz[i] * gyzz[i]) +
                        (Gamxyz[i] * gzzx[i] + Gamyyz[i] * gzzy[i] + Gamzyz[i] * gzzz[i])
                )
                + gupzz[i] * (
                    TWO * (Gamxzz[i] * gxzz[i] + Gamyzz[i] * gyzz[i] + Gamzzz[i] * gzzz[i]) +
                        (Gamxzz[i] * gzzx[i] + Gamyzz[i] * gzzy[i] + Gamzzz[i] * gzzz[i])
                );

            Rxy[i] =
                HALF * (
                    -Rxy[i]
                    + (dxx[i] + ONE) * Gamxy[i] + gxy[i] * Gamyy[i] + gxz[i] * Gamzy[i]
                    + gxy[i] * Gamxx[i] + (dyy[i] + ONE) * Gamyx[i] + gyz[i] * Gamzx[i]
                    + gxa * gxyx[i] + gya * gyyx[i] + gza * gyzx[i]
                    + gxa * gxxy[i] + gya * gxyy[i] + gza * gxzy[i]
                )
                + gupxx[i] * (
                    Gamxxx[i] * gxxy[i] + Gamyxx[i] * gxyy[i] + Gamzxx[i] * gxzy[i]
                    + Gamxxy[i] * gxxx[i] + Gamyxy[i] * gxyx[i] + Gamzxy[i] * gxzx[i]
                    + Gamxxx[i] * gxyx[i] + Gamyxx[i] * gxyy[i] + Gamzxx[i] * gxyz[i]
                )
                + gupxy[i] * (
                    Gamxxx[i] * gxyy[i] + Gamyxx[i] * gyyy[i] + Gamzxx[i] * gyzy[i]
                    + Gamxxy[i] * gxyx[i] + Gamyxy[i] * gyyx[i] + Gamzxy[i] * gyzx[i]
                    + Gamxxy[i] * gxyx[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxyz[i]
                    + Gamxxy[i] * gxxy[i] + Gamyxy[i] * gxyy[i] + Gamzxy[i] * gxzy[i]
                    + Gamxyy[i] * gxxx[i] + Gamyyy[i] * gxyx[i] + Gamzyy[i] * gxzx[i]
                    + Gamxxx[i] * gyyx[i] + Gamyxx[i] * gyyy[i] + Gamzxx[i] * gyyz[i]
                )
                + gupxz[i] * (
                    Gamxxx[i] * gxzy[i] + Gamyxx[i] * gyzy[i] + Gamzxx[i] * gzzy[i]
                    + Gamxxy[i] * gxzx[i] + Gamyxy[i] * gyzx[i] + Gamzxy[i] * gzzx[i]
                    + Gamxxz[i] * gxyx[i] + Gamyxz[i] * gxyy[i] + Gamzxz[i] * gxyz[i]
                    + Gamxxz[i] * gxxy[i] + Gamyxz[i] * gxyy[i] + Gamzxz[i] * gxzy[i]
                    + Gamxyz[i] * gxxx[i] + Gamyyz[i] * gxyx[i] + Gamzyz[i] * gxzx[i]
                    + Gamxxx[i] * gyzx[i] + Gamyxx[i] * gyzy[i] + Gamzxx[i] * gyzz[i]
                )
                + gupyy[i] * (
                    Gamxxy[i] * gxyy[i] + Gamyxy[i] * gyyy[i] + Gamzxy[i] * gyzy[i]
                    + Gamxyy[i] * gxyx[i] + Gamyyy[i] * gyyx[i] + Gamzyy[i] * gyzx[i]
                    + Gamxxy[i] * gyyx[i] + Gamyxy[i] * gyyy[i] + Gamzxy[i] * gyyz[i]
                )
                + gupyz[i] * (
                    Gamxxy[i] * gxzy[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gzzy[i]
                    + Gamxyy[i] * gxzx[i] + Gamyyy[i] * gyzx[i] + Gamzyy[i] * gzzx[i]
                    + Gamxxz[i] * gyyx[i] + Gamyxz[i] * gyyy[i] + Gamzxz[i] * gyyz[i]
                    + Gamxxz[i] * gxyy[i] + Gamyxz[i] * gyyy[i] + Gamzxz[i] * gyzy[i]
                    + Gamxyz[i] * gxyx[i] + Gamyyz[i] * gyyx[i] + Gamzyz[i] * gyzx[i]
                    + Gamxxy[i] * gyzx[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gyzz[i]
                )
                + gupzz[i] * (
                    Gamxxz[i] * gxzy[i] + Gamyxz[i] * gyzy[i] + Gamzxz[i] * gzzy[i]
                    + Gamxyz[i] * gxzx[i] + Gamyyz[i] * gyzx[i] + Gamzyz[i] * gzzx[i]
                    + Gamxxz[i] * gyzx[i] + Gamyxz[i] * gyzy[i] + Gamzxz[i] * gyzz[i]
                );

            Rxz[i] =
                HALF * (
                    -Rxz[i]
                    + (dxx[i] + ONE) * Gamxz[i] + gxy[i] * Gamyz[i] + gxz[i] * Gamzz[i]
                    + gxz[i] * Gamxx[i] + gyz[i] * Gamyx[i] + (dzz[i] + ONE) * Gamzx[i]
                    + gxa * gxzx[i] + gya * gyzx[i] + gza * gzzx[i]
                    + gxa * gxxz[i] + gya * gxyz[i] + gza * gxzz[i]
                )
                + gupxx[i] * (
                    Gamxxx[i] * gxxz[i] + Gamyxx[i] * gxyz[i] + Gamzxx[i] * gxzz[i]
                    + Gamxxz[i] * gxxx[i] + Gamyxz[i] * gxyx[i] + Gamzxz[i] * gxzx[i]
                    + Gamxxx[i] * gxzx[i] + Gamyxx[i] * gxzy[i] + Gamzxx[i] * gxzz[i]
                )
                + gupxy[i] * (
                    Gamxxx[i] * gxyz[i] + Gamyxx[i] * gyyz[i] + Gamzxx[i] * gyzz[i]
                    + Gamxxz[i] * gxyx[i] + Gamyxz[i] * gyyx[i] + Gamzxz[i] * gyzx[i]
                    + Gamxxy[i] * gxzx[i] + Gamyxy[i] * gxzy[i] + Gamzxy[i] * gxzz[i]
                    + Gamxxy[i] * gxxz[i] + Gamyxy[i] * gxyz[i] + Gamzxy[i] * gxzz[i]
                    + Gamxyz[i] * gxxx[i] + Gamyyz[i] * gxyx[i] + Gamzyz[i] * gxzx[i]
                    + Gamxxx[i] * gyzx[i] + Gamyxx[i] * gyzy[i] + Gamzxx[i] * gyzz[i]
                )
                + gupxz[i] * (
                    Gamxxx[i] * gxzz[i] + Gamyxx[i] * gyzz[i] + Gamzxx[i] * gzzz[i]
                    + Gamxxz[i] * gxzx[i] + Gamyxz[i] * gyzx[i] + Gamzxz[i] * gzzx[i]
                    + Gamxxz[i] * gxzx[i] + Gamyxz[i] * gxzy[i] + Gamzxz[i] * gxzz[i]
                    + Gamxxz[i] * gxxz[i] + Gamyxz[i] * gxyz[i] + Gamzxz[i] * gxzz[i]
                    + Gamxzz[i] * gxxx[i] + Gamyzz[i] * gxyx[i] + Gamzzz[i] * gxzx[i]
                    + Gamxxx[i] * gzzx[i] + Gamyxx[i] * gzzy[i] + Gamzxx[i] * gzzz[i]
                )
                + gupyy[i] * (
                    Gamxxy[i] * gxyz[i] + Gamyxy[i] * gyyz[i] + Gamzxy[i] * gyzz[i]
                    + Gamxyz[i] * gxyx[i] + Gamyyz[i] * gyyx[i] + Gamzyz[i] * gyzx[i]
                    + Gamxxy[i] * gyzx[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gyzz[i]
                )
                + gupyz[i] * (
                    Gamxxy[i] * gxzz[i] + Gamyxy[i] * gyzz[i] + Gamzxy[i] * gzzz[i]
                    + Gamxyz[i] * gxzx[i] + Gamyyz[i] * gyzx[i] + Gamzyz[i] * gzzx[i]
                    + Gamxxz[i] * gyzx[i] + Gamyxz[i] * gyzy[i] + Gamzxz[i] * gyzz[i]
                    + Gamxxz[i] * gxyz[i] + Gamyxz[i] * gyyz[i] + Gamzxz[i] * gyzz[i]
                    + Gamxzz[i] * gxyx[i] + Gamyzz[i] * gyyx[i] + Gamzzz[i] * gyzx[i]
                    + Gamxxy[i] * gzzx[i] + Gamyxy[i] * gzzy[i] + Gamzxy[i] * gzzz[i]
                )
                + gupzz[i] * (
                    Gamxxz[i] * gxzz[i] + Gamyxz[i] * gyzz[i] + Gamzxz[i] * gzzz[i]
                    + Gamxzz[i] * gxzx[i] + Gamyzz[i] * gyzx[i] + Gamzzz[i] * gzzx[i]
                    + Gamxxz[i] * gzzx[i] + Gamyxz[i] * gzzy[i] + Gamzxz[i] * gzzz[i]
                );

            Ryz[i] =
                HALF * (
                    -Ryz[i]
                    + gxy[i] * Gamxz[i] + (dyy[i] + ONE) * Gamyz[i] + gyz[i] * Gamzz[i]
                    + gxz[i] * Gamxy[i] + gyz[i] * Gamyy[i] + (dzz[i] + ONE) * Gamzy[i]
                    + gxa * gxzy[i] + gya * gyzy[i] + gza * gzzy[i]
                    + gxa * gxyz[i] + gya * gyyz[i] + gza * gyzz[i]
                )
                + gupxx[i] * (
                    Gamxxy[i] * gxxz[i] + Gamyxy[i] * gxyz[i] + Gamzxy[i] * gxzz[i]
                    + Gamxxz[i] * gxxy[i] + Gamyxz[i] * gxyy[i] + Gamzxz[i] * gxzy[i]
                    + Gamxxy[i] * gxzx[i] + Gamyxy[i] * gxzy[i] + Gamzxy[i] * gxzz[i]
                )
                + gupxy[i] * (
                    Gamxxy[i] * gxyz[i] + Gamyxy[i] * gyyz[i] + Gamzxy[i] * gyzz[i]
                    + Gamxxz[i] * gxyy[i] + Gamyxz[i] * gyyy[i] + Gamzxz[i] * gyzy[i]
                    + Gamxyy[i] * gxzx[i] + Gamyyy[i] * gxzy[i] + Gamzyy[i] * gxzz[i]
                    + Gamxyy[i] * gxxz[i] + Gamyyy[i] * gxyz[i] + Gamzyy[i] * gxzz[i]
                    + Gamxyz[i] * gxxy[i] + Gamyyz[i] * gxyy[i] + Gamzyz[i] * gxzy[i]
                    + Gamxxy[i] * gyzx[i] + Gamyxy[i] * gyzy[i] + Gamzxy[i] * gyzz[i]
                )
                + gupxz[i] * (
                    Gamxxy[i] * gxzz[i] + Gamyxy[i] * gyzz[i] + Gamzxy[i] * gzzz[i]
                    + Gamxxz[i] * gxzy[i] + Gamyxz[i] * gyzy[i] + Gamzxz[i] * gzzy[i]
                    + Gamxyz[i] * gxzx[i] + Gamyyz[i] * gxzy[i] + Gamzyz[i] * gxzz[i]
                    + Gamxyz[i] * gxxz[i] + Gamyyz[i] * gxyz[i] + Gamzyz[i] * gxzz[i]
                    + Gamxzz[i] * gxxy[i] + Gamyzz[i] * gxyy[i] + Gamzzz[i] * gxzy[i]
                    + Gamxxy[i] * gzzx[i] + Gamyxy[i] * gzzy[i] + Gamzxy[i] * gzzz[i]
                )
                + gupyy[i] * (
                    Gamxyy[i] * gxyz[i] + Gamyyy[i] * gyyz[i] + Gamzyy[i] * gyzz[i]
                    + Gamxyz[i] * gxyy[i] + Gamyyz[i] * gyyy[i] + Gamzyz[i] * gyzy[i]
                    + Gamxyy[i] * gyzx[i] + Gamyyy[i] * gyzy[i] + Gamzyy[i] * gyzz[i]
                )
                + gupyz[i] * (
                    Gamxyy[i] * gxzz[i] + Gamyyy[i] * gyzz[i] + Gamzyy[i] * gzzz[i]
                    + Gamxyz[i] * gxzy[i] + Gamyyz[i] * gyzy[i] + Gamzyz[i] * gzzy[i]
                    + Gamxyz[i] * gyzx[i] + Gamyyz[i] * gyzy[i] + Gamzyz[i] * gyzz[i]
                    + Gamxyz[i] * gxyz[i] + Gamyyz[i] * gyyz[i] + Gamzyz[i] * gyzz[i]
                    + Gamxzz[i] * gxyy[i] + Gamyzz[i] * gyyy[i] + Gamzzz[i] * gyzy[i]
                    + Gamxyy[i] * gzzx[i] + Gamyyy[i] * gzzy[i] + Gamzyy[i] * gzzz[i]
                )
                + gupzz[i] * (
                    Gamxyz[i] * gxzz[i] + Gamyyz[i] * gyzz[i] + Gamzyz[i] * gzzz[i]
                    + Gamxzz[i] * gxzy[i] + Gamyzz[i] * gyzy[i] + Gamzzz[i] * gzzy[i]
                    + Gamxyz[i] * gzzx[i] + Gamyyz[i] * gzzy[i] + Gamzyz[i] * gzzz[i]
                );
        }
        RHS_KERNEL_TIMER_ADD(KB_RICCI_METRIC, timer_ricci_metric);

        RHS_KERNEL_TIMER_DECL(timer_chi_lapse);
        // 22.3ms //
            fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);

        // 7ms //
        for (int i=0;i<all;i+=1) {
            const double inv_chin1 = ONE / chin1[i];
            const double half_inv_chin1 = HALF * inv_chin1;
            const double scaled_inv = F3o2 * inv_chin1;
            const double cxx = fxx[i] - Gamxxx[i] * chix[i] - Gamyxx[i] * chiy[i] - Gamzxx[i] * chiz[i];
            const double cxy = fxy[i] - Gamxxy[i] * chix[i] - Gamyxy[i] * chiy[i] - Gamzxy[i] * chiz[i];
            const double cxz = fxz[i] - Gamxxz[i] * chix[i] - Gamyxz[i] * chiy[i] - Gamzxz[i] * chiz[i];
            const double cyy = fyy[i] - Gamxyy[i] * chix[i] - Gamyyy[i] * chiy[i] - Gamzyy[i] * chiz[i];
            const double cyz = fyz[i] - Gamxyz[i] * chix[i] - Gamyyz[i] * chiy[i] - Gamzyz[i] * chiz[i];
            const double czz = fzz[i] - Gamxzz[i] * chix[i] - Gamyzz[i] * chiy[i] - Gamzzz[i] * chiz[i];
            const double ricci_chi =
                gupxx[i] * (cxx - scaled_inv * chix[i] * chix[i])
                + gupyy[i] * (cyy - scaled_inv * chiy[i] * chiy[i])
                + gupzz[i] * (czz - scaled_inv * chiz[i] * chiz[i])
                + TWO * gupxy[i] * (cxy - scaled_inv * chix[i] * chiy[i])
                + TWO * gupxz[i] * (cxz - scaled_inv * chix[i] * chiz[i])
                + TWO * gupyz[i] * (cyz - scaled_inv * chiy[i] * chiz[i]);
            f[i] = ricci_chi;
            Rxx[i] = Rxx[i] + ( cxx - half_inv_chin1 * chix[i] * chix[i] + (dxx[i] + ONE) * ricci_chi ) * half_inv_chin1;
            Ryy[i] = Ryy[i] + ( cyy - half_inv_chin1 * chiy[i] * chiy[i] + (dyy[i] + ONE) * ricci_chi ) * half_inv_chin1;
            Rzz[i] = Rzz[i] + ( czz - half_inv_chin1 * chiz[i] * chiz[i] + (dzz[i] + ONE) * ricci_chi ) * half_inv_chin1;

            Rxy[i] = Rxy[i] + ( cxy - half_inv_chin1 * chix[i] * chiy[i] + gxy[i] * ricci_chi ) * half_inv_chin1;
            Rxz[i] = Rxz[i] + ( cxz - half_inv_chin1 * chix[i] * chiz[i] + gxz[i] * ricci_chi ) * half_inv_chin1;
            Ryz[i] = Ryz[i] + ( cyz - half_inv_chin1 * chiy[i] * chiz[i] + gyz[i] * ricci_chi ) * half_inv_chin1;
        }

        // 24ms //
        fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);

        // 6ms //
        for (int i=0;i<all;i+=1) {
            const double inv_chin1 = ONE / chin1[i];
            const double gchi_x = (gupxx[i] * chix[i] + gupxy[i] * chiy[i] + gupxz[i] * chiz[i]) * inv_chin1;
            const double gchi_y = (gupxy[i] * chix[i] + gupyy[i] * chiy[i] + gupyz[i] * chiz[i]) * inv_chin1;
            const double gchi_z = (gupxz[i] * chix[i] + gupyz[i] * chiy[i] + gupzz[i] * chiz[i]) * inv_chin1;

            /* Christoffel 修正项 */
            Gamxxx[i] = Gamxxx[i] - ( ((chix[i] + chix[i]) * inv_chin1) - (dxx[i] + ONE) * gchi_x ) * HALF;
            Gamyxx[i] = Gamyxx[i] - ( 0.0           - (dxx[i] + ONE) * gchi_y ) * HALF;  /* 原式只有 -gxx*gxxy */
            Gamzxx[i] = Gamzxx[i] - ( 0.0           - (dxx[i] + ONE) * gchi_z ) * HALF;

            Gamxyy[i] = Gamxyy[i] - ( 0.0           - (dyy[i] + ONE) * gchi_x ) * HALF;
            Gamyyy[i] = Gamyyy[i] - ( ((chiy[i] + chiy[i]) * inv_chin1) - (dyy[i] + ONE) * gchi_y ) * HALF;
            Gamzyy[i] = Gamzyy[i] - ( 0.0           - (dyy[i] + ONE) * gchi_z ) * HALF;

            Gamxzz[i] = Gamxzz[i] - ( 0.0           - (dzz[i] + ONE) * gchi_x ) * HALF;
            Gamyzz[i] = Gamyzz[i] - ( 0.0           - (dzz[i] + ONE) * gchi_y ) * HALF;
            Gamzzz[i] = Gamzzz[i] - ( ((chiz[i] + chiz[i]) * inv_chin1) - (dzz[i] + ONE) * gchi_z ) * HALF;

            Gamxxy[i] = Gamxxy[i] - ( ( chiy[i] * inv_chin1)      - gxy[i] * gchi_x ) * HALF;
            Gamyxy[i] = Gamyxy[i] - ( ( chix[i] * inv_chin1)      - gxy[i] * gchi_y ) * HALF;
            Gamzxy[i] = Gamzxy[i] - ( 0.0           - gxy[i] * gchi_z ) * HALF;

            Gamxxz[i] = Gamxxz[i] - ( ( chiz[i] * inv_chin1)      - gxz[i] * gchi_x ) * HALF;
            Gamyxz[i] = Gamyxz[i] - ( 0.0           - gxz[i] * gchi_y ) * HALF;
            Gamzxz[i] = Gamzxz[i] - ( ( chix[i] * inv_chin1)      - gxz[i] * gchi_z ) * HALF;

            Gamxyz[i] = Gamxyz[i] - ( 0.0           - gyz[i] * gchi_x ) * HALF;
            Gamyyz[i] = Gamyyz[i] - ( ( chiz[i] * inv_chin1)      - gyz[i] * gchi_y ) * HALF;
            Gamzyz[i] = Gamzyz[i] - ( ( chiy[i] * inv_chin1)      - gyz[i] * gchi_z ) * HALF;

            /* fxx..fyz 修正：减去 Γ * ∂Lap */
            fxx[i] = fxx[i] - Gamxxx[i] * Lapx[i] - Gamyxx[i] * Lapy[i] - Gamzxx[i] * Lapz[i];
            fyy[i] = fyy[i] - Gamxyy[i] * Lapx[i] - Gamyyy[i] * Lapy[i] - Gamzyy[i] * Lapz[i];
            fzz[i] = fzz[i] - Gamxzz[i] * Lapx[i] - Gamyzz[i] * Lapy[i] - Gamzzz[i] * Lapz[i];

            fxy[i] = fxy[i] - Gamxxy[i] * Lapx[i] - Gamyxy[i] * Lapy[i] - Gamzxy[i] * Lapz[i];
            fxz[i] = fxz[i] - Gamxxz[i] * Lapx[i] - Gamyxz[i] * Lapy[i] - Gamzxz[i] * Lapz[i];
            fyz[i] = fyz[i] - Gamxyz[i] * Lapx[i] - Gamyyz[i] * Lapy[i] - Gamzyz[i] * Lapz[i];

            trK_rhs[i] =  gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                    +  TWO * ( gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i] );
        }
        RHS_KERNEL_TIMER_ADD(KB_CHI_LAPSE, timer_chi_lapse);
        RHS_KERNEL_TIMER_DECL(timer_aij_trk_gauge);
        // 2.5ms //
        for (int i=0;i<all;i+=1) {
            const double divb = betaxx[i] + betayy[i] + betazz[i];

            S[i] = chin1[i] * (
                    gupxx[i] * Sxx[i] + gupyy[i] * Syy[i] + gupzz[i] * Szz[i]
                + TWO * ( gupxy[i] * Sxy[i] + gupxz[i] * Sxz[i] + gupyz[i] * Syz[i] )
            );

            f[i] = F2o3 * trK[i] * trK[i]
                - (
                    gupxx[i] * (
                        gupxx[i] * Axx[i] * Axx[i] + gupyy[i] * Axy[i] * Axy[i] + gupzz[i] * Axz[i] * Axz[i]
                        + TWO * ( gupxy[i] * Axx[i] * Axy[i] + gupxz[i] * Axx[i] * Axz[i] + gupyz[i] * Axy[i] * Axz[i] )
                    )
                    + gupyy[i] * (
                        gupxx[i] * Axy[i] * Axy[i] + gupyy[i] * Ayy[i] * Ayy[i] + gupzz[i] * Ayz[i] * Ayz[i]
                        + TWO * ( gupxy[i] * Axy[i] * Ayy[i] + gupxz[i] * Axy[i] * Ayz[i] + gupyz[i] * Ayy[i] * Ayz[i] )
                    )
                    + gupzz[i] * (
                        gupxx[i] * Axz[i] * Axz[i] + gupyy[i] * Ayz[i] * Ayz[i] + gupzz[i] * Azz[i] * Azz[i]
                        + TWO * ( gupxy[i] * Axz[i] * Ayz[i] + gupxz[i] * Axz[i] * Azz[i] + gupyz[i] * Ayz[i] * Azz[i] )
                    )
                    + TWO * (
                        gupxy[i] * (
                            gupxx[i] * Axx[i] * Axy[i] + gupyy[i] * Axy[i] * Ayy[i] + gupzz[i] * Axz[i] * Ayz[i]
                            + gupxy[i] * ( Axx[i] * Ayy[i] + Axy[i] * Axy[i] )
                            + gupxz[i] * ( Axx[i] * Ayz[i] + Axz[i] * Axy[i] )
                            + gupyz[i] * ( Axy[i] * Ayz[i] + Axz[i] * Ayy[i] )
                        )
                        + gupxz[i] * (
                            gupxx[i] * Axx[i] * Axz[i] + gupyy[i] * Axy[i] * Ayz[i] + gupzz[i] * Axz[i] * Azz[i]
                            + gupxy[i] * ( Axx[i] * Ayz[i] + Axy[i] * Axz[i] )
                            + gupxz[i] * ( Axx[i] * Azz[i] + Axz[i] * Axz[i] )
                            + gupyz[i] * ( Axy[i] * Azz[i] + Axz[i] * Ayz[i] )
                        )
                        + gupyz[i] * (
                            gupxx[i] * Axy[i] * Axz[i] + gupyy[i] * Ayy[i] * Ayz[i] + gupzz[i] * Ayz[i] * Azz[i]
                            + gupxy[i] * ( Axy[i] * Ayz[i] + Ayy[i] * Axz[i] )
                            + gupxz[i] * ( Axy[i] * Azz[i] + Ayz[i] * Axz[i] )
                            + gupyz[i] * ( Ayy[i] * Azz[i] + Ayz[i] * Ayz[i] )
                        )
                    )
                )
                - F16 * PI * rho[i]
                + EIGHT * PI * S[i];

            f[i] = -F1o3 * (
                    gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                + TWO * ( gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i] )
                + (alpn1[i] / chin1[i]) * f[i]
            );

            double l_fxx = alpn1[i] * (Rxx[i] - EIGHT * PI * Sxx[i]) - fxx[i];
            double l_fxy = alpn1[i] * (Rxy[i] - EIGHT * PI * Sxy[i]) - fxy[i];
            double l_fxz = alpn1[i] * (Rxz[i] - EIGHT * PI * Sxz[i]) - fxz[i];
            double l_fyy = alpn1[i] * (Ryy[i] - EIGHT * PI * Syy[i]) - fyy[i];
            double l_fyz = alpn1[i] * (Ryz[i] - EIGHT * PI * Syz[i]) - fyz[i];
            double l_fzz = alpn1[i] * (Rzz[i] - EIGHT * PI * Szz[i]) - fzz[i];

            /* Aij_rhs = fij - gij * f */
            Axx_rhs[i] = l_fxx - (dxx[i] + ONE) * f[i];
            Ayy_rhs[i] = l_fyy - (dyy[i] + ONE) * f[i];
            Azz_rhs[i] = l_fzz - (dzz[i] + ONE) * f[i];
            Axy_rhs[i] = l_fxy - gxy[i] * f[i];
            Axz_rhs[i] = l_fxz - gxz[i] * f[i];
            Ayz_rhs[i] = l_fyz - gyz[i] * f[i];

            /* Now: store A_il A^l_j into fij: */
            fxx[i] =
                gupxx[i] * Axx[i] * Axx[i]
                + gupyy[i] * Axy[i] * Axy[i]
                + gupzz[i] * Axz[i] * Axz[i]
                + TWO * ( gupxy[i] * Axx[i] * Axy[i]
                        + gupxz[i] * Axx[i] * Axz[i]
                        + gupyz[i] * Axy[i] * Axz[i] );

            fyy[i] =
                gupxx[i] * Axy[i] * Axy[i]
                + gupyy[i] * Ayy[i] * Ayy[i]
                + gupzz[i] * Ayz[i] * Ayz[i]
                + TWO * ( gupxy[i] * Axy[i] * Ayy[i]
                        + gupxz[i] * Axy[i] * Ayz[i]
                        + gupyz[i] * Ayy[i] * Ayz[i] );

            fzz[i] =
                gupxx[i] * Axz[i] * Axz[i]
                + gupyy[i] * Ayz[i] * Ayz[i]
                + gupzz[i] * Azz[i] * Azz[i]
                + TWO * ( gupxy[i] * Axz[i] * Ayz[i]
                        + gupxz[i] * Axz[i] * Azz[i]
                        + gupyz[i] * Ayz[i] * Azz[i] );

            fxy[i] =
                gupxx[i] * Axx[i] * Axy[i]
                + gupyy[i] * Axy[i] * Ayy[i]
                + gupzz[i] * Axz[i] * Ayz[i]
                + gupxy[i] * (Axx[i] * Ayy[i] + Axy[i] * Axy[i])
                + gupxz[i] * (Axx[i] * Ayz[i] + Axz[i] * Axy[i])
                + gupyz[i] * (Axy[i] * Ayz[i] + Axz[i] * Ayy[i]);

            fxz[i] =
                gupxx[i] * Axx[i] * Axz[i]
                + gupyy[i] * Axy[i] * Ayz[i]
                + gupzz[i] * Axz[i] * Azz[i]
                + gupxy[i] * (Axx[i] * Ayz[i] + Axy[i] * Axz[i])
                + gupxz[i] * (Axx[i] * Azz[i] + Axz[i] * Axz[i])
                + gupyz[i] * (Axy[i] * Azz[i] + Axz[i] * Ayz[i]);

            fyz[i] =
                gupxx[i] * Axy[i] * Axz[i]
                + gupyy[i] * Ayy[i] * Ayz[i]
                + gupzz[i] * Ayz[i] * Azz[i]
                + gupxy[i] * (Axy[i] * Ayz[i] + Ayy[i] * Axz[i])
                + gupxz[i] * (Axy[i] * Azz[i] + Ayz[i] * Axz[i])
                + gupyz[i] * (Ayy[i] * Azz[i] + Ayz[i] * Ayz[i]);

            /* f = chin1 */
            f[i] = chin1[i];

            /* store D^i D_i Lap in trK_rhs */
            trK_rhs[i] = f[i] * trK_rhs[i];

            /* rhs for Aij */
            Axx_rhs[i] =
                f[i] * Axx_rhs[i]
                + alpn1[i] * ( trK[i] * Axx[i] - TWO * fxx[i] )
                + TWO * ( Axx[i] * betaxx[i] + Axy[i] * betayx[i] + Axz[i] * betazx[i] )
                - F2o3 * Axx[i] * divb;

            Ayy_rhs[i] =
                f[i] * Ayy_rhs[i]
                + alpn1[i] * ( trK[i] * Ayy[i] - TWO * fyy[i] )
                + TWO * ( Axy[i] * betaxy[i] + Ayy[i] * betayy[i] + Ayz[i] * betazy[i] )
                - F2o3 * Ayy[i] * divb;

            Azz_rhs[i] =
                f[i] * Azz_rhs[i]
                + alpn1[i] * ( trK[i] * Azz[i] - TWO * fzz[i] )
                + TWO * ( Axz[i] * betaxz[i] + Ayz[i] * betayz[i] + Azz[i] * betazz[i] )
                - F2o3 * Azz[i] * divb;

            Axy_rhs[i] =
                f[i] * Axy_rhs[i]
                + alpn1[i] * ( trK[i] * Axy[i] - TWO * fxy[i] )
                + Axx[i] * betaxy[i]
                + Axz[i] * betazy[i]
                + Ayy[i] * betayx[i]
                + Ayz[i] * betazx[i]
                + F1o3 * Axy[i] * divb
                - Axy[i] * betazz[i];

            Ayz_rhs[i] =
                f[i] * Ayz_rhs[i]
                + alpn1[i] * ( trK[i] * Ayz[i] - TWO * fyz[i] )
                + Axy[i] * betaxz[i]
                + Ayy[i] * betayz[i]
                + Axz[i] * betaxy[i]
                + Azz[i] * betazy[i]
                + F1o3 * Ayz[i] * divb
                - Ayz[i] * betaxx[i];

            Axz_rhs[i] =
                f[i] * Axz_rhs[i]
                + alpn1[i] * ( trK[i] * Axz[i] - TWO * fxz[i] )
                + Axx[i] * betaxz[i]
                + Axy[i] * betayz[i]
                + Ayz[i] * betayx[i]
                + Azz[i] * betazx[i]
                + F1o3 * Axz[i] * divb
                - Axz[i] * betayy[i];

            /* Compute trace of S_ij */
            S[i] =
                f[i] * ( gupxx[i] * Sxx[i] + gupyy[i] * Syy[i] + gupzz[i] * Szz[i]
                        + TWO * ( gupxy[i] * Sxy[i] + gupxz[i] * Sxz[i] + gupyz[i] * Syz[i] ) );

            /* rhs for trK */
            trK_rhs[i] =
                -trK_rhs[i]
                + alpn1[i] * ( F1o3 * trK[i] * trK[i]
                            + gupxx[i] * fxx[i] + gupyy[i] * fyy[i] + gupzz[i] * fzz[i]
                            + TWO * ( gupxy[i] * fxy[i] + gupxz[i] * fxz[i] + gupyz[i] * fyz[i] )
                            + FOUR * PI * ( rho[i] + S[i] ) );

            /* gauge variable part */
            Lap_rhs[i] = -TWO * alpn1[i] * trK[i];
        }

        // 1ms //
        #if (GAUGE == 6 || GAUGE == 7)
        if (BHN == 2) {
            const double M = Mass[0] + Mass[1];
            const double A = TWO / M;
            const double w1 = 1.2e1;
            const double w2 = w1;
            const double C1 = ONE / Mass[0] - A;
            const double C2 = ONE / Mass[1] - A;
            const double denom =
                (Porg[0] - Porg[3]) * (Porg[0] - Porg[3]) +
                (Porg[1] - Porg[4]) * (Porg[1] - Porg[4]) +
                (Porg[2] - Porg[5]) * (Porg[2] - Porg[5]);

            for (int k0 = 0; k0 < nz; ++k0) {
                for (int j0 = 0; j0 < ny; ++j0) {
                    for (int i0 = 0; i0 < nx; ++i0) {
                        const size_t p = idx_ex(i0, j0, k0, ex);
                        const double dx1 = Porg[0] - X[i0];
                        const double dy1 = Porg[1] - Y[j0];
                        const double dz1 = Porg[2] - Z[k0];
                        const double dx2 = Porg[3] - X[i0];
                        const double dy2 = Porg[4] - Y[j0];
                        const double dz2 = Porg[5] - Z[k0];
                        const double r1 = (dx1 * dx1 + dy1 * dy1 + dz1 * dz1) / denom;
                        const double r2 = (dx2 * dx2 + dy2 * dy2 + dz2 * dz2) / denom;
                        #if (GAUGE == 6)
                        reta[p] = A + C1 / (ONE + w1 * r1) + C2 / (ONE + w2 * r2);
                        #else
                        reta[p] = A + C1 * exp(-w1 * r1) + C2 * exp(-w2 * r2);
                        #endif
                    }
                }
            }
        } else {
            printf("not support BH_num in Jason's form %d %d\n", (GAUGE == 6) ? 1 : 2, BHN);
            for (int i = 0; i < all; ++i) reta[i] = ZEO;
        }
        #endif

        for (int i = 0; i < all; i += 1) {
            #if (GAUGE == 0)
            betax_rhs[i] = FF * dtSfx[i];
            betay_rhs[i] = FF * dtSfy[i];
            betaz_rhs[i] = FF * dtSfz[i];

            dtSfx_rhs[i] = Gamx_rhs[i] - eta * dtSfx[i];
            dtSfy_rhs[i] = Gamy_rhs[i] - eta * dtSfy[i];
            dtSfz_rhs[i] = Gamz_rhs[i] - eta * dtSfz[i];
            #elif (GAUGE == 1)
            betax_rhs[i] = Gamx[i] - eta * betax[i];
            betay_rhs[i] = Gamy[i] - eta * betay[i];
            betaz_rhs[i] = Gamz[i] - eta * betaz[i];

            dtSfx_rhs[i] = ZEO;
            dtSfy_rhs[i] = ZEO;
            dtSfz_rhs[i] = ZEO;
            #elif (GAUGE == 2 || GAUGE == 3)
            betax_rhs[i] = FF * dtSfx[i];
            betay_rhs[i] = FF * dtSfy[i];
            betaz_rhs[i] = FF * dtSfz[i];

            reta[i] =
                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );

            #if (GAUGE == 2)
            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
            #else
            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
            #endif

            dtSfx_rhs[i] = Gamx_rhs[i] - reta[i] * dtSfx[i];
            dtSfy_rhs[i] = Gamy_rhs[i] - reta[i] * dtSfy[i];
            dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
            #elif (GAUGE == 4 || GAUGE == 5)
            reta[i] =
                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );

            #if (GAUGE == 4)
            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
            #else
            reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - chin1[i]), 2.0 );
            #endif

            betax_rhs[i] = FF * Gamx[i] - reta[i] * betax[i];
            betay_rhs[i] = FF * Gamy[i] - reta[i] * betay[i];
            betaz_rhs[i] = FF * Gamz[i] - reta[i] * betaz[i];

            dtSfx_rhs[i] = ZEO;
            dtSfy_rhs[i] = ZEO;
            dtSfz_rhs[i] = ZEO;
            #elif (GAUGE == 6 || GAUGE == 7)
            betax_rhs[i] = FF * dtSfx[i];
            betay_rhs[i] = FF * dtSfy[i];
            betaz_rhs[i] = FF * dtSfz[i];

            dtSfx_rhs[i] = Gamx_rhs[i] - reta[i] * dtSfx[i];
            dtSfy_rhs[i] = Gamy_rhs[i] - reta[i] * dtSfy[i];
            dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
            #endif
        }
        RHS_KERNEL_TIMER_ADD(KB_AIJ_TRK_GAUGE, timer_aij_trk_gauge);
        RHS_KERNEL_TIMER_DECL(timer_ko_constraint);
        // advection + KO dissipation with shared symmetry buffer
        lopsided_kodis(ex,X,Y,Z,dxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
        lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps);
        lopsided_kodis(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps);
        lopsided_kodis(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,eps);
        lopsided_kodis(ex,X,Y,Z,dyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,eps);
        lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps);
        lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
        lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps);
        lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps);
        lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps);
        lopsided_kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,eps);
        lopsided_kodis(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,eps);
        lopsided_kodis(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,eps);
        lopsided_kodis(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,eps);
        lopsided_kodis(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,eps);
        lopsided_kodis(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,eps);
        // 2ms //
        if(co==0){
            for (int i=0;i<all;i+=1) {
                ham_Res[i] =
                    gupxx[i] * Rxx[i] + gupyy[i] * Ryy[i] + gupzz[i] * Rzz[i]
                    + TWO * ( gupxy[i] * Rxy[i] + gupxz[i] * Rxz[i] + gupyz[i] * Ryz[i] );
                ham_Res[i] =
                    chin1[i] * ham_Res[i]
                    + F2o3 * trK[i] * trK[i]
                    - (
                        gupxx[i] * (
                            gupxx[i] * Axx[i] * Axx[i] + gupyy[i] * Axy[i] * Axy[i] + gupzz[i] * Axz[i] * Axz[i]
                        + TWO * ( gupxy[i] * Axx[i] * Axy[i] + gupxz[i] * Axx[i] * Axz[i] + gupyz[i] * Axy[i] * Axz[i] )
                        )
                    + gupyy[i] * (
                            gupxx[i] * Axy[i] * Axy[i] + gupyy[i] * Ayy[i] * Ayy[i] + gupzz[i] * Ayz[i] * Ayz[i]
                        + TWO * ( gupxy[i] * Axy[i] * Ayy[i] + gupxz[i] * Axy[i] * Ayz[i] + gupyz[i] * Ayy[i] * Ayz[i] )
                        )
                    + gupzz[i] * (
                            gupxx[i] * Axz[i] * Axz[i] + gupyy[i] * Ayz[i] * Ayz[i] + gupzz[i] * Azz[i] * Azz[i]
                        + TWO * ( gupxy[i] * Axz[i] * Ayz[i] + gupxz[i] * Axz[i] * Azz[i] + gupyz[i] * Ayz[i] * Azz[i] )
                        )
                    + TWO * (
                            gupxy[i] * (
                                gupxx[i] * Axx[i] * Axy[i] + gupyy[i] * Axy[i] * Ayy[i] + gupzz[i] * Axz[i] * Ayz[i]
                            + gupxy[i] * ( Axx[i] * Ayy[i] + Axy[i] * Axy[i] )
                            + gupxz[i] * ( Axx[i] * Ayz[i] + Axz[i] * Axy[i] )
                            + gupyz[i] * ( Axy[i] * Ayz[i] + Axz[i] * Ayy[i] )
                            )
                        + gupxz[i] * (
                                gupxx[i] * Axx[i] * Axz[i] + gupyy[i] * Axy[i] * Ayz[i] + gupzz[i] * Axz[i] * Azz[i]
                            + gupxy[i] * ( Axx[i] * Ayz[i] + Axy[i] * Axz[i] )
                            + gupxz[i] * ( Axx[i] * Azz[i] + Axz[i] * Axz[i] )
                            + gupyz[i] * ( Axy[i] * Azz[i] + Axz[i] * Ayz[i] )
                            )
                        + gupyz[i] * (
                                gupxx[i] * Axy[i] * Axz[i] + gupyy[i] * Ayy[i] * Ayz[i] + gupzz[i] * Ayz[i] * Azz[i]
                            + gupxy[i] * ( Axy[i] * Ayz[i] + Ayy[i] * Axz[i] )
                            + gupxz[i] * ( Axy[i] * Azz[i] + Ayz[i] * Axz[i] )
                            + gupyz[i] * ( Ayy[i] * Azz[i] + Ayz[i] * Ayz[i] )
                            )
                        )
                    ) - F16 * PI * rho[i];
            }

            // 1ms //
            fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
            fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0);
            fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0);
            fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
            fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0);
            fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
            // 7ms //
            for (int i=0;i<all;i+=1) {
                gxxx[i] = gxxx[i] - (  Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]
                                      + Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]) - chix[i]*Axx[i]/chin1[i];
                gxyx[i] = gxyx[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
                                      + Gamxxx[i] * Axy[i] + Gamyxx[i] * Ayy[i] + Gamzxx[i] * Ayz[i]) - chix[i]*Axy[i]/chin1[i];
                gxzx[i] = gxzx[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
                                      + Gamxxx[i] * Axz[i] + Gamyxx[i] * Ayz[i] + Gamzxx[i] * Azz[i]) - chix[i]*Axz[i]/chin1[i];
                gyyx[i] = gyyx[i] - (  Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]
                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chix[i]*Ayy[i]/chin1[i];
                gyzx[i] = gyzx[i] - (  Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]
                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chix[i]*Ayz[i]/chin1[i];
                gzzx[i] = gzzx[i] - (  Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]
                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chix[i]*Azz[i]/chin1[i];
                gxxy[i] = gxxy[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
                                      + Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]) - chiy[i]*Axx[i]/chin1[i];
                gxyy[i] = gxyy[i] - (  Gamxyy[i] * Axx[i] + Gamyyy[i] * Axy[i] + Gamzyy[i] * Axz[i]
                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chiy[i]*Axy[i]/chin1[i];
                gxzy[i] = gxzy[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chiy[i]*Axz[i]/chin1[i];
                gyyy[i] = gyyy[i] - (  Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]
                                      + Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]) - chiy[i]*Ayy[i]/chin1[i];
                gyzy[i] = gyzy[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
                                      + Gamxyy[i] * Axz[i] + Gamyyy[i] * Ayz[i] + Gamzyy[i] * Azz[i]) - chiy[i]*Ayz[i]/chin1[i];
                gzzy[i] = gzzy[i] - (  Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]
                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiy[i]*Azz[i]/chin1[i];
                gxxz[i] = gxxz[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
                                      + Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]) - chiz[i]*Axx[i]/chin1[i];
                gxyz[i] = gxyz[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
                                      + Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]) - chiz[i]*Axy[i]/chin1[i];
                gxzz[i] = gxzz[i] - (  Gamxzz[i] * Axx[i] + Gamyzz[i] * Axy[i] + Gamzzz[i] * Axz[i]
                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chiz[i]*Axz[i]/chin1[i];
                gyyz[i] = gyyz[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
                                      + Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]) - chiz[i]*Ayy[i]/chin1[i];
                gyzz[i] = gyzz[i] - (  Gamxzz[i] * Axy[i] + Gamyzz[i] * Ayy[i] + Gamzzz[i] * Ayz[i]
                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiz[i]*Ayz[i]/chin1[i];
                gzzz[i] = gzzz[i] - (  Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]
                                      + Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]) - chiz[i]*Azz[i]/chin1[i];

                movx_Res[i] = gupxx[i]*gxxx[i] + gupyy[i]*gxyy[i] + gupzz[i]*gxzz[i]
                            + gupxy[i]*gxyx[i] + gupxz[i]*gxzx[i] + gupyz[i]*gxzy[i]
                            + gupxy[i]*gxxy[i] + gupxz[i]*gxxz[i] + gupyz[i]*gxyz[i];
                movy_Res[i] = gupxx[i]*gxyx[i] + gupyy[i]*gyyy[i] + gupzz[i]*gyzz[i]
                            + gupxy[i]*gyyx[i] + gupxz[i]*gyzx[i] + gupyz[i]*gyzy[i]
                            + gupxy[i]*gxyy[i] + gupxz[i]*gxyz[i] + gupyz[i]*gyyz[i];
                movz_Res[i] = gupxx[i]*gxzx[i] + gupyy[i]*gyzy[i] + gupzz[i]*gzzz[i]
                            + gupxy[i]*gyzx[i] + gupxz[i]*gzzx[i] + gupyz[i]*gzzy[i]
                            + gupxy[i]*gxzy[i] + gupxz[i]*gxzz[i] + gupyz[i]*gyzz[i];

                movx_Res[i] = movx_Res[i] - F2o3*Kx[i] - F8*PI*Sx[i];
                movy_Res[i] = movy_Res[i] - F2o3*Ky[i] - F8*PI*Sy[i];
                movz_Res[i] = movz_Res[i] - F2o3*Kz[i] - F8*PI*Sz[i];
            }
        }
        RHS_KERNEL_TIMER_ADD(KB_KO_CONSTRAINT, timer_ko_constraint);


            
    return 0; // success
}