R devel - Dec 2006 - FW: Passing lists from R to C, extracting elements, and sending lists back again

Thank you to everyone who responded to my previous post regarding the
integration into R of C programs that use external libraries.  I have another
issue that I simply have not been able to figure out using documentation, list
archives, and so forth.

I have data that is stored in R lists.  I would like to pass the list to C code
using the .Call function, perform mathematical operations on the elements of the
lists (which might be other lists), and return results in a new list.

As a simple example, suppose I have a list of n matrices (dimension k-by-k) that
I want to invert, and return the result in another list.  I wrote the C code
below that reads the list as a SEXP.  Then, for each of the i=1..n matrices, I
convert the R vector to a gsl matrix, perform the linear algebra operations and
populate the new list with the inverted matrix.

The problem is that I can't figure out the best way to get the real elements
into the list SEXP to return it to R.  When I run the code below, the compiler
(R CMD SHLIB sexp1.c) gives me:

-bash-3.00$ R CMD SHLIB sexp1.c
gcc -I/usr/lib64/R/include -I/usr/lib64/R/include -I/usr/include
-I/usr/local/include   -fpic  -O2 -g -c sexp1.c -o sexp1.o
sexp1.c: In function `invMatList':
sexp1.c:37: warning: passing arg 3 of `SET_VECTOR_ELT' from incompatible
pointer type gcc -shared -L/usr/local/lib64 -o sexp1.so sexp1.o -L/usr/lib64
-lgsl -lgslcblas -L/usr/lib  -L/usr/lib64/R/lib -lR

So, not only would it be helpful to understand how to solve this specific
problem, but I would greatly appreciate some general guidance on how to best
more data from lists or C datatypes and back again (especially arrays).  Also,
are the functions and macros described in the Rinternals header documented
anywhere?  I can't seem to find the source code, or a description of what
these functions actually do.

Thanks again for any help you can provide.

Best wishes,

Michael Braun
MIT Sloan School of Management
braunm at mit.edu

	SEXP r = PROTECT(allocVector(VECSXP,n));
	double rm[k][k];

	for (i=0; i<n; i=i+1) {

		for (a=0; a<k; ++a){ // convert to gsl matrix
			for (b=0; b<k; ++b){
				ind = b*k + a;
				m = REAL(VECTOR_ELT(list,i))[ind];
				gsl_matrix_set(mat, b, a, m);
			}
		} // end conversion loops

 		gsl_linalg_LU_decomp(mat, p, &s);
		gsl_linalg_LU_invert(mat, p, inv);


		for (a=0; a<k; ++a){ // convert from gsl matrix
			for (b=0; b<k; ++b){
				rm[b][a] = gsl_matrix_get(inv,b,a);
			}
		}
		SET_VECTOR_ELT(r, i, rm);
	}
	UNPROTECT(1);
	return (r);
}

>
> I have data that is stored in R lists.  I would like to pass the list to 
> C code using the .Call function, perform mathematical operations on the 
> elements of the lists (which might be other lists), and return results 
> in a new list.
>
> As a simple example, suppose I have a list of n matrices (dimension 
> k-by-k) that I want to invert, and return the result in another list. 
> I wrote the C code below that reads the list as a SEXP.  Then, for each 
> of the i=1..n matrices, I convert the R vector to a gsl matrix, perform 
> the linear algebra operations and populate the new list with the 
> inverted matrix.
>
> The problem is that I can't figure out the best way to get the real 
> elements into the list SEXP to return it to R.  When I run the code 
> below, the compiler (R CMD SHLIB sexp1.c) gives me:
>
> -bash-3.00$ R CMD SHLIB sexp1.c
> gcc -I/usr/lib64/R/include -I/usr/lib64/R/include -I/usr/include
-I/usr/local/include   -fpic  -O2 -g -c sexp1.c -o sexp1.o
> sexp1.c: In function `invMatList':
> sexp1.c:37: warning: passing arg 3 of `SET_VECTOR_ELT' from
incompatible pointer type gcc -shared -L/usr/local/lib64 -o sexp1.so sexp1.o
-L/usr/lib64 -lgsl -lgslcblas -L/usr/lib  -L/usr/lib64/R/lib -lR
>
<snip>
> 	SEXP r = PROTECT(allocVector(VECSXP,n));
> 	double rm[k][k];
This allocates memory for one k x k C array wherever C allocates temporary 
storage. You need to allocate n vectors of length k*k on the R heap to 
create a list of n matrices in R.   So you need to use allocVector(k*k, 
REALSXP) to allocate space for each element of the list, inside the loop, 
and then copy the data into that rather than into rm.

SET_VECTOR_ELT is giving a warning about incompatible pointer types, but 
the more important issue is that you need to return pointers to n 
different blocks of memory that are not going to go away after your 
function finishes.

 	-thomas
>
>
> 	for (i=0; i<n; i=i+1) {
>
> 		for (a=0; a<k; ++a){ // convert to gsl matrix
> 			for (b=0; b<k; ++b){
> 				ind = b*k + a;
> 				m = REAL(VECTOR_ELT(list,i))[ind];
> 				gsl_matrix_set(mat, b, a, m);
> 			}
> 		} // end conversion loops
>
> 		gsl_linalg_LU_decomp(mat, p, &s);
> 		gsl_linalg_LU_invert(mat, p, inv);
>
>
> 		for (a=0; a<k; ++a){ // convert from gsl matrix
> 			for (b=0; b<k; ++b){
> 				rm[b][a] = gsl_matrix_get(inv,b,a);
> 			}
> 		}
> 		SET_VECTOR_ELT(r, i, rm);
> 	}
> 	UNPROTECT(1);
> 	return (r);
> }
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>
Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle

Dear Michael,

I had to tackle a similar issue while programming the AMORE package.
I recommend you to have a look at the source code of its "copynet.c"
file. You
may find useful the reading of the "copynet_RC" function ( Copies the
SEXP net
to the *ptnet) and the "copynet_CR" function (Copies *ptnet to SEXP
net).

Of course, I do not pretend to establish what best practices should be, just to
show you how I solved a similar issue following my personal understanding of
what was appropriated in that situation.

Best wishes and Merry Christmas!
Manuel