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)
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
