merge with main repository

18 files changed, 52 insertions, 68 deletions

diff --git a/doc/C05_TutorialLinearAlgebra.dox b/doc/C05_TutorialLinearAlgebra.dox
index fbf809d58..c50e9c6bc 100644
--- a/doc/C05_TutorialLinearAlgebra.dox
+++ b/doc/C05_TutorialLinearAlgebra.dox
@@ -55,8 +55,8 @@ matrix with a vector or another matrix: \f$ A^{-1} \mathbf{v} \f$ or \f$ A^{-1}
 This is a general-purpose algorithm which performs well in most cases (provided the matrix \f$ A \f$
 is invertible), so if you are unsure about which algorithm to pick, choose this. The method proceeds
 in two steps. First, the %LU decomposition with partial pivoting is computed using the
-MatrixBase::partialLu() function. This yields an object of the class PartialLU. Then, the
-PartialLU::solve() method is called to compute a solution.
+MatrixBase::partialPivLu() function. This yields an object of the class PartialPivLU. Then, the
+PartialPivLU::solve() method is called to compute a solution.
 
 As an example, suppose we want to solve the following system of linear equations:
 
@@ -69,9 +69,9 @@ As an example, suppose we want to solve the following system of linear equations
 The following program solves this system:
 
 <table class="tutorial_code"><tr><td>
-\include Tutorial_PartialLU_solve.cpp
+\include Tutorial_PartialPivLU_solve.cpp
 </td><td>
-output: \include Tutorial_PartialLU_solve.out
+output: \include Tutorial_PartialPivLU_solve.out
 </td></tr></table>
 
 There are many situations in which we want to solve the same system of equations with different
@@ -91,7 +91,7 @@ problem, and whether you want to solve it at all, after you solved the first pro
 case, it's best to save the %LU decomposition and reuse it to solve the second problem. This is
 worth the effort because computing the %LU decomposition is much more expensive than using it to
 solve the equation. Here is some code to illustrate the procedure. It uses the constructor
-PartialLU::PartialLU(const MatrixType&) to compute the %LU decomposition.
+PartialPivLU::PartialPivLU(const MatrixType&) to compute the %LU decomposition.
 
 <table class="tutorial_code"><tr><td>
 \include Tutorial_solve_reuse_decomposition.cpp
@@ -102,7 +102,7 @@ output: \include Tutorial_solve_reuse_decomposition.out
 \b Warning: All this code presumes that the matrix \f$ A \f$ is invertible, so that the system 
 \f$ A \mathbf{x} = \mathbf{b} \f$ has a unique solution. If the matrix \f$ A \f$ is not invertible,
 then the system \f$ A \mathbf{x} = \mathbf{b} \f$ has either zero or infinitely many solutions. In
-both cases, PartialLU::solve() will give nonsense results. For example, suppose that we want to
+both cases, PartialPivLU::solve() will give nonsense results. For example, suppose that we want to
 solve the same system as above, but with the 10 in the last equation replaced by 9. Then the system
 of equations is inconsistent: adding the first and the third equation gives \f$ 8x + 10y + 12z = 7 \f$,
 which implies \f$ 4x + 5y + 6z = 3\frac12 \f$, in contradiction with the second equation. If we try
@@ -114,10 +114,10 @@ to solve this inconsistent system with Eigen, we find:
 output: \include Tutorial_solve_singular.out
 </td></tr></table>
 
