Subversion Repositories slepc-dev

Rev

Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
6 dsic.upv.es!jroman 1 
 
2 
/*                      
3 
       This implements the power iteration for finding the eigenpair
4 
       corresponding to the eigenvalue with largest magnitude.
5 
*/
6 
#include "src/eps/epsimpl.h"
7 
 
8 
#undef __FUNCT__  
9 
#define __FUNCT__ "EPSSetUp_POWER"
10 
static int EPSSetUp_POWER(EPS eps)
11 
{
12 
  int      ierr, N;
13 
 
14 
  PetscFunctionBegin;
15 
  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
16 
  if (eps->ncv) {
17 
    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");
18 
  }
19 
  else eps->ncv = eps->nev;
20 
  if (!eps->max_it) eps->max_it = PetscMax(2000,100*N);
21 
  if (!eps->tol) eps->tol = 1.e-7;
259 dsic.upv.es!antodo 22 
 
23 
  if (eps->which!=EPS_LARGEST_MAGNITUDE)
24 
    SETERRQ(1,"Wrong value of eps->which");
25 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
26 
  ierr = EPSDefaultGetWork(eps,3);CHKERRQ(ierr);
6 dsic.upv.es!jroman 27 
  PetscFunctionReturn(0);
28 
}
29 
 
30 
#undef __FUNCT__  
31 
#define __FUNCT__ "EPSSolve_POWER"
78 dsic.upv.es!antodo 32 
static int  EPSSolve_POWER(EPS eps)
6 dsic.upv.es!jroman 33 
{
267 dsic.upv.es!antodo 34 
  int         ierr, i, maxit=eps->max_it;
6 dsic.upv.es!jroman 35 
  Vec         v, w, y, e;
36 
  PetscReal   relerr, norm, tol=eps->tol;
37 
  PetscScalar theta, alpha, eta;
38 
  PetscTruth  isSinv;
39 
 
40 
  PetscFunctionBegin;
41 
  v = eps->V[0];
42 
  y = eps->work[0];
43 
  w = eps->work[1];
44 
  e = eps->work[2];
45 
 
46 
  ierr = PetscTypeCompare((PetscObject)eps->OP,STSINV,&isSinv);CHKERRQ(ierr);
47 
 
48 
  ierr = VecCopy(eps->vec_initial,y);CHKERRQ(ierr);
49 
 
252 dsic.upv.es!jroman 50 
  eps->nconv = 0;
51 
  eps->its = 0;
6 dsic.upv.es!jroman 52 
 
281 dsic.upv.es!antodo 53 
  for (i=0;i<eps->ncv;i++) eps->eigi[i]=0.0;
54 
 
252 dsic.upv.es!jroman 55 
  while (eps->its<maxit) {
56 
 
57 
    eps->its = eps->its + 1;
58 
 
316 dsic.upv.es!antodo 59 
    if (isSinv && eps->isgeneralized) {
6 dsic.upv.es!jroman 60 
      /* w = B y */
61 
      ierr = STApplyB(eps->OP,y,w);CHKERRQ(ierr);
62 
 
63 
      /* eta = ||y||_B */
64 
      ierr = VecDot(w,y,&eta);CHKERRQ(ierr);
65 
#if !defined(PETSC_USE_COMPLEX)
66 
      if (eta<0.0) SETERRQ(1,"Negative value of eta");
67 
#endif
68 
      eta = PetscSqrtScalar(eta);
69 
 
70 
      /* normalize y and w */
71 
      ierr = VecCopy(y,v);CHKERRQ(ierr);
72 
      if (eta==0.0) SETERRQ(1,"Zero value of eta");
73 
      alpha = 1.0/eta;
74 
      ierr = VecScale(&alpha,v);CHKERRQ(ierr);
75 
      ierr = VecScale(&alpha,w);CHKERRQ(ierr);
76 
 
77 
      /* y = OP w */
78 
      ierr = STApplyNoB(eps->OP,w,y);CHKERRQ(ierr);
79 
 
252 dsic.upv.es!jroman 80 
      /* deflation of converged eigenvectors */
379 dsic.upv.es!antodo 81 
      ierr = EPSPurge(eps,y);
6 dsic.upv.es!jroman 82 
 
83 
      /* theta = w^* y */
84 
      ierr = VecDot(y,w,&theta);CHKERRQ(ierr);
85 
    }
86 
    else {
87 
      /* v = y/||y||_2 */
88 
      ierr = VecCopy(y,v);CHKERRQ(ierr);
89 
      ierr = VecNorm(y,NORM_2,&norm);CHKERRQ(ierr);
90 
      alpha = 1.0/norm;
91 
      ierr = VecScale(&alpha,v);CHKERRQ(ierr);
92 
 
93 
      /* y = OP v */
94 
      ierr = STApply(eps->OP,v,y);CHKERRQ(ierr);
95 
 
252 dsic.upv.es!jroman 96 
      /* deflation of converged eigenvectors */
379 dsic.upv.es!antodo 97 
      ierr = EPSPurge(eps,y);
6 dsic.upv.es!jroman 98 
 
99 
      /* theta = v^* y */
100 
      ierr = VecDot(y,v,&theta);CHKERRQ(ierr);
101 
    }
102 
 
252 dsic.upv.es!jroman 103 
    /* if ||y-theta v||_2 / |theta| < tol, accept */
6 dsic.upv.es!jroman 104 
    ierr = VecCopy(y,e);CHKERRQ(ierr);
105 
    alpha = -theta;
106 
    ierr = VecAXPY(&alpha,v,e);CHKERRQ(ierr);
107 
    ierr = VecNorm(e,NORM_2,&relerr);CHKERRQ(ierr);
108 
    relerr = relerr / PetscAbsScalar(theta);
109 
    eps->errest[eps->nconv] = relerr;
110 
    eps->eigr[eps->nconv] = theta;
252 dsic.upv.es!jroman 111 
 
6 dsic.upv.es!jroman 112 
    if (relerr<tol) {
316 dsic.upv.es!antodo 113 
      alpha = 1.0/norm;
114 
      ierr = VecScale(&alpha,v);CHKERRQ(ierr);
6 dsic.upv.es!jroman 115 
      eps->nconv = eps->nconv + 1;
116 
      if (eps->nconv==eps->nev) break;
117 
      v = eps->V[eps->nconv];
118 
    }
119 
 
281 dsic.upv.es!antodo 120 
    EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,eps->nconv+1);
6 dsic.upv.es!jroman 121 
  }
122 
 
123 
  if( eps->nconv == eps->nev ) eps->reason = EPS_CONVERGED_TOL;
124 
  else eps->reason = EPS_DIVERGED_ITS;
125 
 
126 
  PetscFunctionReturn(0);
127 
}
128 
 
129 
EXTERN_C_BEGIN
130 
#undef __FUNCT__  
131 
#define __FUNCT__ "EPSCreate_POWER"
132 
int EPSCreate_POWER(EPS eps)
133 
{
134 
  PetscFunctionBegin;
135 
  eps->data                      = (void *) 0;
307 dsic.upv.es!antodo 136 
  eps->ops->setfromoptions       = 0;
6 dsic.upv.es!jroman 137 
  eps->ops->setup                = EPSSetUp_POWER;
138 
  eps->ops->solve                = EPSSolve_POWER;
259 dsic.upv.es!antodo 139 
  eps->ops->destroy              = EPSDestroy_Default;
6 dsic.upv.es!jroman 140 
  eps->ops->view                 = 0;
176 dsic.upv.es!antodo 141 
  eps->ops->backtransform        = EPSBackTransform_Default;
6 dsic.upv.es!jroman 142 
  PetscFunctionReturn(0);
143 
}
144 
EXTERN_C_END