Le 03/11/10 12:16, Kaveh Vakili a écrit :
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 ?
Hi,
You can use a(j,k)
There are a lot of places in this code where Rcpp sugar might help. I
don't want to go through all of it right now, but for example this line:
for(l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
could be replaced by something like this (assuming Rcpp 0.8.8 or later):
a(_,ll) = a(_,ll) - a(_,icol) * dum ;
You can find more examples of using sugar in this file:
> system.file( "unitTests/runit.sugar.R", package = "Rcpp" )
Also, I'm not sure that this code is valid:
throw ("gaussj: Singular Matrix");
Maybe you mean :
throw std::range_error("gaussj: Singular Matrix");
or perhaps using some other exception class instead of range_error.
HTH,
Romain
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
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/czHPM7 : Rcpp Google tech talk on youtube
|- http://bit.ly/9P0eF9 : Google slides
`- http://bit.ly/cVPjPe : Chicago R Meetup slides
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
