]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
include a class to be able to interface to ARPACK and some documentation
authorjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 Sep 2010 16:40:45 +0000 (16:40 +0000)
committerjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 20 Sep 2010 16:40:45 +0000 (16:40 +0000)
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 [new file with mode: 0644]
deal.II/doc/readme.html
deal.II/lac/include/lac/arpack_solver.h [new file with mode: 0644]

diff --git a/deal.II/doc/external-libs/arpack.html b/deal.II/doc/external-libs/arpack.html
new file mode 100644 (file)
index 0000000..b42ca13
--- /dev/null
@@ -0,0 +1,76 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"                  
+         "http://www.w3.org/TR/html4/loose.dtd">                        
+
+<html>                                                                          
+  <head>                            
+    <title>The deal.II Readme on interfacing to ARPACK</title>
+    <link href="screen.css" rel="StyleSheet">
+    <meta name="author" content="the deal.II authors <authors @ dealii.org>"> 
+    <meta name="copyright" content="Copyright (C) 2010 by the deal.II authors">
+    <meta name="date" content="$Date: 2010-05-05 (Thu, 05 May 2010) $">
+    <meta name="svn_id" content="$Id: readme-arpack.html$">
+    <meta name="keywords" content="deal.II"> 
+  </head>                                                                      
+  <body> 
+
+    
+    <h1>Using and Installing <acronym>ARPACK</acronym></h1>
+    <p> <a href="http://www.caam.rice.edu/software/ARPACK/">ARPACK</a> 
+    is a collection of Fortran77 subroutines designed to solve large
+    scale eigenvalue problems.
+
+    After you downloaded the Fortran version of <acronym>ARPACK</acronym>
+    and the patch unzip the files you got. That will create a directory 
+    named <acronym>ARPACK</acronym>. If you need further instructions 
+    please read the README file or the <a href="http://www.caam.rice.edu/software/ARPACK/SRC/instruction.arpack">instructions</a>. We will explain 
+    here in a few steps what has to be done to be able to compile 
+    <acronym>ARPACK</acronym>.
+    </p>
+
+    <ul>
+      <li>edit ARmake.inc</li>
+      <ul>
+        <li>change home to the correct directory</li>
+        <li>choose which <acronym>BLAS</acronym> and 
+        <acronym>LAPACK</acronym> you would like to use</li>
+      </ul>
+      <li>change the file second.f in the UTIL directory</li>
+      <li>do <code>make lib</code> in the current directory 
+      to build the standard library <code>libarpack_$(PLAT).a</code>
+      </li>
+    </ul>
+
+    <p>
+    Note: For compilation of <acronym>ARPACK</acronym> we emphasise
+    adding the compiler flag <code>-fPIC</code>. This is a definate
+    requirement if we are compiling <acronym>deal.II</acronym> with
+    shared libraries (which is the default). If we had preferred to be
+    compiling <acronym>deal.II</acronym> 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 <code>-fPIC</code> flag from the scheme.
+    </p>
+
+    <p>
+    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 <code>EXAMPLES</code> directory.
+    </p> 
+
+    <p>
+    If that output you produced looks like it should you can 
+    proceed to comile <acronym>deal.II</acronym> with 
+    <acronym>ARPACK</acronym>.
+    </p>
+
+    <h1>Interfacing <acronym>deal.II</acronym>
+      to <acronym>ARPACK</acronym></h1>
+
+    <p> To be able to use the interface of <acronym>deal.II</acronym>
+    for <acronym>ARPACK</acronym> we then
+      configure <acronym>deal.II</acronym> with the following
+      additional option:
+       <code>--with-arpack=path/to/arpack</code>
+    </p>
+  </body>
+</html>
index fd626a5a8c35d2dceb4a8c1092720ab9510930b6..571822225109c68ad560a5e233566e56e5e0cc0b 100644 (file)
       <a href="http://www.students.tut.fi/~warp/">Juha Nieminen</a>.
     </dd>
 
+    <a name="ARPACK"></a>
+    <dt>ARPACK</dt>
+    <dd>
+      There is a wrapper for the 
+      <a href="http://www.caam.rice.edu/software/ARPACK/">ARPACK</a>
+      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
+      <a href="http://www.caam.rice.edu/software/ARPACK/RiceBSD.txt">here</a>.
+      For a detailed description of how to compile ARPACK and linking with deal.II, see
+      <a href="external-libs/arpack.html" target="body">this page</a>.
+    </dd>
+
+
 
     </dl>
 
