From 9abc705db5177a878f5206d8f10f0c84cef8fb78 Mon Sep 17 00:00:00 2001 From: janssen Date: Mon, 20 Sep 2010 16:40:45 +0000 Subject: [PATCH] include a class to be able to interface to ARPACK and some documentation how to install ARPACK and run it with deal.II git-svn-id: https://svn.dealii.org/trunk@22085 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/external-libs/arpack.html | 76 ++++ deal.II/doc/readme.html | 14 + deal.II/lac/include/lac/arpack_solver.h | 451 ++++++++++++++++++++++++ 3 files changed, 541 insertions(+) create mode 100644 deal.II/doc/external-libs/arpack.html create mode 100644 deal.II/lac/include/lac/arpack_solver.h diff --git a/deal.II/doc/external-libs/arpack.html b/deal.II/doc/external-libs/arpack.html new file mode 100644 index 0000000000..b42ca13c63 --- /dev/null +++ b/deal.II/doc/external-libs/arpack.html @@ -0,0 +1,76 @@ + + + + + The deal.II Readme on interfacing to ARPACK + + + + + + + + + + +

Using and Installing ARPACK

+

ARPACK + is a collection of Fortran77 subroutines designed to solve large + scale eigenvalue problems. + + After you downloaded the Fortran version of ARPACK + and the patch unzip the files you got. That will create a directory + named ARPACK. If you need further instructions + please read the README file or the instructions. We will explain + here in a few steps what has to be done to be able to compile + ARPACK. +

+ + + +

+ Note: For compilation of ARPACK we emphasise + adding the compiler flag -fPIC. This is a definate + requirement if we are compiling deal.II with + shared libraries (which is the default). If we had preferred to be + compiling deal.II without shared libraries, + that's ok too; in that case we would do exactly the same thing + as described above, but this time omitting + the -fPIC flag from the scheme. +

+ +

+ Try to run one of the examples and compare the output. + How the output should look like is stated in the README + that can be found in the EXAMPLES directory. +

+ +

+ If that output you produced looks like it should you can + proceed to comile deal.II with + ARPACK. +

+ +

Interfacing deal.II + to ARPACK

+ +

To be able to use the interface of deal.II + for ARPACK we then + configure deal.II with the following + additional option: + --with-arpack=path/to/arpack +