-The %LU decomposition with \b full pivoting (class LU) and the singular value decomposition (class
+The %LU decomposition with \b full pivoting (class FullPivLU) and the singular value decomposition (class
 SVD) may be helpful in this case, as explained in the section \ref TutorialAdvSolvers_Misc below.
 
-\sa LU_Module, MatrixBase::partialLu(), PartialLU::solve(), class PartialLU.
+\sa LU_Module, MatrixBase::partialPivLu(), PartialPivLU::solve(), class PartialPivLU.
 
 
 \subsection TutorialAdvSolvers_Cholesky Cholesky decomposition
@@ -228,7 +228,7 @@ Note that the function inverse() is defined in the \ref LU_Module.
 
 Finally, Eigen also offer solvers based on a singular value decomposition (%SVD) or the %LU
 decomposition with full pivoting. These have the same API as the solvers based on the %LU
-decomposition with partial pivoting (PartialLU). 
+decomposition with partial pivoting (PartialPivLU). 
 
 The solver based on the %SVD uses the class SVD. It can handle singular matrices. Here is an example
 of its use:
@@ -245,7 +245,7 @@ svdOfA.solve(b, &x);
 \endcode
 
 %LU decomposition with full pivoting has better numerical stability than %LU decomposition with
-partial pivoting. It is defined in the class LU. The solver can also handle singular matrices.
+partial pivoting. It is defined in the class FullPivLU. The solver can also handle singular matrices.
 
 \code
 #include <Eigen/LU>
@@ -254,13 +254,13 @@ MatrixXf A = MatrixXf::Random(20,20);
 VectorXf b = VectorXf::Random(20);
 VectorXf x;
 A.lu().solve(b, &x);
-LU<MatrixXf> luOfA(A);
+FullPivLU<MatrixXf> luOfA(A);
 luOfA.solve(b, &x);
 \endcode
 
 See the section \ref TutorialAdvLU below.
 
-\sa class SVD, SVD::solve(), SVD_Module, class LU, LU::solve(), LU_Module.
+\sa class SVD, SVD::solve(), SVD_Module, class FullPivLU, LU::solve(), LU_Module.
 
 
 
@@ -281,7 +281,7 @@ Alternatively, you can construct a named LU decomposition, which allows you to r
 \code
 #include <Eigen/LU>
 MatrixXf A = MatrixXf::Random(20,20);
-Eigen::LU<MatrixXf> lu(A);
+Eigen::FullPivLU<MatrixXf> lu(A);
 cout << "The rank of A is" << lu.rank() << endl;
 if(lu.isInvertible()) {
   cout << "A is invertible, its inverse is:" << endl << lu.inverse() << endl;
@@ -292,7 +292,7 @@ else {
 }
 \endcode
 
-\sa LU_Module, LU::solve(), class LU
+\sa LU_Module, LU::solve(), class FullPivLU
 
 <a href="#" class="top">top</a>\section TutorialAdvCholesky Cholesky
 todo
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 5b055ed11..412831e14 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -215,7 +215,10 @@ ALIASES                = "only_for_vectors=This is only for vectors (either row-
                          "eigenvalues_module=This is defined in the %Eigenvalues module. \code #include <Eigen/Eigenvalues> \endcode" \
                          "label=\bug" \
                          "redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \
-                         "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\""
+                         "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \
+                         "note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \
+                         "note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \
+                         "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values."
 
 # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
 # sources only. Doxygen will then generate output that is more tailored for C.
diff --git a/doc/examples/Tutorial_PartialLU_solve.cpp b/doc/examples/Tutorial_PartialLU_solve.cpp
index 80c393f9a..98f9d6dac 100644
--- a/doc/examples/Tutorial_PartialLU_solve.cpp
+++ b/doc/examples/Tutorial_PartialLU_solve.cpp
@@ -12,7 +12,6 @@ int main(int, char *[])
    b << 3, 3, 4;
    cout << "Here is the matrix A:" << endl << A << endl;
    cout << "Here is the vector b:" << endl << b << endl;
-   Vector3f x;
-   A.partialLu().solve(b, &x);
+   Vector3f x = A.lu().solve(b);
    cout << "The solution is:" << endl << x << endl;
 }
diff --git a/doc/snippets/LU_image.cpp b/doc/snippets/FullPivLU_image.cpp
index 0d1088a2f..817bc1e2d 100644
--- a/doc/snippets/LU_image.cpp
+++ b/doc/snippets/FullPivLU_image.cpp
@@ -6,4 +6,4 @@ cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Notice that the middle column is the sum of the two others, so the "
      << "columns are linearly dependent." << endl;
 cout << "Here is a matrix whose columns have the same span but are linearly independent:"
-     << endl << m.lu().image() << endl;
+     << endl << m.fullPivLu().image(m) << endl;
diff --git a/doc/snippets/LU_kernel.cpp b/doc/snippets/FullPivLU_kernel.cpp
index e01186d38..7086e01e2 100644
--- a/doc/snippets/LU_kernel.cpp
+++ b/doc/snippets/FullPivLU_kernel.cpp
@@ -1,6 +1,6 @@
 MatrixXf m = MatrixXf::Random(3,5);
 cout << "Here is the matrix m:" << endl << m << endl;
-MatrixXf ker = m.lu().kernel();
+MatrixXf ker = m.fullPivLu().kernel();
 cout << "Here is a matrix whose columns form a basis of the kernel of m:"
      << endl << ker << endl;
 cout << "By definition of the kernel, m*ker is zero:"
diff --git a/doc/snippets/LU_solve.cpp b/doc/snippets/FullPivLU_solve.cpp
index 7323338c3..c1f88235e 100644
--- a/doc/snippets/LU_solve.cpp
+++ b/doc/snippets/FullPivLU_solve.cpp
@@ -1,15 +1,11 @@
-typedef Matrix<float,2,3> Matrix2x3;
-typedef Matrix<float,3,2> Matrix3x2;
-Matrix2x3 m = Matrix2x3::Random();
+Matrix<float,2,3> m = Matrix<float,2,3>::Random();
 Matrix2f y = Matrix2f::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is the matrix y:" << endl << y << endl;
-Matrix3x2 x;
-if(m.lu().solve(y, &x))
+Matrix<float,3,2> x = m.fullPivLu().solve(y);
+if((m*x).isApprox(y))
 {
-  assert(y.isApprox(m*x));
   cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
 }
 else
   cout << "The equation mx=y does not have any solution." << endl;
-
diff --git a/doc/snippets/HouseholderQR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp
index 429bd81e3..8cce6ce6c 100644
--- a/doc/snippets/HouseholderQR_solve.cpp
+++ b/doc/snippets/HouseholderQR_solve.cpp
@@ -4,6 +4,6 @@ Matrix3f y = Matrix3f::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "Here is the matrix y:" << endl << y << endl;
 Matrix3f x;
-m.householderQr().solve(y, &x);
+x = m.householderQr().solve(y);
 assert(y.isApprox(m*x));
 cout << "Here is a solution x to the equation mx=y:" << endl << x << endl;
diff --git a/doc/snippets/LLT_solve.cpp b/doc/snippets/LLT_solve.cpp
index 76ab09ec5..7095d2cc3 100644
--- a/doc/snippets/LLT_solve.cpp
+++ b/doc/snippets/LLT_solve.cpp
@@ -3,6 +3,6 @@ typedef Matrix<float,Dynamic,2> DataMatrix;
 DataMatrix samples = DataMatrix::Random(12,2);
 VectorXf elevations = 2*samples.col(0) + 3*samples.col(1) + VectorXf::Random(12)*0.1;
 // and let's solve samples * [x y]^T = elevations in least square sense:
-Matrix<float,2,1> xy;
-(samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations), &xy);
+Matrix<float,2,1> xy
+ = (samples.adjoint() * samples).llt().solve((samples.adjoint()*elevations));
 cout << xy << endl;
