00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
00222
00223
00224
00225 const unsigned int n = system_matrix.m();
00226 const unsigned int n_inside_arpack = system_matrix.m();
00227
00228
00229
00230
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
00243 int mode = 3;
00244
00245
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
00299 const double tol = control().tolerance();
00300
00301
00302 std::vector<double> resid(n, 1.);
00303
00304
00305
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
00312 std::vector<int> iparam (11, 0);
00313
00314 iparam[0] = 1;
00315
00316
00317 iparam[2] = control().max_steps();
00318
00327 iparam[6] = mode;
00328 std::vector<int> ipntr (14, 0);
00329
00330
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
00340 int info = 1;
00341
00342 const unsigned int nev = n_eigenvalues;
00343 while(ido != 99)
00344 {
00345
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
00372 mass_matrix.vmult(tmp, src);
00373
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
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
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
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;
00455 double sigmai = 0.0;
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
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