Ändern von Eingaben mit RcppArmadillo

Question

Ich versuche, LU-Dekomposition mit vollständigem Pivotieren unter Verwendung von RcppArmadillo zu implementieren. Glücklicherweise habe ich this Matlab Code, der tut, was ich will, aber ich habe einige Herausforderungen, die es in Armadillo umwandeln.

Ich wollte meine gecpLU Funktion arbeiten wie arma::LU wo Sie Eingang L, U und P und die arma::LU Funktion modifiziert die L, U und P-Matrizen machen, die L, U und P. anstatt Rückkehr

eingegeben werden

ich weiß, dass mit regelmäßigen Rcpp Sie Eingaben leicht wie so ändern können:

NumericVector example(NumericVector X) { 
    X = 2 * X; 
    return X; 
} 

Dies würde Rückkehr ein Vektor zweimal die Eingabe, und auch die Eingabe in gleich zwei mal den ursprünglichen Wert ändern. Ich habe jedoch schnell festgestellt, dass dies bei RcppArmadillo nicht funktioniert.

arma::colvec example(arma::colvec X) { 
    X = 2 * X; 
    return X; 
} 

Ich verstehe das nicht eine Eingabe ändern, wenn R ausgesetzt, weil die arma Objekte Kopien der R sind Objekte, so dass Sie nicht direkt auf das R-Objekt ändern, aber ich fühle mich wie ich noch sein soll Lage, eine Funktion wie Armadillo LU zu schreiben wie so: ich

#include <RcppArmadillo.h> 
using namespace Rcpp; 
// [[Rcpp::depends(RcppArmadillo)]] 
int gecpLU(arma::mat L, arma::mat U, arma::mat P, arma::mat Q, arma::mat A) { 
    // Take A and overwrite LUPQ such that A=P*L*U*Q 
    int n=A.n_rows; 
    P.eye(n,n); 
    Q.eye(n,n); 
    arma::mat AA=A; 
    // for (int i=0;i<(n-1);i++) { 
    // delete a whole bunch of stuff not relevant to question 
    // } 
    L.eye(n,n); 
    arma::mat tempmat=arma::trimatl(AA); 
    tempmat.diag()*=0; 
    L=L-tempmat; 
    U=arma::trimatu(AA); 

    return 0; 
} 

// [[Rcpp::export]] 
List test(arma::mat A) { 
    arma::mat L1,U1,P1,L2,U2,P2,Q; 
    arma::lu(L1,U1,P1,A); 
    gecpLU(L2,U2,P2,Q,A); 
    return List::create(_["L1"]=L1, 
         _["U1"]=U1, 
         _["P1"]=P1, 
         _["L2"]=L2, 
         _["U2"]=U2, 
         _["P2"]=P2, 
         _["Q"]=Q); 
} 

In diesem Fall vorbei R bin nicht Matrizen meiner gecpLU Funktion, aber arma::mat so sollte es in der Lage sein, die Eingaben zu ändern.

Wenn ich test ausführen bekomme ich Matrizen für L1, U1 und P1, aber 0x0 Matrizen für L2, U2, P2 und Q. Ich fühle mich wie ich etwas falsch verstanden werden muss. Können Eingaben mit RcppArmadillo geändert werden? Wenn nicht, was ist der beste Weg, um 4 Matrizen auszugeben? Eine Liste?

Answer 1

Under:

int gecpLU(arma::mat L, arma::mat U, arma::mat P, arma::mat Q, arma::mat A) 

Sie schaffen neue Kopien von jedem dieser Matrizen und dann am Ende der Funktion, die sie zerstört werden.

Was Sie erwarten, ist das Objekt wird geändert. Dazu müssen Sie am Ende des Objekttyps eine & anhängen, damit der Compiler die Referenz des Objekts ermitteln kann.

void gecpLU(arma::mat& L, arma::mat& U, arma::mat& P, arma::mat& Q, arma::mat& A) 

Hinweis, ich änderte auch den Rückgabetyp gecpLU von int zu void. Siehe:

#include <RcppArmadillo.h> 
using namespace Rcpp; 
// [[Rcpp::depends(RcppArmadillo)]] 
void gecpLU(arma::mat& L, arma::mat& U, arma::mat& P, arma::mat& Q, arma::mat& A) { 
    // Take A and overwrite LUPQ such that A=P*L*U*Q 
    int n=A.n_rows; 
    P.eye(n,n); 
    Q.eye(n,n); 
    arma::mat AA=A; 
    // for (int i=0;i<(n-1);i++) { 
    // delete a whole bunch of stuff not relevant to question 
    // } 
    L.eye(n,n); 
    arma::mat tempmat=arma::trimatl(AA); 
    tempmat.diag()*=0; 
    L=L-tempmat; 
    U=arma::trimatu(AA); 
} 

// [[Rcpp::export]] 
List test(arma::mat A) { 
    arma::mat L1,U1,P1,L2,U2,P2,Q; 
    arma::lu(L1,U1,P1,A); 
    gecpLU(L2,U2,P2,Q,A); 
    return List::create(_["L1"]=L1, 
         _["U1"]=U1, 
         _["P1"]=P1, 
         _["L2"]=L2, 
         _["U2"]=U2, 
         _["P2"]=P2, 
         _["Q"]=Q); 
} 

Einfach behelfsmäßigen Beispiel den Pass anhand zeigen (das ermöglicht Modifikation) vs. durch Kopie übergeben (das zerstört das Objekt am Ende)

#include <RcppArmadillo.h> 
// [[Rcpp::depends(RcppArmadillo)]] 
using namespace Rcpp; 

void reference_obj(arma::vec& y){ 
    y.fill(1); 
} 

void copy_obj(arma::vec y){ 
    y.fill(0); 
} 

// [[Rcpp::export]] 
arma::vec on_reference_mod(arma::vec x) { 

    reference_obj(x); 

    return x; 
} 


// [[Rcpp::export]] 
arma::vec on_copy_mod(arma::vec x) { 

    copy_obj(x); 

    return x; 
} 

/*** R 
# Should get a vector of 1's 
on_reference_mod(1:10) 

# Should get a vector of 1:10 
on_copy_mod(1:10) 
*/