diff --git a/doc/snippets/LU_computeImage.cpp b/doc/snippets/LU_computeImage.cpp
deleted file mode 100644
index 5c812cc4c..000000000
--- a/doc/snippets/LU_computeImage.cpp
+++ /dev/null
@@ -1,13 +0,0 @@
-MatrixXd m(3,3);
-m << 1,1,0,
-     1,3,2,
-     0,1,1;
-cout << "Here is the matrix m:" << endl << m << endl;
-LU<Matrix3d> lu(m);
-// allocate the matrix img with the correct size to avoid reallocation
-MatrixXd img(m.rows(), lu.rank());
-lu.computeImage(&img);
-cout << "Notice that the middle column is the sum of the two others, so the "
-     << "columns are linearly dependent." << endl;
-cout << "Here is a matrix whose columns have the same span but are linearly independent:"
-     << endl << img << endl;
diff --git a/doc/snippets/LU_computeKernel.cpp b/doc/snippets/LU_computeKernel.cpp
deleted file mode 100644
index b08f7f1ea..000000000
--- a/doc/snippets/LU_computeKernel.cpp
+++ /dev/null
@@ -1,10 +0,0 @@
-MatrixXf m = MatrixXf::Random(3,5);
-cout << "Here is the matrix m:" << endl << m << endl;
-LU<MatrixXf> lu(m);
-// allocate the matrix ker with the correct size to avoid reallocation
-MatrixXf ker(m.rows(), lu.dimensionOfKernel());
-lu.computeKernel(&ker);
-cout << "Here is a matrix whose columns form a basis of the kernel of m:"
-     << endl << ker << endl;
-cout << "By definition of the kernel, m*ker is zero:"
-     << endl << m*ker << endl;
diff --git a/doc/snippets/MatrixBase_computeInverse.cpp b/doc/snippets/MatrixBase_computeInverse.cpp
deleted file mode 100644
index ba7377a4a..000000000
--- a/doc/snippets/MatrixBase_computeInverse.cpp
+++ /dev/null
@@ -1,5 +0,0 @@
-Matrix3d m = Matrix3d::Random();
-cout << "Here is the matrix m:" << endl << m << endl;
-Matrix3d inv;
-m.computeInverse(&inv);
-cout << "Its inverse is:" << endl << inv << endl;
diff --git a/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp
new file mode 100644
index 000000000..a7b084fd0
--- /dev/null
+++ b/doc/snippets/MatrixBase_computeInverseAndDetWithCheck.cpp
@@ -0,0 +1,13 @@
+Matrix3d m = Matrix3d::Random();
+cout << "Here is the matrix m:" << endl << m << endl;
+Matrix3d inverse;
+bool invertible;
+double determinant;
+m.computeInverseAndDetWithCheck(inverse,determinant,invertible);
+cout << "Its determinant is " << determinant << endl;
+if(invertible) {
+  cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
+}
+else {
+  cout << "It is not invertible." << endl;
+}
diff --git a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
index 19e24c90b..873a9f870 100644
--- a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
+++ b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp
@@ -1,8 +1,10 @@
 Matrix3d m = Matrix3d::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
