<div dir="ltr"><div>Dear all,<br><br></div>I need to solve a linear system using only the lower triangular matrix 
from the Cholesky decomposition of a SPD matrix.  My goal is to compute 
C^{-1/2} V^{-1/2} where V^{-1/2} and C are known matrices.  <br>For 
general matrix class MatrixXd I can do using the code below. But, now I 
would like to do the same using the SparseMatrix class. I tried to adapt
 my code, but it did not work. Any help will be welcome.<div><div><br><br>src <- '<br>using namespace Rcpp;<br>using namespace Eigen;<br>const Map<MatrixXd> Omega(as<Map<MatrixXd> >(Omegas));<br>const Map<MatrixXd> invSqrtV(as<Map<MatrixXd> >(invSqrtVs));<br>MatrixXd LT = Omega.llt().matrixL().solve(invSqrtV);<br>return wrap(LT);'<br><br><br>srcsparse <- '<br>using Eigen::Map; <br>using Eigen::SparseMatrix;<br>using Eigen::LLT;<br><br>const SparseMatrix<double> Omega(as<SparseMatrix<double> >(Omegas));<br>const SparseMatrix<double> invSqrtV(as<SparseMatrix<double> >(invSqrtVs));<br>Eigen::SimplicialLLT<SparseMatrix<double>, Eigen::Lower> solver;<br>SparseMatrix<double>  LT = solver.compute(Omega).matrixL()<b>.solve(invSqrtV);</b><br>return wrap(LT);'<br><br></div><div>The problem is that there is no .solve() method associated with the matrixL(). <br></div><div>It is just a intermediary step when computing<br><br>solver.compute(Omega).solve(invSqrtV)<br><br></div><div>This code computing C^{-1}V^{-1/2}, but I want only C^{-1/2}V^{-1/2}. It will be used on the Generalized Kronecker product in the context of multivariate spatial models.<br><br></div><div>Thank you.<br></div><div>-- <br><div class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>Wagner Hugo Bonat<br>----------------------------------------------------------------------------------------------<br>Department of Mathematics and Computer Science (IMADA)<br>University of Southern Denmark (SDU) and<br>Laboratório de Estatística e Geoinformação (LEG)<br>Universidade Federal do Paraná (UFPR)<br><br></div></div></div></div></div>
</div></div></div>