 		sc[i] = 1.0 / (ATOL + fabs(y0[i]) * RTOL);
 	}
 }
-
-#ifdef DEBUG
-#include <stdio.h>
-#define PRINT(X, Y)  \
-	if (threadIdx.x == 0) {\
-		printf("Block:%d\t", blockIdx.x); \
-		printf(X, Y); \
-	} \
-	else \
-		(void)0
-#else
-#define PRINT(X, Y)
-#endif
 
 #ifdef NEWTON_UNROLL
 	#define NEWTON_UNROLLER	#pragma unroll NEWTON_UNROLL


            NewtonDone = (NewtonRate * NewtonIncrement <= NewtonTol);
#ifndef NEWTON_UNROLL
            PRINT("Done:%d\n", ALL(NewtonDone));
            if (ALL (NewtonDone)) break;
#else //only break if it's at the end of the unroll
            if(ALL (NewtonDone && (NewtonIter + 1) % NewtonIter == 0)) break;
#endif
            if (NewtonIter >= NewtonMaxit)
            {

			if (NewtonIter > 0) 
			{
				Theta = NewtonIncrement / NewtonIncrementOld;
				if(Theta >= 0.99) //! Non-convergence of Newton: Theta too large
					break;
				else
					NewtonRate = Theta / (ONE - Theta);
				//Predict error at the end of Newton process 
				double NewtonPredictedErr = (NewtonIncrement * pow(Theta, (NewtonMaxit - NewtonIter - 1))) / (ONE - Theta);
				if(NewtonPredictedErr >= NewtonTol) {
					//Non-convergence of Newton: predicted error too large
					double Qnewton = fmin(10.0, NewtonPredictedErr / NewtonTol);
                    Fac = 0.8 * pow(Qnewton, -ONE/((double)(NewtonMaxit-NewtonIter)));