-Matrix3d inv;
-if(m.computeInverseWithCheck(&inv)) {
-  cout << "It is invertible, and its inverse is:" << endl << inv << endl;
+Matrix3d inverse;
+bool invertible;
+m.computeInverseWithCheck(inverse,invertible);
+if(invertible) {
+  cout << "It is invertible, and its inverse is:" << endl << inverse << endl;
 }
 else {
   cout << "It is not invertible." << endl;
diff --git a/doc/snippets/PartialLU_solve.cpp b/doc/snippets/PartialLU_solve.cpp
index 12441437f..fa3570ab8 100644
--- a/doc/snippets/PartialLU_solve.cpp
+++ b/doc/snippets/PartialLU_solve.cpp
@@ -2,7 +2,6 @@ MatrixXd A = MatrixXd::Random(3,3);
 MatrixXd B = MatrixXd::Random(3,2);
 cout << "Here is the invertible matrix A:" << endl << A << endl;
 cout << "Here is the matrix B:" << endl << B << endl;
-MatrixXd X;
-if(A.lu().solve(B, &X))
+MatrixXd X = A.lu().solve(B);
 cout << "Here is the (unique) solution X to the equation AX=B:" << endl << X << endl;
 cout << "Relative error: " << (A*X-B).norm() / B.norm() << endl;
diff --git a/doc/snippets/Tutorial_solve_multiple_rhs.cpp b/doc/snippets/Tutorial_solve_multiple_rhs.cpp
index fbb15165a..72dab9075 100644
--- a/doc/snippets/Tutorial_solve_multiple_rhs.cpp
+++ b/doc/snippets/Tutorial_solve_multiple_rhs.cpp
@@ -3,7 +3,7 @@ A << 1,2,3,  4,5,6,  7,8,10;
 Matrix<float,3,2> B;
 B << 3,1, 3,1, 4,1;
 Matrix<float,3,2> X;
-A.partialLu().solve(B, &X);
+X = A.lu().solve(B);
 cout << "The solution with right-hand side (3,3,4) is:" << endl;
 cout << X.col(0) << endl;
 cout << "The solution with right-hand side (1,1,1) is:" << endl;
diff --git a/doc/snippets/Tutorial_solve_reuse_decomposition.cpp b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp
index b4112adc4..3ca06453a 100644
--- a/doc/snippets/Tutorial_solve_reuse_decomposition.cpp
+++ b/doc/snippets/Tutorial_solve_reuse_decomposition.cpp
@@ -1,13 +1,13 @@
 Matrix3f A(3,3);
 A << 1,2,3,  4,5,6,  7,8,10;
-PartialLU<Matrix3f> luOfA(A); // compute LU decomposition of A
+PartialPivLU<Matrix3f> luOfA(A); // compute LU decomposition of A
 Vector3f b;
 b << 3,3,4;
 Vector3f x;
-luOfA.solve(b, &x);
+x = luOfA.solve(b);
 cout << "The solution with right-hand side (3,3,4) is:" << endl;
 cout << x << endl;
 b << 1,1,1;
-luOfA.solve(b, &x);
+x = luOfA.solve(b);
 cout << "The solution with right-hand side (1,1,1) is:" << endl;
 cout << x << endl;
diff --git a/doc/snippets/Tutorial_solve_singular.cpp b/doc/snippets/Tutorial_solve_singular.cpp
index da94ad445..abff1ef73 100644
--- a/doc/snippets/Tutorial_solve_singular.cpp
+++ b/doc/snippets/Tutorial_solve_singular.cpp
@@ -5,5 +5,5 @@ b << 3, 3, 4;
 cout << "Here is the matrix A:" << endl << A << endl;
 cout << "Here is the vector b:" << endl << b << endl;
 Vector3f x;
-A.partialLu().solve(b, &x);
+x = A.lu().solve(b);
 cout << "The solution is:" << endl << x << endl;
diff --git a/doc/snippets/class_LU.cpp b/doc/snippets/class_FullPivLU.cpp
index 9958368f1..40d76e8e6 100644
--- a/doc/snippets/class_LU.cpp
+++ b/doc/snippets/class_FullPivLU.cpp
@@ -2,7 +2,7 @@ typedef Matrix<double, 5, 3> Matrix5x3;
 typedef Matrix<double, 5, 5> Matrix5x5;
 Matrix5x3 m = Matrix5x3::Random();
 cout << "Here is the matrix m:" << endl << m << endl;
-Eigen::LU<Matrix5x3> lu(m);
+Eigen::FullPivLU<Matrix5x3> lu(m);
 cout << "Here is, up to permutations, its LU decomposition matrix:"
      << endl << lu.matrixLU() << endl;
 cout << "Here is the L part:" << endl;