 }
 
 __device__ void DAXPY3(double DA1, double DA2, double DA3, const double* DX, double* DY1, double* DY2, double* DY3) {
-	#pragma unroll 1
+	#pragma unroll 8
 	for (int i = 0; i < NSP; i++) {
 		DY1[i] += DA1 * DX[i];
 		DY2[i] += DA2 * DX[i];

    double F1[NSP];
    double F2[NSP];
    double TMP[NSP];
    #pragma unroll 1
    for (int i = 0; i < NSP; ++i) {
    	F2[i] = HrkE1 * Z1[i] + HrkE2 * Z2[i] + HrkE3 * Z3[i];
    }
    #pragma unroll 1
    for (int i = 0; i < NSP; ++i) {
    	TMP[i] = rkE[0] * F0[i] + F2[i];
    }
    dgetrs(E1, TMP, ipiv1);
    double Err = RK_ErrorNorm(scale, TMP);
    if (Err >= 1.0 && (FirstStep || Reject)) {
        #pragma unroll 1
    	for (int i = 0; i < NSP; i++) {
        	TMP[i] += Y[i];
        }
    	dydt(t, pr, TMP, F1);
    	#pragma unroll 1
    	for (int i = 0; i < NSP; i++) {
        	TMP[i] = F1[i] + F2[i];
        }

            }
            NewtonIncrementOld = NewtonIncrement;
            //! Update solution: Z4 <-- Z4 + DZ4
            #pragma unroll 8
            for (int i = 0; i < NSP; i++) {
            	Z4[i] += DZ4[i];
            }
