From: Denis Davydov Date: Wed, 6 Dec 2017 13:44:22 +0000 (+0100) Subject: lapack: add vmult() Tvmult() for triangular matrices X-Git-Tag: v9.0.0-rc1~684^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bcf3958aea24bded6d96642c70e5550fb2093bcc;p=dealii.git lapack: add vmult() Tvmult() for triangular matrices --- diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index 344f09ec57..3cfc34d741 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -206,6 +206,35 @@ LAPACKFullMatrix::vmult ( const number beta = (adding ? 1. : 0.); const number null = 0.; + // use trmv for triangular matrices + if ((property == upper_triangular || + property == lower_triangular) && + (mm==nn) && + state == matrix) + { + Assert (adding == false, + ExcNotImplemented()); + + AssertDimension(v.size(), this->n_cols()); + AssertDimension(w.size(), this->n_rows()); + + const char diag = 'N'; + const char trans = 'N'; + const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L); + + w = v; + + const int N = mm; + const int lda = N; + const int incx = 1; + + trmv (&uplo, &trans, &diag, + &N, &this->values[0], &lda, + &w[0], &incx); + + return; + } + switch (state) { case matrix: @@ -266,6 +295,36 @@ LAPACKFullMatrix::Tvmult ( const number beta = (adding ? 1. : 0.); const number null = 0.; + // use trmv for triangular matrices + if ((property == upper_triangular || + property == lower_triangular) && + (mm==nn) && + state == matrix) + { + Assert (adding == false, + ExcNotImplemented()); + + AssertDimension(v.size(), this->n_cols()); + AssertDimension(w.size(), this->n_rows()); + + const char diag = 'N'; + const char trans = 'T'; + const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L); + + w = v; + + const int N = mm; + const int lda = N; + const int incx = 1; + + trmv (&uplo, &trans, &diag, + &N, &this->values[0], &lda, + &w[0], &incx); + + return; + } + + switch (state) { case matrix: diff --git a/tests/lapack/full_matrix_23.cc b/tests/lapack/full_matrix_23.cc new file mode 100644 index 0000000000..2356fb5843 --- /dev/null +++ b/tests/lapack/full_matrix_23.cc @@ -0,0 +1,89 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::vmult() and Tvmult() for triangular matrices + +/* MWE for size=3 in Octave: +R = [1,2,3; 0, 5, 6; 0, 0, 9] +x = [2; -7; 1] + +> R*x +ans = + + -9 + -29 + 9 + +R' * x +ans = + + 2 + -31 + -27 +*/ + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test() +{ + const unsigned int size = 3; + LAPACKFullMatrix M(size); + M.set_property(LAPACKSupport::upper_triangular); + + M = 0.; + unsigned int counter = 1; + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = 0; j < size; ++j) + { + if (j >= i) + M(i,j) = counter; + + counter++; + } + + Vector x(size), y(size); + x[0] = 2; + x[1] = -7; + x[2] = 1; + + M.vmult(y,x); + y.print(deallog.get_file_stream(), 6, false); + + M.Tvmult(y,x); + y.print(deallog.get_file_stream(), 6, false); + +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + test(); + +} diff --git a/tests/lapack/full_matrix_23.output b/tests/lapack/full_matrix_23.output new file mode 100644 index 0000000000..6d295e6154 --- /dev/null +++ b/tests/lapack/full_matrix_23.output @@ -0,0 +1,3 @@ + +-9.000000 -29.000000 9.000000 +2.000000 -31.000000 -27.000000