+ + diff --git a/deal.II/doc/readme.html b/deal.II/doc/readme.html index fd626a5a8c..5718222251 100644 --- a/deal.II/doc/readme.html +++ b/deal.II/doc/readme.html @@ -751,6 +751,20 @@ Juha Nieminen. + +
ARPACK
+
+ There is a wrapper for the + ARPACK + library, which has its own license; if you want to use it with deal.II, + please read it and make sure that you agree with it. + You can find the license of ARPACK + here. + For a detailed description of how to compile ARPACK and linking with deal.II, see + this page. +
+ + diff --git a/deal.II/lac/include/lac/arpack_solver.h b/deal.II/lac/include/lac/arpack_solver.h new file mode 100644 index 0000000000..38b7923a48 --- /dev/null +++ b/deal.II/lac/include/lac/arpack_solver.h @@ -0,0 +1,451 @@ +//--------------------------------------------------------------------------- +// $Id: $ +// Version: $Name$ +// Authors: Bärbel Janssen, University of Heidelberg, +// Agnieszka Miedler, TU Berlin, 2010 +// +// Copyright (C) 2010 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__arpack_solver_h +#define __deal2__arpack_solver_h + +#include + +#ifdef DEAL_II_USE_ARPACK + +DEAL_II_NAMESPACE_OPEN + +extern "C" void dnaupd_(int *ido, char *bmat, const unsigned int *n, char *which, + const unsigned int *nev, const double *tol, double *resid, int *ncv, + double *v, int *ldv, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, + int *info); + +extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d, + double *di, double *z, int *ldz, double *sigmar, + double *sigmai, double* workev, char *bmat,const unsigned int *n, char *which, + const unsigned int *nev, const double *tol, double *resid, int *ncv, + double *v, int *ldv, int *iparam, int *ipntr, + double *workd, double *workl, int *lworkl, int *info); + + +class ArpackSolver : public Subscriptor +{ + public: + enum which + { + algebraically_largest, + algebraically_smallest, + largest_magnitude, + smallest_magnitude, + largest_real_part, + smallest_real_part, + largest_imaginary_part, + smallest_imaginary_part, + both_ends + }; + + struct AdditionalData + { + const unsigned int number_of_arnoldi_vectors; + const which eigenvalue_of_interest; + AdditionalData( + const unsigned int number_of_arnoldi_vectors = 15, + const which eigenvalue_of_interest = largest_magnitude); + }; + SolverControl &control () const; + + ArpackSolver(SolverControl& control, + const AdditionalData& data = AdditionalData()); + + template + void solve( + const MATRIX &A, + const MATRIX &B, + const INVERSE& inverse, + std::vector &eigenvalues_r, + std::vector &eigenvalues_i, + std::vector &eigenvectors, + const unsigned int n_eigenvalues); + + template + void solve( + const MATRIX &A, + const INVERSE& inverse, + std::vector &eigenvalues_r, + std::vector &eigenvalues_i, + std::vector &eigenvectors, + const unsigned int n_eigenvalues); + protected: + + /** + * Reference to the object that + * controls convergence of the + * iterative solver. + */ + SolverControl &solver_control; + + /** + * Store a copy of the flags for + * this particular solver. + */ + const AdditionalData additional_data; + private: + DeclException2 (ExcInvalidNumberofEigenvalues, int, int, + << "Number of wanted eigenvalues " << arg1 + << " is larger that the size of the matrix " << arg2); + + DeclException2 (ExcInvalidNumberofArnoldiVectors, int, int, + << "Number of Arnoldi vectors " << arg1 + << " is larger that the size of the matrix " << arg2); + + DeclException2 (ExcSmallNumberofArnoldiVectors, int, int, + << "Number of Arnoldi vectors " << arg1 + << " is too small to obtain " << arg2 + << " eigenvalues"); + + DeclException1 (ExcArpackIdo, int, << "This ido " << arg1 + << " is not supported. Check documentation of ARPACK"); + + DeclException1 (ExcArpackMode, int, << "This mode " << arg1 + << " is not supported. Check documentation of ARPACK"); + + DeclException1 (ExcArpackInfodsaupd, int, + << "Error with dsaupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (ExcArpackInfodneupd, int, + << "Error with dneupd, info " << arg1 + << ". Check documentation of ARPACK"); + + DeclException1 (ExcArpackInfoMaxIt, int, + << "Maximum number " << arg1 + << " of iterations reached."); + + DeclException1 (ExcArpackNoShifts, int, + << "No shifts could be applied during implicit" + << " Arnoldi update, try increasing the number of" + << " Arnoldi vectors."); +}; + +ArpackSolver::AdditionalData:: +AdditionalData (const unsigned int number_of_arnoldi_vectors, + const which eigenvalue_of_interest) + : + number_of_arnoldi_vectors(number_of_arnoldi_vectors), + eigenvalue_of_interest(eigenvalue_of_interest) +{} + +ArpackSolver::ArpackSolver (SolverControl& control, + const AdditionalData& data) + : + solver_control (control), + additional_data (data) + +{} + +template +void ArpackSolver::solve ( + const MATRIX &system_matrix, + const MATRIX &mass_matrix, + const INVERSE& inverse, + std::vector &eigenvalues_real, + std::vector &eigenvalues_im, + std::vector &eigenvectors, + const unsigned int n_eigenvalues) +{ + //inside the routines of ARPACK the + //values change magically, so store + //them here + + const unsigned int n = system_matrix.m(); + const unsigned int n_inside_arpack = system_matrix.m(); + + /* + if(n < 0 || nev <0 || p < 0 || maxit < 0 ) + std:cout << "All input parameters have to be positive.\n"; + */ + Assert (n_eigenvalues < n, + ExcInvalidNumberofEigenvalues(n_eigenvalues, n)); + + Assert (additional_data.number_of_arnoldi_vectors < n, + ExcInvalidNumberofArnoldiVectors( + additional_data.number_of_arnoldi_vectors, n)); + + Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1, + ExcSmallNumberofArnoldiVectors( + additional_data.number_of_arnoldi_vectors, n_eigenvalues)); + // ARPACK mode for dnaupd, here only mode 3 + int mode = 3; + + // reverse communication parameter + int ido = 0; + + /** + * 'G' generalized eigenvalue problem + * 'I' standard eigenvalue problem + */ + char bmat[2] = "G"; + + /** Specify the eigenvalues of interest, + * possible parameters + * "LA" algebraically largest + * "SA" algebraically smallest + * "LM" largest magnitude + * "SM" smallest magnitude + * "LR" largest real part + * "SR" smallest real part + * "LI" largest imaginary part + * "SI" smallest imaginary part + * "BE" both ends of spectrum simultaneous + */ + char which[3]; + switch(additional_data.eigenvalue_of_interest) + { + case algebraically_largest: + strcpy (which, "LA"); + break; + case algebraically_smallest: + strcpy (which, "SA"); + break; + case largest_magnitude: + strcpy (which, "LM"); + break; + case smallest_magnitude: + strcpy (which, "SM"); + break; + case largest_real_part: + strcpy (which, "LR"); + break; + case smallest_real_part: + strcpy (which, "SR"); + break; + case largest_imaginary_part: + strcpy (which, "LI"); + break; + case smallest_imaginary_part: + strcpy (which, "SI"); + break; + case both_ends: + strcpy (which, "BE"); + break; + } + + // tolerance for ARPACK + const double tol = control().tolerance(); + + // if the starting vector is used it has to be in resid + std::vector resid(n, 1.); + + // number of Arnoldi basis vectors specified by user in p + int ncv = additional_data.number_of_arnoldi_vectors; + + int ldv = n; + std::vector v (ldv*ncv, 0.0); + + //information to the routines + std::vector iparam (11, 0); + + iparam[0] = 1; // shift strategy + + // maximum number of iterations + iparam[2] = control().max_steps(); + + /** Sets the mode of dsaupd. + * 1 is exact shifting, + * 2 is user-supplied shifts, + * 3 is shift-invert mode, + * 4 is buckling mode, + * 5 is Cayley mode. + */ + + iparam[6] = mode; + std::vector ipntr (14, 0); + + // work arrays for ARPACK + double *workd; + workd = new double[3*n]; + + for(unsigned int i=0; i<3*n; ++i) + workd[i] = 0.0; + + int lworkl = 3*ncv*(ncv + 2); + std::vector workl (lworkl, 0.); + //information out of the iteration + int info = 1; + + const unsigned int nev = n_eigenvalues; + while(ido != 99) + { + // call of ARPACK dnaupd routine + dnaupd_(&ido, bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0], + workd, &workl[0], &lworkl, &info); + + if(ido == 99) + break; + + switch(mode) + { + case 3: + { + switch(ido) + { + case -1: + { + + VECTOR src,dst,tmp; + src.reinit(eigenvectors[0]); + dst.reinit(src); + tmp.reinit(src); + + + for(unsigned int i=0; i select (ncv, 0); + + // Arnoldi basis vectors + int ldz = n; + + std::vector z (ldz*ncv, 0.); + + double sigmar = 0.0; // real part of the shift + double sigmai = 0.0; // imaginary part of the shift + + int lworkev = 3*ncv; + std::vector workev (lworkev, 0.); + + // call of ARPACK dneupd routine + dneupd_(&rvec, howmany, &select[0], &eigenvalues_real[0], + &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai, + &workev[0], bmat, &n_inside_arpack, which, &nev, &tol, + &resid[0], &ncv, &v[0], &ldv, + &iparam[0], &ipntr[0], workd, &workl[0], &lworkl, &info); + + if (info == 1) + { + Assert (false, ExcArpackInfoMaxIt(control().max_steps())); + } + else if (info == 3) + { + Assert (false, ExcArpackNoShifts(1));; + } + else if (info!=0) + { + Assert (false, ExcArpackInfodneupd(info)); + } + + for (unsigned int i=0; i +void ArpackSolver::solve ( + const MATRIX &system_matrix, + const INVERSE& inverse, + std::vector &eigenvalues_real, + std::vector &eigenvalues_im, + std::vector &eigenvectors, + const unsigned int n_eigenvalues) +{ + solve (system_matrix, IdentityMatrix(system_matrix.n()), + inverse, eigenvalues_real, eigenvalues_im, n_eigenvalues); +} + +SolverControl& ArpackSolver::control () const +{ + return solver_control; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif +#endif -- 2.39.5