2016-04-11 20 views
0

RcppArmadilloを使用して完全なピボットを使用してLU分解を実装しようとしています。幸いにも私はthisのMatlabコードを持っていますが、私はArmadilloに変換するいくつかの課題を抱えています。 RcppArmadilloで入力を変更する

私はあなたの入力L、U、およびP、および arma::LU関数はL、UおよびP.

を返すのではなく、入力されたL、U、およびP行列を修正 arma::LUように私の gecpLU機能の仕事を作りたかったです

は、私は定期的にRcppとあなたがそうのように簡単に入力を変更することができることを知っている:

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

これは、ベクトルを2回入力を返す、とも等しい2倍元の値に入力を変更します。しかし、私はすぐにこれがRcppArmadilloではうまくいかないことを発見しました。

私はRにさらされたときに、あなたが直接Rオブジェクトを変更することはできませんので armaオブジェクトはRのオブジェクトのコピーであるため、これは、入力を変更しません理解し、私はまだする必要がありますように私は感じ
arma::colvec example(arma::colvec X) { 
    X = 2 * X; 
    return X; 
} 

そのようアルマジロのLUのような関数を書くことができる:

#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); 
} 

この場合、私は私のgecpLU関数にR行列を渡していないのですが、arma::matので、入力を変更することができるはずです。

私がtestを実行すると、L1、U1、P1の行列が得られますが、L2、U2、P2、Qの行列は0x0です。私は何かを誤解しているように感じます。 RcppArmadilloで入力を変更することは可能ですか?そうでない場合、4つの行列を出力する最良の方法は何ですか?リスト?アンダー

+2

最も簡単な方法は、中に転送することである(すなわち、末尾にオブジェクトを破棄)コピーを通過Rcppベクタを生成し、補助コンストラクタを使用して 'arma'オブジェクトを作成します(詳細はhttp://arma.sourceforge.net/docs.html#Colの' Advanced constructors'セクションを参照してください)。 –

答えて

2

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

は、あなたは、彼らが破壊されている関数の最後にそれらの行列の各の新しいコピーを作成しています。

オブジェクトは変更されています。これを行うには、オブジェクト型の最後に&を追加して、コンパイラがオブジェクトの参照を取得できるようにする必要があります。

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

注、私はまた、voidintからgecpLUの戻り値の型を変更しました。参照:

#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); 
} 

簡易仮設例(つまり変形可能)参照によってパスを実証するために対

#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) 
*/ 
+0

ありがとう!私は何かが足りないことを知っていた。 – Carl

関連する問題