Hi Charles, lu_factorize factors the matrix A into a lower triangular matrix L (with unit diagonal) and an upper triangular matrix U. The values in A are overwritten with these values.
If you want to obtain the inverse, you have to call viennacl::linalg::lu_substitute(vcl_A, vcl_B); where vcl_B is the unit matrix. The inverse will be then stored in vcl_B. Best regards, Karli On 09/08/2016 07:45 PM, Charles Determan wrote: > I am trying to calculate the inverse of a matrix taking the advice from > a previous post > (https://sourceforge.net/p/viennacl/discussion/1143678/thread/ba394d35/) > suggesting the use of LU factorization. So I do the following: > > I have vcl_A matrix > > viennacl::vector<T> vcl_lu_rhs(vcl_A.size1()); > > // solution of a full system right into the load vector vcl_rhs: > viennacl::linalg::lu_factorize(vcl_A); > viennacl::linalg::lu_substitute(vcl_A, vcl_lu_rhs); > > std::cout << "matrix A" << std::endl; > std::cout << vcl_A << std::endl; > > std::cout << "vector" << std::endl; > std::cout << vcl_lu_rhs << std::endl; > > However, neither of these outputs is remotely close to the output I > expect to see for the inverse of a matrix. In R the output would be: > >># mat = vcl_A >> mat > [,1] [,2] [,3] [,4] > [1,] -1.0099356 0.19566691 0.47349181 2.2673060 > [2,] -0.7398383 0.81302435 0.34390506 0.4029221 > [3,] 1.0020811 -0.06548085 -0.09373213 0.3257177 > [4,] 1.1549178 0.87441621 1.53483119 -0.5862660 >> solve(mat) > [,1] [,2] [,3] [,4] > [1,] -0.09714492 0.004242602 0.8125165 0.07863873 > [2,] -0.45207626 1.583809306 0.8987234 -0.16052995 > [3,] 0.46074427 -0.886299611 -0.9398235 0.65059443 > [4,] 0.34057494 0.050298145 0.4806311 -0.08698459 > > but with the above viennacl code I see: > > matrix A > [4,4]((-1.00994,0.195667,0.473492,2.26731),(0.73256,0.669687,-0.00295601,-1.25802),(-0.992223,0.192126,0.376645,2.81709),(-1.14356,1.63983,5.52547,-11.4963)) > vector > [4](-0,0,0,-0) > > Did I miss something here? > > Thanks, > Charles > > > ------------------------------------------------------------------------------ > > > > _______________________________________________ > ViennaCL-devel mailing list > [email protected] > https://lists.sourceforge.net/lists/listinfo/viennacl-devel > ------------------------------------------------------------------------------ _______________________________________________ ViennaCL-devel mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/viennacl-devel
