 	register double x1 = 1.0 + rkC[0] * r;
 	register double x2 = 1.0 + rkC[1] * r;
 	register double x3 = 1.0 + rkC[2] * r;
-	#pragma unroll 1
+	#pragma unroll 8
 	for (int i = 0; i < NSP; i++) {
 		Z1[i] = CONT[i] + x1 * (CONT[NSP + i] + x1 * CONT[NSP + NSP + i]);
 		Z2[i] = CONT[i] + x2 * (CONT[NSP + i] + x2 * CONT[NSP + NSP + i]);

*/
__device__ void RK_Decomp(double H, double* E1, cuDoubleComplex* E2, const double* Jac, int* ipiv1, int* ipiv2, int* info) {
	cuDoubleComplex temp = make_cuDoubleComplex(rkAlpha/H, rkBeta/H);
	#pragma unroll 1
	for (int i = 0; i < NSP; i++)
	{
		#pragma unroll 1
		for(int j = 0; j < NSP; j++)
		{
			E1[i + j * NSP] = -Jac[i + j * NSP];

		}
	}
	else{
		#pragma unroll 8
		for (int k = 0; k < NSP; k++) {
			if (cuCabs(b[k]) > 0) {
				if (nounit) {
					b[k] = cuCdiv(b[k], A[k + k * NSP]);
				}
				#pragma unroll 8
				for (int i = k + 1; i < NSP; i++)
				{
					b[i] = cuCsub(b[i], cuCmul(b[k], A[i + k * NSP]));
