include/deal.II/lac/arpack_solver.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: arpack_solver.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //    Authors: Bärbel Janssen, University of Heidelberg,
00004 //    Agnieszka Miedler, TU Berlin, 2010
00005 //
00006 //    Copyright (C) 2010, 2012 by the deal.II authors
00007 //
00008 //    This file is subject to QPL and may not be  distributed
00009 //    without copyright and license information. Please refer
00010 //    to the file deal.II/doc/license.html for the  text  and
00011 //    further information on this license.
00012 //
00013 //---------------------------------------------------------------------------
00014 #ifndef __deal2__arpack_solver_h
00015 #define __deal2__arpack_solver_h
00016 
00017 #include <deal.II/base/config.h>
00018 
00019 #ifdef DEAL_II_USE_ARPACK
00020 
00021 DEAL_II_NAMESPACE_OPEN
00022 
00059 extern "C" void dnaupd_(int *ido, char *bmat, const unsigned int *n, char *which,
00060                         const unsigned int *nev, const double *tol, double *resid, int *ncv,
00061                         double *v, int *ldv, int *iparam, int *ipntr,
00062                         double *workd, double *workl, int *lworkl,
00063                         int *info);
00064 
00065 extern "C" void dneupd_(int *rvec, char *howmany, int *select, double *d,
00066                         double *di, double *z, int *ldz, double *sigmar,
00067                         double *sigmai, double* workev, char *bmat,const unsigned int *n, char *which,
00068                         const unsigned int *nev, const double *tol, double *resid, int *ncv,
00069                         double *v, int *ldv, int *iparam, int *ipntr,
00070                         double *workd, double *workl, int *lworkl, int *info);
00071 
00076 class ArpackSolver : public Subscriptor
00077 {
00078   public:
00084     enum WhichEigenvalues
00085     {
00086       algebraically_largest,
00087       algebraically_smallest,
00088       largest_magnitude,
00089       smallest_magnitude,
00090       largest_real_part,
00091       smallest_real_part,
00092       largest_imaginary_part,
00093       smallest_imaginary_part,
00094       both_ends
00095     };
00096 
00102     struct AdditionalData
00103     {
00104       const unsigned int number_of_arnoldi_vectors;
00105       const WhichEigenvalues eigenvalue_of_interest;
00106       AdditionalData(
00107           const unsigned int number_of_arnoldi_vectors = 15,
00108           const WhichEigenvalues eigenvalue_of_interest = largest_magnitude);
00109     };
00110 
00115     SolverControl &control () const;
00116 
00120     ArpackSolver(SolverControl& control,
00121        const AdditionalData& data = AdditionalData());
00122 
00128     template <typename VECTOR, typename MATRIX1,
00129              typename MATRIX2, typename INVERSE>
00130     void solve(
00131         const MATRIX1 &A,
00132         const MATRIX2 &B,
00133         const INVERSE& inverse,
00134         std::vector<std::complex<double> > &eigenvalues,
00135         std::vector<VECTOR> &eigenvectors,
00136         const unsigned int n_eigenvalues);
00137 
00138   protected:
00139 
00145     SolverControl &solver_control;
00146 
00151     const AdditionalData additional_data;
00152 
00153   private:
00154 
00158     DeclException2 (ExcInvalidNumberofEigenvalues, int, int,
00159         << "Number of wanted eigenvalues " << arg1
00160         << " is larger that the size of the matrix " << arg2);
00161 
00162     DeclException2 (ExcInvalidNumberofArnoldiVectors, int, int,
00163         << "Number of Arnoldi vectors " << arg1
00164         << " is larger that the size of the matrix " << arg2);
00165 
00166     DeclException2 (ExcSmallNumberofArnoldiVectors, int, int,
00167         << "Number of Arnoldi vectors " << arg1
00168         << " is too small to obtain " << arg2
00169         << " eigenvalues");
00170 
00171     DeclException1 (ExcArpackIdo, int, << "This ido " << arg1
00172         << " is not supported. Check documentation of ARPACK");
00173 
00174     DeclException1 (ExcArpackMode, int, << "This mode " << arg1
00175         << " is not supported. Check documentation of ARPACK");
00176 
00177     DeclException1 (ExcArpackInfodsaupd, int,
00178         << "Error with dsaupd, info " << arg1
00179         << ". Check documentation of ARPACK");
00180 
00181     DeclException1 (ExcArpackInfodneupd, int,
00182         << "Error with dneupd, info " << arg1
00183         << ". Check documentation of ARPACK");
00184 
00185     DeclException1 (ExcArpackInfoMaxIt, int,
00186         << "Maximum number " << arg1
00187         << " of iterations reached.");
00188 
00189     DeclException1 (ExcArpackNoShifts, int,
00190         << "No shifts could be applied during implicit"
00191         << " Arnoldi update, try increasing the number of"
00192         << " Arnoldi vectors.");
00193 };
00194 
00195 ArpackSolver::AdditionalData::
00196 AdditionalData (const unsigned int number_of_arnoldi_vectors,
00197     const WhichEigenvalues eigenvalue_of_interest)
00198   :
00199     number_of_arnoldi_vectors(number_of_arnoldi_vectors),
00200     eigenvalue_of_interest(eigenvalue_of_interest)
00201 {}
00202 
00203 ArpackSolver::ArpackSolver (SolverControl& control,
00204     const AdditionalData& data)
00205     :
00206       solver_control (control),
00207       additional_data (data)
00208 
00209 {}
00210 
00211 template <typename VECTOR, typename MATRIX1,
00212          typename MATRIX2, typename INVERSE>
00213 void ArpackSolver::solve (
00214     const MATRIX1 &system_matrix,
00215     const MATRIX2 &mass_matrix,
00216     const INVERSE& inverse,
00217     std::vector<std::complex<double> > &eigenvalues,
00218     std::vector<VECTOR> &eigenvectors,
00219     const unsigned int n_eigenvalues)
00220 {
00221   //inside the routines of ARPACK the
00222   //values change magically, so store
00223   //them here
00224 
00225   const unsigned int n = system_matrix.m();
00226   const unsigned int n_inside_arpack = system_matrix.m();
00227 
00228    /*
00229    if(n < 0 || nev <0 || p < 0 || maxit < 0 )
00230         std:cout << "All input parameters have to be positive.\n";
00231         */
00232   Assert (n_eigenvalues < n,
00233         ExcInvalidNumberofEigenvalues(n_eigenvalues, n));
00234 
00235   Assert (additional_data.number_of_arnoldi_vectors < n,
00236         ExcInvalidNumberofArnoldiVectors(
00237           additional_data.number_of_arnoldi_vectors, n));
00238 
00239   Assert (additional_data.number_of_arnoldi_vectors > 2*n_eigenvalues+1,
00240         ExcSmallNumberofArnoldiVectors(
00241           additional_data.number_of_arnoldi_vectors, n_eigenvalues));
00242   // ARPACK mode for dnaupd, here only mode 3
00243   int mode = 3;
00244 
00245   // reverse communication parameter
00246   int ido = 0;
00247 
00252   char bmat[2] = "G";
00253 
00266   char which[3];
00267   switch(additional_data.eigenvalue_of_interest)
00268   {
00269     case algebraically_largest:
00270       strcpy (which, "LA");
00271       break;
00272     case algebraically_smallest:
00273       strcpy (which, "SA");
00274       break;
00275     case largest_magnitude:
00276       strcpy (which, "LM");
00277       break;
00278     case smallest_magnitude:
00279       strcpy (which, "SM");
00280       break;
00281     case largest_real_part:
00282       strcpy (which, "LR");
00283       break;
00284     case smallest_real_part:
00285       strcpy (which, "SR");
00286       break;
00287     case largest_imaginary_part:
00288       strcpy (which, "LI");
00289       break;
00290     case smallest_imaginary_part:
00291       strcpy (which, "SI");
00292       break;
00293     case both_ends:
00294       strcpy (which, "BE");
00295       break;
00296   }
00297 
00298   // tolerance for ARPACK
00299   const double tol = control().tolerance();
00300 
00301   // if the starting vector is used it has to be in resid
00302   std::vector<double> resid(n, 1.);
00303 
00304   // number of Arnoldi basis vectors specified
00305   // in additional_data
00306   int ncv = additional_data.number_of_arnoldi_vectors;
00307 
00308   int ldv = n;
00309   std::vector<double> v (ldv*ncv, 0.0);
00310 
00311   //information to the routines
00312   std::vector<int> iparam (11, 0);
00313 
00314   iparam[0] = 1;        // shift strategy
00315 
00316   // maximum number of iterations
00317   iparam[2] = control().max_steps();
00318 
00327   iparam[6] = mode;
00328   std::vector<int> ipntr (14, 0);
00329 
00330   // work arrays for ARPACK
00331   double *workd;
00332   workd = new double[3*n];
00333 
00334   for(unsigned int i=0; i<3*n; ++i)
00335     workd[i] = 0.0;
00336 
00337   int lworkl = 3*ncv*(ncv + 6);
00338   std::vector<double> workl (lworkl, 0.);
00339   //information out of the iteration
00340   int info = 1;
00341 
00342   const unsigned int nev = n_eigenvalues;
00343   while(ido != 99)
00344   {
00345     // call of ARPACK dnaupd routine
00346     dnaupd_(&ido, bmat, &n_inside_arpack, which, &nev, &tol,
00347         &resid[0], &ncv, &v[0], &ldv, &iparam[0], &ipntr[0],
00348         workd, &workl[0], &lworkl, &info);
00349 
00350     if(ido == 99)
00351       break;
00352 
00353     switch(mode)
00354     {
00355       case 3:
00356         {
00357           switch(ido)
00358           {
00359             case -1:
00360               {
00361 
00362                 VECTOR src,dst,tmp;
00363                 src.reinit(eigenvectors[0]);
00364                 dst.reinit(src);
00365                 tmp.reinit(src);
00366 
00367 
00368                 for(unsigned int i=0; i<src.size(); ++i)
00369                   src(i) = *(workd+ipntr[0]-1+i);
00370 
00371                 // multiplication with mass matrix M
00372                 mass_matrix.vmult(tmp, src);
00373                 // solving linear system
00374                 inverse.vmult(dst,tmp);
00375 
00376                 for(unsigned int i=0; i<dst.size(); ++i)
00377                   *(workd+ipntr[1]-1+i) = dst(i);
00378               }
00379               break;
00380 
00381             case  1:
00382               {
00383 
00384                 VECTOR src,dst,tmp, tmp2;
00385                 src.reinit(eigenvectors[0]);
00386                 dst.reinit(src);
00387                 tmp.reinit(src);
00388                 tmp2.reinit(src);
00389 
00390                 for(unsigned int i=0; i<src.size(); ++i)
00391                 {
00392                   src(i) = *(workd+ipntr[2]-1+i);
00393                   tmp(i) = *(workd+ipntr[0]-1+i);
00394                 }
00395                 // solving linear system
00396                 inverse.vmult(dst,src);
00397 
00398                 for(unsigned int i=0; i<dst.size(); ++i)
00399                   *(workd+ipntr[1]-1+i) = dst(i);
00400               }
00401               break;
00402 
00403             case  2:
00404               {
00405 
00406                 VECTOR src,dst;
00407                 src.reinit(eigenvectors[0]);
00408                 dst.reinit(src);
00409 
00410                 for(unsigned int i=0; i<src.size(); ++i)
00411                   src(i) = *(workd+ipntr[0]-1+i);
00412 
00413                 // Multiplication with mass matrix M
00414                 mass_matrix.vmult(dst, src);
00415 
00416                 for(unsigned int i=0; i<dst.size(); ++i)
00417                   *(workd+ipntr[1]-1+i) = dst(i);
00418 
00419               }
00420               break;
00421 
00422             default:
00423               Assert (false, ExcArpackIdo(ido));
00424               break;
00425           }
00426         }
00427         break;
00428       default:
00429         Assert (false, ExcArpackMode(mode));
00430         break;
00431     }
00432   }
00433 
00434   if (info<0)
00435   {
00436     Assert (false, ExcArpackInfodsaupd(info));
00437   }
00438   else
00439   {
00443     int rvec = 1;
00444 
00445     // which eigenvectors
00446     char howmany[4] = "All";
00447 
00448     std::vector<int> select (ncv, 0);
00449 
00450     int ldz = n;
00451 
00452     std::vector<double> z (ldz*ncv, 0.);
00453 
00454     double sigmar = 0.0; // real part of the shift
00455     double sigmai = 0.0; // imaginary part of the shift
00456 
00457     int lworkev = 3*ncv;
00458     std::vector<double> workev (lworkev, 0.);
00459 
00460     std::vector<double> eigenvalues_real (n_eigenvalues, 0.);
00461     std::vector<double> eigenvalues_im (n_eigenvalues, 0.);
00462 
00463     // call of ARPACK dneupd routine
00464     dneupd_(&rvec, howmany, &select[0], &eigenvalues_real[0],
00465         &eigenvalues_im[0], &z[0], &ldz, &sigmar, &sigmai,
00466         &workev[0], bmat, &n_inside_arpack, which, &nev, &tol,
00467         &resid[0], &ncv, &v[0], &ldv,
00468         &iparam[0], &ipntr[0], workd, &workl[0], &lworkl, &info);
00469 
00470     if (info == 1)
00471     {
00472       Assert (false, ExcArpackInfoMaxIt(control().max_steps()));
00473     }
00474     else if (info == 3)
00475     {
00476       Assert (false, ExcArpackNoShifts(1));;
00477     }
00478     else if (info!=0)
00479     {
00480       Assert (false, ExcArpackInfodneupd(info));
00481     }
00482 
00483     for (unsigned int i=0; i<eigenvectors.size(); ++i)
00484       for (unsigned int j=0; j<n; ++j)
00485         eigenvectors[i](j) = v[i*n+j];
00486 
00487     delete[] workd;
00488 
00489     AssertDimension (eigenvalues.size(), eigenvalues_real.size());
00490     AssertDimension (eigenvalues.size(), eigenvalues_im.size());
00491 
00492     for(unsigned int i=0; i<eigenvalues.size(); ++i)
00493     {
00494       eigenvalues[i].real() = eigenvalues_real[i];
00495       eigenvalues[i].imag() = eigenvalues_im[i];
00496     }
00497   }
00498 }
00499 
00500 SolverControl& ArpackSolver::control () const
00501 {
00502   return solver_control;
00503 }
00504 
00505 DEAL_II_NAMESPACE_CLOSE
00506 
00507 #endif
00508 #endif
00509 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:39 by doxygen 1.7.3