From b1a8baa9177b6e49bdd407474733bd8775b56b51 Mon Sep 17 00:00:00 2001 From: kanschat Date: Tue, 26 Sep 2006 16:25:00 +0000 Subject: [PATCH] add eigenvector computation to compute_eigenvalues access to eigenvectors is missing yet and needs thinking git-svn-id: https://svn.dealii.org/trunk@13962 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/lapack_full_matrix.h | 21 ++++++++++++++++++-- deal.II/lac/source/lapack_full_matrix.cc | 20 +++++++++++++------ 2 files changed, 33 insertions(+), 8 deletions(-) diff --git a/deal.II/lac/include/lac/lapack_full_matrix.h b/deal.II/lac/include/lac/lapack_full_matrix.h index ef953b3cc2..4eb12002b2 100644 --- a/deal.II/lac/include/lac/lapack_full_matrix.h +++ b/deal.II/lac/include/lac/lapack_full_matrix.h @@ -218,10 +218,15 @@ class LAPACKFullMatrix : public TransposeTable * LAPACKSupport::unusable after * this operation. * + * The optional arguments allow + * to compute left and right + * eigenvectors as well. + * * @note Calls the LAPACK * function Xgeev. */ - void compute_eigenvalues (); + void compute_eigenvalues (const bool right_eigenvectors = false, + const bool left_eigenvectors = false); /** * Retrieve eigenvalue after @@ -234,7 +239,7 @@ class LAPACKFullMatrix : public TransposeTable /** * Print the matrix and allow * formatting of entries. - * + * * The parameters allow for a * flexible setting of the output * format: @@ -302,6 +307,7 @@ class LAPACKFullMatrix : public TransposeTable * some LAPACK functions. */ mutable std::vector work; + /** * Real parts of * eigenvalues. Filled by @@ -316,6 +322,17 @@ class LAPACKFullMatrix : public TransposeTable */ std::vector wi; + /** + * Space where left eigenvectors + * can be stored. + */ + std::vector vl; + + /** + * Space where right eigenvectors + * can be stored. + */ + std::vector vr; }; /**@}*/ diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index ab1ac8f74f..723d0b96a5 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2005 by the deal.II authors +// Copyright (C) 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -154,15 +154,23 @@ LAPACKFullMatrix::Tvmult ( template void -LAPACKFullMatrix::compute_eigenvalues() +LAPACKFullMatrix::compute_eigenvalues( + const bool right, + const bool left) { const int nn = this->n_cols(); wr.resize(nn); wi.resize(nn); + if (right) vr.resize(nn*nn); + if (left) vl.resize(nn*nn); + number* values = const_cast (this->data()); int info = 0; int lwork = 1; + const char * const jobvr = (right) ? (&V) : (&N); + const char * const jobvl = (left) ? (&V) : (&N); + // Optimal workspace query: // The LAPACK routine DGEEV requires @@ -186,9 +194,9 @@ LAPACKFullMatrix::compute_eigenvalues() lwork = -1; work.resize(1); - geev(&N, &N, &nn, values, &nn, + geev(jobvl, jobvr, &nn, values, &nn, &wr[0], &wi[0], - 0, &one, 0, &one, + &vl[0], &nn, &vr[0], &nn, &work[0], &lwork, &info); // geev returns info=0 on // success. Since we only queried @@ -206,9 +214,9 @@ LAPACKFullMatrix::compute_eigenvalues() work.resize((unsigned int) lwork); // Finally compute the eigenvalues. - geev(&N, &N, &nn, values, &nn, + geev(jobvl, jobvr, &nn, values, &nn, &wr[0], &wi[0], - 0, &one, 0, &one, + &vl[0], &nn, &vr[0], &nn, &work[0], &lwork, &info); // Negative return value implies a // wrong argument. This should be -- 2.39.5