[Rcpp-devel] Matrix assignment with inline

[Kaveh Vakili]

Hi list,

I'm trying to learn inline by translating snippets of c++ code from 
the book "Numerical Recipes".

in page 44-45, there is a c++ code to perform Gauss-Jordan elimination of a matrix (see below). I understand that they have a large template of additional operator/function/class (it is online: www.nr.com/codefile.php?nr3 ) which means the codes in the book have to be slightly 'translated' (correct me if i'm wrong) unto proper inline inputs. I think i have a problem with translating the matrix functions. I have searched the devel archives, but have only found codes to assign full row/column (not individual matrix entries). 

Anyhow, my use of 'a[j][k]' in line 44 (and similar threats later) seem to be an issue: none of the examples in Rcpp seems to use it. My question is how do you translate this construct (i.e. a[j][k]) so that inline can 'eat' it ?


library(Rcpp)
library(inline)


fun3<-'
	NumericMatrix a(x) ;
	NumericMatrix b = clone<NumericMatrix>(a);
	int i,icol,irow,j,k,l,ll,n=a.nrow(),m=a.ncol();
	std::vector<int> indxc(n); 
	std::vector<int> indxr(n);
	std::vector<int> ipiv(n);
	double big,dum,piinv;
	for(j=0;j<n;j++) ipiv[j]=0;
	for(i=0;i<n;i++ ){
		big=0.0;
		for(j=0;j<n;j++){
			if(ipiv[j] != 1)
				for(k=0;k<n;k++){
					if(ipiv[k] == 0) {
						if(abs(a[j][k]) >= big) {
							big = abs(a[j][k]);
							irow = j;
							icol = k;
						}
					}
				}
			++(ipiv[icol]);
			if(irow != icol){
				for (l=0;l<n;l++){
					dum = a[irow][l];
					a[irow][l] = a[icol][l];
					a[icol][l] = dum;
				}
				for (l=0; l<m ; l++){
					dum = b[irow][l];
					b[irow][l] = b[icol][l];
					b[icol][l] = dum;
				}
			}
			idxr[i] = irow;
			idxc[i] = icol;
			if(a[icol][icol]==0.0) throw ("gaussj: Singular Matrix");
			pivinv=1.0/a[icol][icol];
			a[icol][icol]=1.0;
			for(l=0;l<n;l++) a[icol][l] *= pivinv;
			for(l=0;l<m;l++) b[icol][l] *= pivinv;
			for(ll=0;ll<n,ll++)
				if(ll != icol) {
					dum = a[ll][icol];
					a[ll][icol] = 0.0;
					for(l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
					for(l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
				}
		}
		for(l=n-1;l>=0;l--){
			if(indxr[l] != indxc[l])
				for(k=0;k<n;k++)
					dum = a[k][indxr[l]];
					a[k][indxr[l]] = a[k][indxc[l]];
					a[k][indxc[l]] = dum;
		}
}
'
f.Rcpp<-cxxfunction(signature(x="matrix"),fun3,plugin="Rcpp",verbose=T)