//HELPER CODE #define mat_1d(i,j,mat,dim) mat[i+dim*j] __device__ MYFTYPE rhs(MYFTYPE* y, MYFTYPE* dydx, MYFTYPE* q, MYDTYPE nkinStates) { MYFTYPE k; for (MYDTYPE i = 0; i < nkinStates; i++) { k = 0; for (MYDTYPE j = 0; j < nkinStates; j++) { k += mat_1d(i, j, q, nkinStates) * y[j]; } dydx[i] = k; } } __device__ void Cubackwards_euler(MYFTYPE dt, MYDTYPE N, MYDTYPE nkinStates, MYFTYPE *y, MYFTYPE* matq, MYFTYPE *yout) { MYFTYPE top, bot, dw; MYFTYPE h = dt / 1; MYFTYPE* dydx = (MYFTYPE*)malloc(nkinStates*sizeof(MYFTYPE)); for (MYDTYPE i = 0; i < nkinStates; i++) { yout[i] = y[i]; } for (MYDTYPE j = 0; j < N; j++) { rhs(yout, dydx, matq,nkinStates); for (MYDTYPE state = 0; state < nkinStates; state++) { top = (yout[state] - y[state]) - h * dydx[state]; bot = 1 - h * mat_1d(state, state, matq, nkinStates); dw = top / bot; yout[state] = yout[state] - dw; } } free(dydx); } //END HELPER CODE ////__device__ MYDTYPE CuDerivModel_CO(MYFTYPE dt, MYFTYPE v, MYFTYPE &c1, MYFTYPE &o, MYFTYPE gbar_CO, MYFTYPE a12_CO, MYFTYPE a21_CO, MYFTYPE z12_CO, MYFTYPE z21_CO) { //DECLARES MYFTYPE yout[MOD_STATES]; MYFTYPE q[MOD_STATES*MOD_STATES]; MYFTYPE ks[MOD_STATES*MOD_STATES]; for (int i = 0; i < MOD_STATES; i++) { for (int j = 0; j < MOD_STATES; j++) { mat_1d(i, j, ks, MOD_STATES) = 0; } } //END DECLARES //FUNCTION CALLS//mat_1d(0, 1, ks, 2) = k12; mat_1d(1, 0, ks, 2) = k21; //KINBODY for (MYDTYPE i = 0; i < MOD_STATES; i++) {//This can defenitely be more efficient!! for (MYDTYPE j = 0; j < MOD_STATES; j++) { if (i != j) { mat_1d(i, j, q, MOD_STATES) = mat_1d(j, i, ks, MOD_STATES); } else { for (MYDTYPE l = 0; l < MOD_STATES; l++) { mat_1d(i, j, q, MOD_STATES) = mat_1d(i, j, q, MOD_STATES) - mat_1d(i, l, ks, MOD_STATES); } } } } Cubackwards_euler(dt, 3, MOD_STATES, y, q,yout); MYFTYPE sum=0; //conserve 1 for (MYDTYPE i = 0; i < 2; i += 1) { sum += yout[i]; } //yout = [curry/sum for curry in you] //END KINBODY //END DERIV