From cfac335eced54ecab3e7ca01024d5f7f3619b969 Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Thu, 21 May 2015 09:27:24 +0200 Subject: [PATCH] documentation indentation and an assert layout a small error documentation small bugs + rebase small change in the comment notation a small comment renaming the function a -> block_matrix --- include/deal.II/lac/linear_operator.h | 282 +++++++++++------- ...r_operator_15.cc => linear_operator_11.cc} | 18 +- tests/lac/linear_operator_12.cc | 18 +- tests/lac/linear_operator_13.cc | 3 +- tests/lac/linear_operator_14.cc | 27 +- .../linear_operator_14.with_cxx11=on.output | 40 +++ 6 files changed, 226 insertions(+), 162 deletions(-) rename tests/lac/{linear_operator_15.cc => linear_operator_11.cc} (88%) diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index 6cadd62cd2..73b0598c63 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -2127,234 +2127,286 @@ operator*(const PackagedOperation &comp, } - +/** + * @relates LinearOperator + * + * A function that takes a block matrix @p a and returns its associated + * lower triangular matrix operator (diagonal is not included). + * + * @code + * a00 | a01 | a02 | | + * --------------- --------------- + * a10 | a11 | a12 -> a10 | | + * --------------- --------------- + * a20 | a21 | a22 a20 | a21 | + * @endcode + * + * @ingroup LAOperators + */ template , typename Domain = Range, typename BlockMatrix> LinearOperator -lower_triangular_operator(BlockMatrix &a) +lower_triangular_operator(const BlockMatrix &block_matrix) { - Assert( a.n_block_rows() == a.n_block_cols(), - ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) ); + Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(), + ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) ); LinearOperator return_op; - return_op.reinit_range_vector = [&a](Range &v, bool fast) + return_op.reinit_range_vector = [&block_matrix](Range &v, bool fast) { // Reinitialize the block vector to have the number of blocks - // equal to the number of row blocks of the matrix a. - v.reinit(a.n_block_rows()); + // equal to the number of row blocks of the matrix block_matrix. + v.reinit(block_matrix.n_block_rows()); // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < a.n_block_rows(); ++i) - linear_operator(a.block(i,0)).reinit_range_vector(v.block(i), fast); + for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i) + linear_operator(block_matrix.block(i,0)).reinit_range_vector(v.block(i), fast); v.collect_sizes(); }; - return_op.reinit_domain_vector = [&a](Domain &v, bool fast) + return_op.reinit_domain_vector = [&block_matrix](Domain &v, bool fast) { // Reinitialize the block vector to have the number of blocks - // equal to the number of coloumn blocks of the matrix a. - v.reinit(a.n_block_cols()); + // equal to the number of coloumn blocks of the matrix block_matrix. + v.reinit(block_matrix.n_block_cols()); // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - linear_operator(a.block(0,i)).reinit_domain_vector(v.block(i), fast); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + linear_operator(block_matrix.block(0,i)).reinit_domain_vector(v.block(i), fast); v.collect_sizes(); }; - return_op.vmult = [&a](Range &v, const Domain &u) + return_op.vmult = [&block_matrix](Range &v, const Domain &u) { - Assert( v.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); - Assert( u.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + Assert( v.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); + Assert( u.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); v.block(0) *= 0; - for (unsigned int i = 1; i < a.n_block_rows(); ++i) + for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i) { - a.block(i,0).vmult(v.block(i), u.block(0)); + block_matrix.block(i,0).vmult(v.block(i), u.block(0)); for (unsigned int j = 1; j < i; ++j) - a.block(i,j).vmult_add(v.block(i), u.block(j)); + block_matrix.block(i,j).vmult_add(v.block(i), u.block(j)); } }; - return_op.vmult_add = [&a](Range &v, const Domain &u) + return_op.vmult_add = [&block_matrix](Range &v, const Domain &u) { - Assert( v.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); - Assert( u.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + Assert( v.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); + Assert( u.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - for (unsigned int i = 1; i < a.n_block_rows(); ++i) + for (unsigned int i = 1; i < block_matrix.n_block_rows(); ++i) { - a.block(i,0).vmult_add(v.block(i), u.block(0)); + block_matrix.block(i,0).vmult_add(v.block(i), u.block(0)); for (unsigned int j = 1; j < i; ++j) - a.block(i,j).vmult_add(v.block(i), u.block(j)); + block_matrix.block(i,j).vmult_add(v.block(i), u.block(j)); } }; - - return_op.Tvmult = [&a](Domain &v, const Range &u) + return_op.Tvmult = [&block_matrix](Domain &v, const Range &u) { - Assert( v.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); - Assert( u.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); - + Assert( v.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); + Assert( u.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - for (unsigned int i = 0; i < a.n_block_cols(); ++i) + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) { v.block(i) *= 0; - for (unsigned int j = i + 1; j < a.n_block_rows(); ++j) - a.block(j,i).Tvmult_add(v.block(i), u.block(j)); + for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j) + block_matrix.block(j,i).Tvmult_add(v.block(i), u.block(j)); } }; - return_op.Tvmult_add = [&a](Domain &v, const Range &u) + return_op.Tvmult_add = [&block_matrix](Domain &v, const Range &u) { - Assert( v.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); - Assert( u.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); + Assert( v.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); + Assert( u.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - for (unsigned int j = i + 1; j < a.n_block_rows(); ++j) - a.block(j,i).Tvmult_add(v.block(i), u.block(j)); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j) + block_matrix.block(j,i).Tvmult_add(v.block(i), u.block(j)); }; return return_op; } +/** + * @relates LinearOperator + * + * A function that takes a block matrix @p a and returns its associated + * upper triangular matrix operator (diagonal is not included). + * + * @code + * a00 | a01 | a02 | a01 | a02 + * --------------- --------------- + * a10 | a11 | a12 -> | | a12 + * --------------- --------------- + * a20 | a21 | a22 | | + * @endcode + * + * @ingroup LAOperators + */ template , typename Domain = Range, typename BlockMatrix> LinearOperator -upper_triangular_operator(BlockMatrix &a) +upper_triangular_operator(const BlockMatrix &block_matrix) { - Assert( a.n_block_rows() == a.n_block_cols(), - ExcDimensionMismatch(a.n_block_rows(),a.n_block_cols()) ); + Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(), + ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) ); LinearOperator return_op; - return_op.reinit_range_vector = [&a](Range &v, bool fast) + return_op.reinit_range_vector = [&block_matrix](Range &v, bool fast) { // Reinitialize the block vector to have the number of blocks - // equal to the number of row blocks of the matrix a. - v.reinit(a.n_block_rows()); + // equal to the number of row blocks of the matrix block_matrix. + v.reinit(block_matrix.n_block_rows()); // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < a.n_block_rows(); ++i) - linear_operator(a.block(i,0)).reinit_range_vector(v.block(i), fast); + for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i) + linear_operator(block_matrix.block(i,0)).reinit_range_vector(v.block(i), fast); v.collect_sizes(); }; - return_op.reinit_domain_vector = [&a](Domain &v, bool fast) + return_op.reinit_domain_vector = [&block_matrix](Domain &v, bool fast) { // Reinitialize the block vector to have the number of blocks - // equal to the number of coloumn blocks of the matrix a. - v.reinit(a.n_block_cols()); + // equal to the number of coloumn blocks of the matrix block_matrix. + v.reinit(block_matrix.n_block_cols()); // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - linear_operator(a.block(0,i)).reinit_domain_vector(v.block(i), fast); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + linear_operator(block_matrix.block(0,i)).reinit_domain_vector(v.block(i), fast); v.collect_sizes(); }; - return_op.vmult = [&a](Range &v, const Domain &u) + return_op.vmult = [&block_matrix](Range &v, const Domain &u) { - Assert( v.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); - Assert( u.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + Assert( v.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); + Assert( u.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - for (unsigned int i = 0; i < a.n_block_rows(); ++i) + for (unsigned int i = 0; i < block_matrix.n_block_rows(); ++i) { v.block(i) *= 0; - for (unsigned int j = i + 1; j < a.n_block_cols(); ++j) - a.block(i,j).Tvmult_add(v.block(i), u.block(j)); + for (unsigned int j = i + 1; j < block_matrix.n_block_cols(); ++j) + block_matrix.block(i,j).Tvmult_add(v.block(i), u.block(j)); } }; - return_op.vmult_add = [&a](Range &v, const Domain &u) + return_op.vmult_add = [&block_matrix](Range &v, const Domain &u) { - Assert( v.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_rows())); - Assert( u.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_cols())); + Assert( v.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_rows())); + Assert( u.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_cols())); - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - for (unsigned int j = i + 1; j < a.n_block_rows(); ++j) - a.block(i,j).Tvmult_add(v.block(i), u.block(j)); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + for (unsigned int j = i + 1; j < block_matrix.n_block_rows(); ++j) + block_matrix.block(i,j).Tvmult_add(v.block(i), u.block(j)); }; + return_op.Tvmult = [&block_matrix](Domain &v, const Range &u) + { + Assert( v.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); + Assert( u.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + { + v.block(i) *= 0; + for (unsigned int j = 0; j < i; ++j) + block_matrix.block(j,i).vmult_add(v.block(i), u.block(j)); + } + }; - return_op.Tvmult = [&a](Domain &v, const Range &u) - { - Assert( v.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); - Assert( u.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); - - - - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - { - v.block(i) *= 0; - for (unsigned int j = 0; j < i; ++j) - a.block(j,i).vmult_add(v.block(i), u.block(j)); - } - }; - - - return_op.Tvmult_add = [&a](Domain &v, const Range &u) + return_op.Tvmult_add = [&block_matrix](Domain &v, const Range &u) { - Assert( v.n_blocks() == a.n_block_cols(), - ExcDimensionMismatch(v.n_blocks(), a.n_block_cols())); - Assert( u.n_blocks() == a.n_block_rows(), - ExcDimensionMismatch(u.n_blocks(), a.n_block_rows())); + Assert( v.n_blocks() == block_matrix.n_block_cols(), + ExcDimensionMismatch(v.n_blocks(), block_matrix.n_block_cols())); + Assert( u.n_blocks() == block_matrix.n_block_rows(), + ExcDimensionMismatch(u.n_blocks(), block_matrix.n_block_rows())); - for (unsigned int i = 0; i < a.n_block_cols(); ++i) - for (unsigned int j = 0; j < i; ++j) - a.block(j,i).vmult_add(v.block(i), u.block(j)); + for (unsigned int i = 0; i < block_matrix.n_block_cols(); ++i) + for (unsigned int j = 0; j < i; ++j) + block_matrix.block(j,i).vmult_add(v.block(i), u.block(j)); }; return return_op; } +/** + * @relates LinearOperator + * + * Let M be a block matrix of the form Id + T made of nxn blocks + * and where T is a lower / upper triangular (without diagonal). + * Then, its inverse is of the form: + * @code + * Id + sum_{i=1}^{n-1} (-1)^i T^i + * @endcode + * This formula can be used to invert all triangular matrices + * (diagonal included of course). + * + * This function takes a block matrix @p block_matrix (possibly full) and + * a linear block operator @p inverse_diagonal made of the inverse + * of the diagonal blocks inverses. + * The output is the inverse of the matrix in the case of a triangular matrix and + * as inverse_diagonal its diagonal blocks inverses. + * Otherwise, the result is a preconditioner. + * + * The parameter @p lower is a bool that allows to specify if we want + * to use lower triangular part of @p block_matrix (true, this is the default value) or + * to use upper triangular part of @p block_matrix (false). + * @ingroup LAOperators + */ template , typename Domain = Range, typename BlockMatrix> LinearOperator -inverse_operator( BlockMatrix &a, - const LinearOperator &inverse_diagonal, - bool lower = true) +block_triangular_inverse( const BlockMatrix &block_matrix, + const LinearOperator &inverse_diagonal, + bool lower = true) { + Assert( block_matrix.n_block_rows() == block_matrix.n_block_cols(), + ExcDimensionMismatch(block_matrix.n_block_rows(),block_matrix.n_block_cols()) ); LinearOperator op_a; - if(lower) - { - op_a = lower_triangular_operator(a); - } + if (lower) + { + op_a = lower_triangular_operator(block_matrix); + } else - { - op_a = upper_triangular_operator(a); - } + { + op_a = upper_triangular_operator(block_matrix); + } auto id = identity_operator(op_a.reinit_range_vector); auto result = identity_operator(op_a.reinit_range_vector); - - for(unsigned int i = 0; i < a.n_block_cols() - 1; ++i) + // Notice that the following formula is recursive. + // we are evaluating: + // Id - T + T^2 - T^3 ... (- T)^block_matrix.n_block_cols() + for (unsigned int i = 0; i < block_matrix.n_block_cols() - 1; ++i) result = id - inverse_diagonal * op_a * result; return result * inverse_diagonal; diff --git a/tests/lac/linear_operator_15.cc b/tests/lac/linear_operator_11.cc similarity index 88% rename from tests/lac/linear_operator_15.cc rename to tests/lac/linear_operator_11.cc index 18491944f4..02c674bd7b 100644 --- a/tests/lac/linear_operator_15.cc +++ b/tests/lac/linear_operator_11.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- -// Test non symmetric variants: +// Test lower_triangular_operator: #include "../tests.h" @@ -24,7 +24,8 @@ #include #include -#define PRINTME(name, var) \ +#define PRINTME(title, name, var) \ + deallog << title << std::endl; \ deallog << "Block vector: " name << ":" << std::endl; \ for (unsigned int i = 0; i < var.n_blocks(); ++i) \ deallog << "[block " << i << " ] " << var.block(i); @@ -58,7 +59,6 @@ int main() a.block(i,j).set(0, 0, 2 + i - j); auto op_a = linear_operator>(a); - auto triangular_block_op = lower_triangular_operator< BlockVector, BlockVector, BlockSparseMatrix >(a); BlockVector u; @@ -75,31 +75,27 @@ int main() } triangular_block_op.vmult(v, u); - deallog << " -- vmult --" << std::endl; - PRINTME("v", v); + PRINTME(" -- vmult --", "v", v); // vmult_add: for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = 1; triangular_block_op.vmult_add(v, u); - deallog << " -- vmult_add --" << std::endl; - PRINTME("v", v); + PRINTME(" -- vmult_add --", "v", v); // Tvmult for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = i+1; triangular_block_op.Tvmult(u, v); - deallog << " -- Tvmult --" << std::endl; - PRINTME("u", u); + PRINTME(" -- Tvmult --", "u", u); // Tvmult_add for(unsigned int i = 0; i<3; ++i) u.block(i)[0] = 1; triangular_block_op.Tvmult_add(u, v); - deallog << " -- Tvmult_add --" << std::endl; - PRINTME("u", u); + PRINTME(" -- Tvmult_add --", "u", u); } } diff --git a/tests/lac/linear_operator_12.cc b/tests/lac/linear_operator_12.cc index f76d503b95..1b6601d624 100644 --- a/tests/lac/linear_operator_12.cc +++ b/tests/lac/linear_operator_12.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- -// Test non symmetric variants: +// Test upper_triangular_operator: #include "../tests.h" @@ -24,12 +24,12 @@ #include #include -#define PRINTME(name, var) \ +#define PRINTME(title, name, var) \ + deallog << title << std::endl; \ deallog << "Block vector: " name << ":" << std::endl; \ for (unsigned int i = 0; i < var.n_blocks(); ++i) \ deallog << "[block " << i << " ] " << var.block(i); - using namespace dealii; int main() @@ -75,31 +75,27 @@ int main() } triangular_block_op.vmult(v, u); - deallog << " -- vmult --" << std::endl; - PRINTME("v", v); + PRINTME(" -- vmult --","v", v); // vmult_add: for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = 1; triangular_block_op.vmult_add(v, u); - deallog << " -- vmult_add --" << std::endl; - PRINTME("v", v); + PRINTME(" -- vmult_add --", "v", v); // Tvmult for(unsigned int i = 0; i<3; ++i) v.block(i)[0] = i+1; triangular_block_op.Tvmult(u, v); - deallog << " -- Tvmult --" << std::endl; - PRINTME("u", u); + PRINTME(" -- Tvmult --", "u", u); // Tvmult_add for(unsigned int i = 0; i<3; ++i) u.block(i)[0] = 1; triangular_block_op.Tvmult_add(u, v); - deallog << " -- Tvmult_add --" << std::endl; - PRINTME("u", u); + PRINTME(" -- Tvmult_add --", "u", u); } } diff --git a/tests/lac/linear_operator_13.cc b/tests/lac/linear_operator_13.cc index 199bd702d3..61488e2f89 100644 --- a/tests/lac/linear_operator_13.cc +++ b/tests/lac/linear_operator_13.cc @@ -13,7 +13,8 @@ // // --------------------------------------------------------------------- -// Test non symmetric variants: +// Test upper_triangular_block_op and lower_triangular_operator +// simultaneously: #include "../tests.h" diff --git a/tests/lac/linear_operator_14.cc b/tests/lac/linear_operator_14.cc index a722ad1d7a..7a49f4574b 100644 --- a/tests/lac/linear_operator_14.cc +++ b/tests/lac/linear_operator_14.cc @@ -13,7 +13,7 @@ // // --------------------------------------------------------------------- -// Test non symmetric variants: +// Test preconditioner_from_diagonal_inverse: #include "../tests.h" @@ -60,13 +60,11 @@ int main() BlockSparseMatrix d(sparsity_pattern); for (unsigned int i = 0; i < 3; ++i) d.set(i,i, 1.0 / (i+i +1) ); - // std::array diag {{a.block(0,0), a.block(1,1),a.block(2,2)}}; - auto op_a = linear_operator< BlockVector >(a); auto diagonal_inv = linear_operator< BlockVector >(d); // auto diagonal_inv = block_diagonal_operator(diag); - auto inverse_op_a = inverse_operator< BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); + auto inverse_op_a = block_triangular_inverse< BlockVector, BlockVector, BlockSparseMatrix >(a, diagonal_inv); auto identity = inverse_op_a * op_a; @@ -78,7 +76,6 @@ int main() op_a.reinit_range_vector(v, false); for(unsigned int j = 0; j<3; ++j) { - // check 1 for(unsigned int i = 0; i<3; ++i) { u.block(i)[0] = 0;; @@ -86,9 +83,7 @@ int main() } u.block(j)[0] = 1;; - op_a.vmult(v, u); - // identity.Tvmult(v2, u2); PRINTME("v", v); } @@ -98,17 +93,14 @@ int main() inverse_op_a.reinit_range_vector(v, false); for(unsigned int j = 0; j<3; ++j) { - // check 1 for(unsigned int i = 0; i<3; ++i) { - u.block(i)[0] = 0;; + u.block(i)[0] = 0; v.block(i)[0] = 0; } u.block(j)[0] = 1;; - inverse_op_a.vmult(v, u); - // identity.Tvmult(v2, u2); PRINTME("v", v); } @@ -118,7 +110,6 @@ int main() identity.reinit_range_vector(v, false); for(unsigned int j = 0; j<3; ++j) { - // check 1 for(unsigned int i = 0; i<3; ++i) { u.block(i)[0] = 0;; @@ -126,21 +117,9 @@ int main() } u.block(j)[0] = 1;; - identity.vmult(v, u); - // identity.Tvmult(v2, u2); PRINTME("v", v); } - - // op_a.vmult(v, u); - // // identity.Tvmult(v2, u2); - // - // PRINTME("v", v); - // - // inverse_op_a.vmult(v, u); - // // identity.Tvmult(v2, u2); - // - // PRINTME("v", v); } } diff --git a/tests/lac/linear_operator_14.with_cxx11=on.output b/tests/lac/linear_operator_14.with_cxx11=on.output index e69de29bb2..f4633ca2e2 100644 --- a/tests/lac/linear_operator_14.with_cxx11=on.output +++ b/tests/lac/linear_operator_14.with_cxx11=on.output @@ -0,0 +1,40 @@ + +DEAL:: -- Matrix -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 10.0000000000 +DEAL::[block 2 ] 10.0000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 3.00000000000 +DEAL::[block 2 ] 10.0000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 5.00000000000 +DEAL:: -- Inverse -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] -3.33333333333 +DEAL::[block 2 ] 4.66666666667 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.333333333333 +DEAL::[block 2 ] -0.666666666667 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.200000000000 +DEAL:: -- Identity -- +DEAL::Block vector: v: +DEAL::[block 0 ] 1.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 1.00000000000 +DEAL::[block 2 ] 0.00000000000 +DEAL::Block vector: v: +DEAL::[block 0 ] 0.00000000000 +DEAL::[block 1 ] 0.00000000000 +DEAL::[block 2 ] 1.00000000000 -- 2.39.5