| 12 | |
| 13 | |
| 14 | template<typename MatrixType> void triangular_square(const MatrixType& m) |
| 15 | { |
| 16 | typedef typename MatrixType::Scalar Scalar; |
| 17 | typedef typename NumTraits<Scalar>::Real RealScalar; |
| 18 | typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; |
| 19 | |
| 20 | RealScalar largerEps = 10*test_precision<RealScalar>(); |
| 21 | |
| 22 | typename MatrixType::Index rows = m.rows(); |
| 23 | typename MatrixType::Index cols = m.cols(); |
| 24 | |
| 25 | MatrixType m1 = MatrixType::Random(rows, cols), |
| 26 | m2 = MatrixType::Random(rows, cols), |
| 27 | m3(rows, cols), |
| 28 | m4(rows, cols), |
| 29 | r1(rows, cols), |
| 30 | r2(rows, cols); |
| 31 | VectorType v2 = VectorType::Random(rows); |
| 32 | |
| 33 | MatrixType m1up = m1.template triangularView<Upper>(); |
| 34 | MatrixType m2up = m2.template triangularView<Upper>(); |
| 35 | |
| 36 | if (rows*cols>1) |
| 37 | { |
| 38 | VERIFY(m1up.isUpperTriangular()); |
| 39 | VERIFY(m2up.transpose().isLowerTriangular()); |
| 40 | VERIFY(!m2.isLowerTriangular()); |
| 41 | } |
| 42 | |
| 43 | // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); |
| 44 | |
| 45 | // test overloaded operator+= |
| 46 | r1.setZero(); |
| 47 | r2.setZero(); |
| 48 | r1.template triangularView<Upper>() += m1; |
| 49 | r2 += m1up; |
| 50 | VERIFY_IS_APPROX(r1,r2); |
| 51 | |
| 52 | // test overloaded operator= |
| 53 | m1.setZero(); |
| 54 | m1.template triangularView<Upper>() = m2.transpose() + m2; |
| 55 | m3 = m2.transpose() + m2; |
| 56 | VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1); |
| 57 | |
| 58 | // test overloaded operator= |
| 59 | m1.setZero(); |
| 60 | m1.template triangularView<Lower>() = m2.transpose() + m2; |
| 61 | VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1); |
| 62 | |
| 63 | VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(), |
| 64 | m3.conjugate().template triangularView<Lower>().toDenseMatrix()); |
| 65 | |
| 66 | m1 = MatrixType::Random(rows, cols); |
| 67 | for (int i=0; i<rows; ++i) |
| 68 | while (numext::abs2(m1(i,i))<RealScalar(1e-1)) m1(i,i) = internal::random<Scalar>(); |
| 69 | |
| 70 | Transpose<MatrixType> trm4(m4); |
| 71 | // test back and forward subsitution with a vector as the rhs |