lapack: add vmult() Tvmult() for triangular matrices

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index 344f09ec5711ffbc5bda79d70d6f967b89572497..3cfc34d7417164a5ca3e083844f847d26de837d6 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -206,6 +206,35 @@ LAPACKFullMatrix<number>::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<number>::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 (file)
index 0000000..2356fb5
--- /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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test()
+{
+  const unsigned int size = 3;
+  LAPACKFullMatrix<NumberType> 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<NumberType> 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<double>();
+
+}

diff --git a/tests/lapack/full_matrix_23.output b/tests/lapack/full_matrix_23.output
new file mode 100644 (file)
index 0000000..6d295e6
--- /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