diff --git a/deal.II/lac/include/lac/arpack_solver.h b/deal.II/lac/include/lac/arpack_solver.h
new file mode 100644 (file)
index 0000000..38b7923
--- /dev/null
@@ -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 <base/config.h>
+
+#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 <typename VECTOR, typename MATRIX, typename INVERSE>
+    void solve(
+        const MATRIX &A,
+        const MATRIX &B,
+        const INVERSE& inverse,
+        std::vector<double> &eigenvalues_r, 
+        std::vector<double> &eigenvalues_i, 
+        std::vector<VECTOR> &eigenvectors,
+        const unsigned int n_eigenvalues);
+
+    template <typename VECTOR, typename MATRIX, typename INVERSE>
+    void solve(
+        const MATRIX &A,
+        const INVERSE& inverse,
+        std::vector<double> &eigenvalues_r, 
+        std::vector<double> &eigenvalues_i, 
+        std::vector<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 <typename VECTOR, typename MATRIX, typename INVERSE>
+void ArpackSolver::solve (
+    const MATRIX &system_matrix,
+    const MATRIX &mass_matrix,
+    const INVERSE& inverse,
+    std::vector<double> &eigenvalues_real, 
+    std::vector<double> &eigenvalues_im, 
+    std::vector<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<double> 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<double> v (ldv*ncv, 0.0);
+
+  //information to the routines
+  std::vector<int> 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<int> 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<double> 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<src.size(); ++i)
+                  src(i) = *(workd+ipntr[0]-1+i);
+                 
+                // multiplication with mass matrix M
+                mass_matrix.vmult(tmp, src);
+                // solving linear system
+                inverse.vmult(dst,tmp);
+
+                for(unsigned int i=0; i<dst.size(); ++i)
+                  *(workd+ipntr[1]-1+i) = dst(i);
+              }
+              break;
+
+            case  1:
+              {
+
+                VECTOR src,dst,tmp, tmp2;
+                src.reinit(eigenvectors[0]);
+                dst.reinit(src);
+                tmp.reinit(src);
+                tmp2.reinit(src);
+          
+                for(unsigned int i=0; i<src.size(); ++i)
+                {
+                  src(i) = *(workd+ipntr[2]-1+i);
+                  tmp(i) = *(workd+ipntr[0]-1+i);
+                }                              
+                // solving linear system
+                inverse.vmult(dst,src);
+
+                for(unsigned int i=0; i<dst.size(); ++i)
+                  *(workd+ipntr[1]-1+i) = dst(i);
+              }
+              break;
+
+            case  2: 
+              {
+
+                VECTOR src,dst;
+                src.reinit(eigenvectors[0]);
+                dst.reinit(src);
+
+                for(unsigned int i=0; i<src.size(); ++i)
+                  src(i) = *(workd+ipntr[0]-1+i);
+                
+                // Multiplication with mass matrix M
+                mass_matrix.vmult(dst, src);
+
+                for(unsigned int i=0; i<dst.size(); ++i)
+                  *(workd+ipntr[1]-1+i) = dst(i);
+
+              }
+              break;
+
+            default:
+              Assert (false, ExcArpackIdo(ido));
+              break;
+          }
+        }
+        break;
+      default: 
+        Assert (false, ExcArpackMode(mode));
+        break;
+    }
+  }
+
+  if (info<0) 
+  {
+    Assert (false, ExcArpackInfodsaupd(info));
+  }
+  else
+  {
+    /** 1 - compute eigenvectors,
+     *  0 - only eigenvalues
+     */
+    int rvec = 1; 
+
+    // which eigenvectors
+    char howmany[4] = "All"; 
+
+    std::vector<int> select (ncv, 0);
+
+    // Arnoldi basis vectors
+    int ldz = n;
+
+    std::vector<double> 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<double> 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<eigenvectors.size(); ++i)
+      for (unsigned int j=0; j<n; ++j)
+        eigenvectors[i](j) = v[i*n+j];
+
+    delete[] workd;
+  }
+}
+
+template <typename VECTOR, typename MATRIX, typename INVERSE>
+void ArpackSolver::solve (
+    const MATRIX &system_matrix,
+    const INVERSE& inverse,
+    std::vector<double> &eigenvalues_real, 
+    std::vector<double> &eigenvalues_im, 
+    std::vector<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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.