[Rcpp-devel] Matrix assignment with inline

[Romain Francois]

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
> Rcpp-devel at lists.r-forge.r-project.org
> 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