#include <lac/lapack_full_matrix.h>
#include <lac/lapack_templates.h>
#include <lac/lapack_support.h>
+#include <lac/full_matrix.h>
#include <lac/vector.h>
using namespace LAPACKSupport;
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator = (const FullMatrix<number2>& M)
{
- Assert (n() == M.m(), ExcDimensionMismatch(n(), M.m()));
- Assert (n() == M.n(), ExcDimensionMismatch(n(), M.n()));
- for (unsigned int i=0;i<m;++i)
- for (unsigned int j=0;j<n;++j)
+ Assert (n_rows() == M.m(), ExcDimensionMismatch(n_rows(), M.m()));
+ Assert (n_cols() == M.n(), ExcDimensionMismatch(n_cols(), M.n()));
+ for (unsigned int i=0;i<n_rows();++i)
+ for (unsigned int j=0;j<n_cols();++j)
(*this)(i,j) = M(i,j);
return *this;
}
template LAPACKFullMatrix<double>;
+template LAPACKFullMatrix<double> &
+LAPACKFullMatrix<double>::operator = (const FullMatrix<double>& M);
+template LAPACKFullMatrix<double> &
+LAPACKFullMatrix<double>::operator = (const FullMatrix<float>& M);
--- /dev/null
+//--------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//--------------------------------------------------------------------
+
+// Tests all aspects of LAPACKFullMatrix
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <lac/lapack_full_matrix.h>
+#include <lac/full_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+const double symm[] =
+{
+ 4., 1., -1., 1.,
+ 1., 4., 1., -1.,
+ -1., 1., 4., 1.,
+ 1., -1., 4., 1.
+};
+
+const double rect[] =
+{
+ 4., 3., 2., 1.,
+ 5., 8., 1., -2.,
+ 11., 13., -4., -5
+};
+
+
+int main()
+{
+ std::ofstream logfile("lapack.output");
+ logfile.precision(4);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+#ifdef HAVE_LIBLAPACK
+ FullMatrix<double> A(3,4,rect);
+ LAPACKFullMatrix<double> LA(3,4);
+ LA = A;
+
+ Vector<double> u(4);
+ Vector<double> v1(3);
+ Vector<double> v2(3);
+
+ for (unsigned int i=0;i<u.size();++i)
+ u(i) = i*i;
+
+ // Test rectangular vmult. All
+ // results compare with same
+ // operation for FullMatrix.
+ deallog.push("Rect");
+
+ A.vmult(v1,u);
+ LA.vmult(v2,u);
+ v1 -= v2;
+ if (v1.l2_norm() < 1.e-14)
+ deallog << "vmult ok" << std::endl;
+ v1 = v2;
+
+ A.vmult_add(v1,u);
+ LA.vmult_add(v2,u);
+ v1 -= v2;
+ if (v1.l2_norm() < 1.e-14)
+ deallog << "vmult_add ok" << std::endl;
+
+ LA.Tvmult(u, v2);
+ u *= -1;
+ A.Tvmult_add(u, v2);
+ if (u.l2_norm() < 1.e-14)
+ deallog << "Tvmult ok" << std::endl;
+
+ A.Tvmult(u, v2);
+ u *= -1;
+ LA.Tvmult_add(u, v2);
+ if (u.l2_norm() < 1.e-14)
+ deallog << "Tvmult_add ok" << std::endl;
+
+ deallog.pop();
+
+ // Test symmetric system
+ A.reinit(4,4);
+ LA.reinit(4,4);
+ A.fill(symm);
+ LA = A;
+ v1.reinit(4);
+ v2.reinit(4);
+
+#else
+ deallog.push("Rect");
+ deallog << "vmult ok" << std::endl;
+ deallog << "vmult_add ok" << std::endl;
+ deallog << "Tvmult ok" << std::endl;
+ deallog << "Tvmult_add ok" << std::endl;
+ deallog.pop();
+#endif
+}