SEXP eta = PROTECT(allocVector(REALSXP,H_c)); n_prot++; double *eta_c; eta_c = REAL(eta); for (i=0;i<N_c;i++) { start_c = i * H_c; for (j=0; j<H_c; j++) { eta_c[start_c + j] = ...
It looks you expect to be able to write into the N_c*H_c element of eta, but you only allocated H_c elements for it. Bill Dunlap TIBCO Software wdunlap tibco.com On Fri, Feb 2, 2018 at 6:22 AM, Jesper Hybel Pedersen <jes...@dtu.dk> wrote: > Hi > > I'm trying to develop some C code to find the fixpoint of a contraction > mapping, the code compiles and gives the right results when executed in R. > However R-gui session is frequently terminated. I'm suspecting some access > violation error due to the exception code 0xc0000005 > In the error report windows 10 gives me. > > It is the first time I'm writing any C-code so I'm guessing I have done > something really stupid. I have been trying to debug the code for a couple > of days now, > But I simply can't figure out what generates the problem. Could it be > something particular to my windows set up and security stuff? > > > I'm in the process of reading Writing R Extensions and Hadley Wickham's > Advanced R but might have missed something. > > The windows error report: > > Faulting application name: Rgui.exe, version: 3.33.6774.0, time stamp: > 0x58bd6d26 > Faulting module name: R.dll, version: 3.33.6774.0, time stamp: 0x58bd6d0b > Exception code: 0xc0000005 > Fault offset: 0x000000000010b273 > Faulting process id: 0x1d14 > Faulting application start time: 0x01d39aede45c96e9 > Faulting application path: C:\Program Files\R\R-3.3.3\bin\x64\Rgui.exe > Faulting module path: C:\Program Files\R\R-3.3.3\bin\x64\R.dll > Report Id: c78d7c52-72c5-40f3-a3cc-927323d2af07 > Faulting package full name: > Faulting package-relative application ID: > > > ####### How I call the C-function in R > > dyn.load("C://users//jeshyb//desktop//myC//lce_fixpoint_cc.dll") > > > N = 10 > H = 3 > v <- rnorm(N*H) > mu <- 0.1 > psi <- matrix(c(1,0,0.5,0.5,0,1),nrow=2) > K <- dim(psi)[1] > p <- rep(1/H,N*H) > error <- 1e-10 > > > f<-function(p,v,mu,psi,N,H,K) > { > > .Call("lce_fixpoint_cc",p, v, mu, psi, as.integer(N), as.integer(H), > as.integer(K),error) > } > > > for (i in 1:100) > { > v <- rnorm(N*H) > p <- rep(1/H,N*H) > > a<-f(p,v,mu,psi,N,H,K) > } > > > a<-f(p,v,mu,psi,N,H,K) > plot(a) > > > > ######## The C- function > > > > #include <R.h> > #include <Rinternals.h> > > > SEXP lce_fixpoint_cc(SEXP q, SEXP v, SEXP mu, SEXP psi, SEXP N,SEXP H, > SEXP K, SEXP err) > { > > int n_prot = 0; > /* Make ready integer and double constants */ > PROTECT(N); n_prot++; > PROTECT(H); n_prot++; > PROTECT(K); n_prot++; > int N_c = asInteger(N); > int H_c = asInteger(H); > int K_c = asInteger(K); > > double mu_c = asReal(mu); > double mu2_c = 1.0 - mu_c; > double error_c = asReal(err); > double lowest_double = 1e-15; > double tmp_c; > double denom; > double error_temp; > double error_i_c; > > > /* Make ready vector froms input */ > PROTECT(q); n_prot++; > PROTECT(v); n_prot++; > PROTECT(psi); n_prot++; > double *v_c; v_c = REAL(v); > double *psi_c; psi_c = REAL(psi); > > /* Initialize new vectors not given as input */ > SEXP q_copy = PROTECT(duplicate(q)); n_prot++; > double *q_c; q_c = REAL(q_copy); > > SEXP q_new = > PROTECT(allocVector(REALSXP,length(q))); > n_prot++; > double *q_new_c; q_new_c = REAL(q_new); > > SEXP eta = PROTECT(allocVector(REALSXP,H_c)); > n_prot++; > double *eta_c; eta_c = REAL(eta); > > SEXP exp_eta = PROTECT(allocVector(REALSXP,H_c)); > n_prot++; > double *exp_eta_c; exp_eta_c = REAL(exp_eta); > > SEXP psi_ln_psiq = > PROTECT(allocVector(REALSXP,H_c)); n_prot++; > double *psi_ln_psiq_c; psi_ln_psiq_c = > REAL(psi_ln_psiq); > > int not_converged; > int maxIter = 10000; > int iter; > int start_c; > > /* loop indeces */ > int i; > int j; > int k; > > /* loop over observational units to find choice > probabilities for i=1,...,N */ > for (i=0;i<N_c;i++) > { > > start_c = i * H_c; > not_converged = 1; > iter = 0; > > while(iter < maxIter > && not_converged) > { > > /* make v_ij > + (1-mu)*ln(q_ij) */ > > for (j=0; j<H_c; j++) > > { > > eta_c[start_c + j] = v_c[start_c + j] + > mu2_c * log(q_c[start_c + j]); > > psi_ln_psiq_c[j] = 0.0; > > } > > > /* Make psi_ln_psiq_c vector for individual i */ > > for (k=0;k<K_c;k++) > > { > > tmp_c = 0.0; > > > /* Calculate row k of psi %*% q */ > > for (j=0;j<H_c;j++) > > { > > tmp_c += > psi_c[k + j*K_c] * q_c[start_c +j]; > > } > > > tmp_c = mu2_c > * log(tmp_c); > > > for (j=0;j<H_c;j++) > > { > > > psi_ln_psiq_c[j] += psi_c[k + j*K_c] * tmp_c; > > } > > } > > > denom = 0.0; > > for (j=0;j<H_c;j++) > > { > > exp_eta_c[start_c + j] = exp( > eta_c[start_c + j] - psi_ln_psiq_c[j] ) + lowest_double; > > denom += exp_eta_c[start_c + j]; > > } > > > error_i_c = 0.0; > > error_temp = 0.0; > > > /* calculate error and update choice prob */ > > for (j=0;j<H_c;j++) > > { > > q_new_c[start_c + j] = exp_eta_c[start_c > + j]/denom; > > error_temp = fabs(q_new_c[start_c + j] - > q_c[start_c + j]); > > if (error_temp>error_i_c) > > { > > > error_i_c = error_temp; > > } > > q_c[start_c + j] = q_new_c[start_c + j]; > > } > > > not_converged = error_i_c > error_c; > > iter++; > } /* End while loop > for individual i to solve for q_i */ > > > } /* End loop over individuals */ > > > UNPROTECT(n_prot); > return(q_new); > } > > > > Best regards > Jesper Hybel > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > [[alternative HTML version deleted]] ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel