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
647 dsic.upv.es!antodo 1 
/*                      
2 
 
3 
   SLEPc eigensolver: "lanczos"
4 
 
655 dsic.upv.es!jroman 5 
   Method: Explicitly Restarted Symmetric/Hermitian Lanczos
647 dsic.upv.es!antodo 6 
 
7 
   Description:
8 
 
655 dsic.upv.es!jroman 9 
       This solver implements the Lanczos method for symmetric (Hermitian)
10 
       problems, with explicit restart and deflation. When building the
11 
       Lanczos factorization, several reorthogonalization strategies can
12 
       be selected.
13 
 
647 dsic.upv.es!antodo 14 
   Algorithm:
15 
 
655 dsic.upv.es!jroman 16 
       The implemented algorithm builds a Lanczos factorization of order
17 
       ncv. Converged eigenpairs are locked and the iteration is restarted
18 
       with the rest of the columns being the active columns for the next
19 
       Lanczos factorization. Currently, no filtering is applied to the
20 
       vector used for restarting.
21 
 
22 
       The following reorthogonalization schemes are currently implemented:
23 
 
24 
       - Full reorthogonalization: at each Lanczos step, the corresponding
25 
       Lanczos vector is orthogonalized with respect to all the previous
26 
       vectors.
27 
 
647 dsic.upv.es!antodo 28 
   References:
29 
 
655 dsic.upv.es!jroman 30 
       [1] B.N. Parlett, "The Symmetric Eigenvalue Problem", SIAM Classics in
31 
       Applied Mathematics (1998), ch. 13.
647 dsic.upv.es!antodo 32 
 
655 dsic.upv.es!jroman 33 
       [2] L. Komzsik, "The Lanczos Method. Evolution and Application", SIAM
34 
       (2003).
35 
 
36 
       [3] B.N. Parlett and D.S. Scott, The Lanczos algorithm with selective
37 
       orthogonalization, Math. Comp., 33 (1979), no. 145, 217-238.
38 
 
39 
       [4] H.D. Simon, The Lanczos algorithm with partial reorthogonalization,
40 
       Math. Comp., 42 (1984), no. 165, 115-142.
41 
 
42 
   Last update: October 2004
43 
 
647 dsic.upv.es!antodo 44 
*/
707 dsic.upv.es!antodo 45 
#include "src/eps/epsimpl.h"                /*I "slepceps.h" I*/
647 dsic.upv.es!antodo 46 
#include "slepcblaslapack.h"
47 
 
671 dsic.upv.es!antodo 48 
typedef struct {
906 dsic.upv.es!antodo 49 
  EPSLanczosReorthogType reorthog;
671 dsic.upv.es!antodo 50 
} EPS_LANCZOS;
51 
 
647 dsic.upv.es!antodo 52 
#undef __FUNCT__  
53 
#define __FUNCT__ "EPSSetUp_LANCZOS"
54 
PetscErrorCode EPSSetUp_LANCZOS(EPS eps)
55 
{
56 
  PetscErrorCode ierr;
57 
  int            N;
58 
 
59 
  PetscFunctionBegin;
60 
  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
651 dsic.upv.es!antodo 61 
  if (eps->nev > N) eps->nev = N;
647 dsic.upv.es!antodo 62 
  if (eps->ncv) {
63 
    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");
64 
  }
651 dsic.upv.es!antodo 65 
  else eps->ncv = PetscMin(N,PetscMax(2*eps->nev,eps->nev+15));
647 dsic.upv.es!antodo 66 
  if (!eps->max_it) eps->max_it = PetscMax(100,N);
67 
  if (!eps->tol) eps->tol = 1.e-7;
906 dsic.upv.es!antodo 68 
  if (eps->solverclass==EPS_ONE_SIDE) {
69 
    if (eps->which == EPS_LARGEST_IMAGINARY || eps->which == EPS_SMALLEST_IMAGINARY)
70 
      SETERRQ(1,"Wrong value of eps->which");
71 
    if (!eps->ishermitian)
72 
      SETERRQ(PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
73 
  } else {
74 
    if (eps->which != EPS_LARGEST_MAGNITUDE)
75 
      SETERRQ(1,"Wrong value of eps->which");
76 
  }
647 dsic.upv.es!antodo 77 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
78 
  if (eps->T) { ierr = PetscFree(eps->T);CHKERRQ(ierr); }  
79 
  ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&eps->T);CHKERRQ(ierr);
885 dsic.upv.es!jroman 80 
  if (eps->solverclass==EPS_TWO_SIDE) {
81 
    if (eps->Tl) { ierr = PetscFree(eps->Tl);CHKERRQ(ierr); }  
82 
    ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&eps->Tl);CHKERRQ(ierr);
83 
    ierr = EPSDefaultGetWork(eps,2);CHKERRQ(ierr);
84 
  }
85 
  else { ierr = EPSDefaultGetWork(eps,1);CHKERRQ(ierr); }
647 dsic.upv.es!antodo 86 
  PetscFunctionReturn(0);
87 
}
88 
 
89 
#undef __FUNCT__  
655 dsic.upv.es!jroman 90 
#define __FUNCT__ "EPSFullLanczos"
647 dsic.upv.es!antodo 91 
/*
655 dsic.upv.es!jroman 92 
   EPSFullLanczos - Full reorthogonalization.
93 
 
94 
   In this variant, at each Lanczos step, the corresponding Lanczos vector
95 
   is orthogonalized with respect to all the previous Lanczos vectors.
647 dsic.upv.es!antodo 96 
*/
891 dsic.upv.es!antodo 97 
static PetscErrorCode EPSFullLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta,PetscTruth *breakdown)
647 dsic.upv.es!antodo 98 
{
99 
  PetscErrorCode ierr;
772 dsic.upv.es!antodo 100 
  int            j,m = *M;
647 dsic.upv.es!antodo 101 
  PetscReal      norm;
102 
 
103 
  PetscFunctionBegin;
891 dsic.upv.es!antodo 104 
  for (j=k;j<m;j++) {
105 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
704 dsic.upv.es!antodo 106 
    eps->its++;
891 dsic.upv.es!antodo 107 
    ierr = EPSOrthogonalize(eps,eps->nds,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
108 
    ierr = EPSOrthogonalize(eps,j+1,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);
109 
    if (*breakdown) {
110 
      *M = j+1;
111 
      break;
112 
    }
113 
    if (j<m-1) {
772 dsic.upv.es!antodo 114 
      T[m*j+j+1] = norm;
891 dsic.upv.es!antodo 115 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
116 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
647 dsic.upv.es!antodo 117 
    }
118 
  }
891 dsic.upv.es!antodo 119 
  *beta = norm;
647 dsic.upv.es!antodo 120 
  PetscFunctionReturn(0);
121 
}
122 
 
655 dsic.upv.es!jroman 123 
#undef __FUNCT__  
891 dsic.upv.es!antodo 124 
#define __FUNCT__ "EPSSimpleLanczos"
683 dsic.upv.es!jroman 125 
/*
891 dsic.upv.es!antodo 126 
   EPSSimpleLanczos - Local reorthogonalization.
683 dsic.upv.es!jroman 127 
 
128 
   This is the simplest variant. At each Lanczos step, the corresponding Lanczos vector
129 
   is orthogonalized with respect to the two previous Lanczos vectors, according to
130 
   the three term Lanczos recurrence. WARNING: This variant does not track the loss of
131 
   orthogonality that occurs in finite-precision arithmetic and, therefore, the
132 
   generated vectors are not guaranteed to be (semi-)orthogonal.
133 
*/
891 dsic.upv.es!antodo 134 
static PetscErrorCode EPSSimpleLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta,PetscTruth *breakdown)
671 dsic.upv.es!antodo 135 
{
136 
  PetscErrorCode ierr;
772 dsic.upv.es!antodo 137 
  int            j,m = *M;
671 dsic.upv.es!antodo 138 
  PetscReal      norm;
139 
 
140 
  PetscFunctionBegin;
891 dsic.upv.es!antodo 141 
  for (j=k;j<m;j++) {
142 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
704 dsic.upv.es!antodo 143 
    eps->its++;
897 dsic.upv.es!antodo 144 
    ierr = EPSOrthogonalize(eps,eps->nds+k,eps->DSV,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
891 dsic.upv.es!antodo 145 
    if (j == k) {
146 
      ierr = EPSOrthogonalize(eps,1,V+j,f,T+m*j+j,&norm,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 147 
    } else {
891 dsic.upv.es!antodo 148 
      ierr = EPSOrthogonalize(eps,2,V+j-1,f,T+m*j+j-1,&norm,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 149 
    }
891 dsic.upv.es!antodo 150 
    if (*breakdown) {
151 
      *M = j+1;
152 
      break;
153 
    }
154 
    if (j<m-1) {
155 
      T[m*j+j+1] = norm;
156 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
157 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
158 
    }
159 
  }
160 
  *beta = norm;
161 
  PetscFunctionReturn(0);
162 
}
163 
 
164 
#undef __FUNCT__  
165 
#define __FUNCT__ "EPSSelectiveLanczos"
166 
static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta,PetscTruth *breakdown,PetscReal anorm)
167 
{
902 dsic.upv.es!antodo 168 
#if defined(SLEPC_MISSING_LAPACK_DSTEVR) || defined(SLEPC_MISSING_LAPACK_DLAMCH)
169 
  PetscFunctionBegin;
170 
  SETERRQ(PETSC_ERR_SUP,"DSTEVR - Lapack routine is unavailable.");
171 
#else
891 dsic.upv.es!antodo 172 
  PetscErrorCode ierr;
895 dsic.upv.es!antodo 173 
  int            i,j,m = *M,n,vl,vu,mout,*isuppz,*iwork,lwork,liwork,info;
174 
  PetscReal      *D,*E,*ritz,*Y,*work,abstol,norm;
891 dsic.upv.es!antodo 175 
  PetscTruth     conv;
176 
 
177 
  PetscFunctionBegin;
895 dsic.upv.es!antodo 178 
  ierr = PetscMalloc(m*sizeof(PetscReal),&D);CHKERRQ(ierr);
179 
  ierr = PetscMalloc(m*sizeof(PetscReal),&E);CHKERRQ(ierr);
891 dsic.upv.es!antodo 180 
  ierr = PetscMalloc(m*sizeof(PetscReal),&ritz);CHKERRQ(ierr);
181 
  ierr = PetscMalloc(m*m*sizeof(PetscReal),&Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 182 
  ierr = PetscMalloc(2*m*sizeof(int),&isuppz);CHKERRQ(ierr);
183 
  lwork = 20*m;
184 
  ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr);
185 
  liwork = 10*m;
186 
  ierr = PetscMalloc(liwork*sizeof(int),&iwork);CHKERRQ(ierr);
187 
  abstol = 2.0*LAPACKlamch_("S",1);
891 dsic.upv.es!antodo 188 
 
189 
  for (j=k;j<m;j++) {
190 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
191 
    eps->its++;
897 dsic.upv.es!antodo 192 
    ierr = EPSOrthogonalize(eps,eps->nds+k,eps->DSV,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
891 dsic.upv.es!antodo 193 
    if (j == k) {
194 
      ierr = EPSOrthogonalize(eps,1,V+j,f,T+m*j+j,&norm,breakdown);CHKERRQ(ierr);
704 dsic.upv.es!antodo 195 
    } else {
891 dsic.upv.es!antodo 196 
      ierr = EPSOrthogonalize(eps,2,V+j-1,f,T+m*j+j-1,&norm,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 197 
    }
891 dsic.upv.es!antodo 198 
 
199 
    if (*breakdown) {
200 
      *M = j+1;
201 
      break;
202 
    }
203 
 
204 
    n = j-k+1;
205 
    for (i=0;i<n-1;i++) {
206 
      ritz[i] = PetscRealPart(T[(i+k)*(m+1)]);
207 
      E[i] = PetscRealPart(T[(i+k)*(m+1)+1]);
208 
    }
209 
    ritz[n-1] = PetscRealPart(T[(n-1+k)*(m+1)]);
210 
 
211 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
895 dsic.upv.es!antodo 212 
    LAPACKstevr_("V","A",&n,D,E,&vl,&vu,PETSC_NULL,PETSC_NULL,&abstol,&mout,ritz,Y,&n,isuppz,work,&lwork,iwork,&liwork,&info,1,1);
891 dsic.upv.es!antodo 213 
    if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xSTEQR %i",info);
214 
 
215 
    /* Estimate ||A|| */
216 
    for (i=0;i<n;i++)
217 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
218 
 
219 
    /* Exit if residual norm is small [Parlett, page 300] */
220 
    conv = PETSC_FALSE;
221 
    for (i=0;i<n && !conv;i++)
222 
      if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm)
223 
        conv = PETSC_TRUE;
224 
 
225 
    if (conv) {
226 
      *M = j+1;
227 
      break;
228 
    }
229 
 
230 
    if (j<m-1) {
231 
      T[m*j+j+1] = norm;
232 
      ierr = VecScale(f,1.0 / norm);CHKERRQ(ierr);
233 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
234 
    }
671 dsic.upv.es!antodo 235 
  }
891 dsic.upv.es!antodo 236 
  *beta = norm;
895 dsic.upv.es!antodo 237 
  ierr = PetscFree(D);CHKERRQ(ierr);
238 
  ierr = PetscFree(E);CHKERRQ(ierr);
891 dsic.upv.es!antodo 239 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 240 
  ierr = PetscFree(Y);CHKERRQ(ierr);
241 
  ierr = PetscFree(isuppz);CHKERRQ(ierr);
891 dsic.upv.es!antodo 242 
  ierr = PetscFree(work);CHKERRQ(ierr);
895 dsic.upv.es!antodo 243 
  ierr = PetscFree(iwork);CHKERRQ(ierr);
671 dsic.upv.es!antodo 244 
  PetscFunctionReturn(0);
902 dsic.upv.es!antodo 245 
#endif
671 dsic.upv.es!antodo 246 
}
247 
 
248 
#undef __FUNCT__  
891 dsic.upv.es!antodo 249 
#define __FUNCT__ "EPSPeriodicLanczos"
250 
static PetscErrorCode EPSPeriodicLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta, PetscTruth *breakdown)
772 dsic.upv.es!antodo 251 
{
252 
  PetscErrorCode ierr;
891 dsic.upv.es!antodo 253 
  int            i,j,m = *M;
254 
  PetscReal      norm;
255 
  PetscScalar    *omega;
256 
  PetscTruth     reorthog;
746 dsic.upv.es!antodo 257 
 
772 dsic.upv.es!antodo 258 
  PetscFunctionBegin;
891 dsic.upv.es!antodo 259 
  reorthog = PETSC_FALSE;
260 
  ierr = PetscMalloc(m*sizeof(PetscScalar),&omega);CHKERRQ(ierr);
261 
  for (j=k;j<m;j++) {
772 dsic.upv.es!antodo 262 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
263 
    eps->its++;
264 
    ierr = EPSOrthogonalize(eps,eps->nds,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
891 dsic.upv.es!antodo 265 
 
266 
    if (reorthog) {
267 
      ierr = EPSOrthogonalize(eps,j+1,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);
772 dsic.upv.es!antodo 268 
    } else {
891 dsic.upv.es!antodo 269 
      for (i=0;i<j-1;i++) T[m*j+i] = 0.0;
270 
      if (j == k) {
271 
        ierr = EPSOrthogonalize(eps,1,V+j,f,T+m*j+j,&norm,breakdown);CHKERRQ(ierr);
272 
      } else {
273 
        ierr = EPSOrthogonalize(eps,2,V+j-1,f,T+m*j+j-1,&norm,breakdown);CHKERRQ(ierr);
274 
      }
772 dsic.upv.es!antodo 275 
    }
891 dsic.upv.es!antodo 276 
 
277 
    if (*breakdown) {
278 
      *M = j+1;
279 
      break;
772 dsic.upv.es!antodo 280 
    }
281 
 
891 dsic.upv.es!antodo 282 
    if (reorthog) {
283 
      reorthog = PETSC_FALSE;
284 
    } else if (j>1) {
285 
      ierr = VecMDot(j-1,f,eps->V,omega);CHKERRQ(ierr);
286 
      for (i=0;i<j-1 && !reorthog;i++)
897 dsic.upv.es!antodo 287 
        if (PetscAbsScalar(omega[i]) > PETSC_SQRT_MACHINE_EPSILON/sqrt(j)) {
891 dsic.upv.es!antodo 288 
          reorthog = PETSC_TRUE;
289 
        }
290 
      if (reorthog) {
291 
        ierr = EPSOrthogonalize(eps,j-1,V,f,T+m*j,&norm,PETSC_NULL);CHKERRQ(ierr);
772 dsic.upv.es!antodo 292 
      }
293 
    }
891 dsic.upv.es!antodo 294 
 
295 
    if (j<m-1) {
296 
      T[m*j+j+1] = norm;
297 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
298 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
772 dsic.upv.es!antodo 299 
    }
300 
  }
891 dsic.upv.es!antodo 301 
  *beta = norm;
772 dsic.upv.es!antodo 302 
  ierr = PetscFree(omega);CHKERRQ(ierr);
303 
  PetscFunctionReturn(0);
304 
}
305 
 
746 dsic.upv.es!antodo 306 
#undef __FUNCT__  
891 dsic.upv.es!antodo 307 
#define __FUNCT__ "EPSPartialLanczos"
308 
static PetscErrorCode EPSPartialLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta, PetscTruth *breakdown)
746 dsic.upv.es!antodo 309 
{
310 
  PetscErrorCode ierr;
891 dsic.upv.es!antodo 311 
  int            i,j,l,m = *M;
895 dsic.upv.es!antodo 312 
  PetscReal      norm,nu;
313 
  PetscScalar    *omega;
891 dsic.upv.es!antodo 314 
  PetscTruth     reorthog,*which;
746 dsic.upv.es!antodo 315 
 
316 
  PetscFunctionBegin;
891 dsic.upv.es!antodo 317 
  nu = sqrt(PETSC_MACHINE_EPSILON*PETSC_SQRT_MACHINE_EPSILON);
746 dsic.upv.es!antodo 318 
  reorthog = PETSC_FALSE;
891 dsic.upv.es!antodo 319 
  ierr = PetscMalloc(m*sizeof(PetscScalar),&omega);CHKERRQ(ierr);
320 
  ierr = PetscMalloc(m*sizeof(PetscTruth),&which);CHKERRQ(ierr);
321 
  for (j=k;j<m;j++) {
322 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
746 dsic.upv.es!antodo 323 
    eps->its++;
891 dsic.upv.es!antodo 324 
    ierr = EPSOrthogonalize(eps,eps->nds,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
325 
 
326 
    if (j == k) {
327 
      ierr = EPSOrthogonalize(eps,1,V+j,f,T+m*j+j,&norm,breakdown);CHKERRQ(ierr);
746 dsic.upv.es!antodo 328 
    } else {
891 dsic.upv.es!antodo 329 
      ierr = EPSOrthogonalize(eps,2,V+j-1,f,T+m*j+j-1,&norm,breakdown);CHKERRQ(ierr);
746 dsic.upv.es!antodo 330 
    }
891 dsic.upv.es!antodo 331 
 
332 
    if (*breakdown) {
333 
      *M = j+1;
334 
      break;
335 
    }
336 
 
337 
    if (reorthog) {
338 
      for (i=0;i<j-1;i++)
339 
        if (which[i]) {
340 
          ierr = EPSOrthogonalize(eps,1,V+i,f,PETSC_NULL,&norm,PETSC_NULL);CHKERRQ(ierr);
341 
        }
342 
      reorthog = PETSC_FALSE;
343 
    } else if (j>1) {
344 
      ierr = VecMDot(j-1,f,eps->V,omega);CHKERRQ(ierr);
345 
      for (i=0;i<j-1;i++) {
346 
        omega[i] /= norm;
347 
        which[i] = PETSC_FALSE;
746 dsic.upv.es!antodo 348 
      }
891 dsic.upv.es!antodo 349 
      for (i=0;i<j-1;i++)
350 
        if (PetscAbsScalar(omega[i]) > PETSC_SQRT_MACHINE_EPSILON) {
351 
          reorthog = PETSC_TRUE;
352 
          which[i] = PETSC_TRUE;
895 dsic.upv.es!antodo 353 
          for (l=i-1;l>0 && PetscAbsScalar(omega[l]) > nu;l--) which[l] = PETSC_TRUE;
354 
          for (l=i+1;l<j-1 && PetscAbsScalar(omega[l]) > nu;l++) which[l] = PETSC_TRUE;
891 dsic.upv.es!antodo 355 
        }
356 
      if (reorthog)
357 
        for (i=0;i<j-1;i++)
358 
          if (which[i]) {
359 
            ierr = EPSOrthogonalize(eps,1,V+i,f,PETSC_NULL,&norm,PETSC_NULL);CHKERRQ(ierr);
360 
          }
746 dsic.upv.es!antodo 361 
    }
891 dsic.upv.es!antodo 362 
 
363 
    if (j<m-1) {
746 dsic.upv.es!antodo 364 
      T[m*j+j+1] = norm;
891 dsic.upv.es!antodo 365 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
366 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
746 dsic.upv.es!antodo 367 
    }
368 
  }
891 dsic.upv.es!antodo 369 
  *beta = norm;
746 dsic.upv.es!antodo 370 
  ierr = PetscFree(omega);CHKERRQ(ierr);
891 dsic.upv.es!antodo 371 
  ierr = PetscFree(which);CHKERRQ(ierr);
746 dsic.upv.es!antodo 372 
  PetscFunctionReturn(0);
373 
}
374 
 
375 
#undef __FUNCT__  
655 dsic.upv.es!jroman 376 
#define __FUNCT__ "EPSBasicLanczos"
377 
/*
378 
   EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
379 
   columns are assumed to be locked and therefore they are not modified. On
380 
   exit, the following relation is satisfied:
647 dsic.upv.es!antodo 381 
 
655 dsic.upv.es!jroman 382 
                    OP * V - V * T = f * e_m^T
383 
 
384 
   where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
385 
   f is the residual vector and e_m is the m-th vector of the canonical basis.
386 
   The Lanczos vectors (together with vector f) are B-orthogonal (to working
387 
   accuracy) if full reorthogonalization is being used, otherwise they are
388 
   (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
389 
   Lanczos vector can be computed as v_{m+1} = f / beta.
390 
 
391 
   This function simply calls another function which depends on the selected
392 
   reorthogonalization strategy.
393 
*/
891 dsic.upv.es!antodo 394 
static PetscErrorCode EPSBasicLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *m,Vec f,PetscReal *beta,PetscTruth *breakdown,PetscReal anorm)
655 dsic.upv.es!jroman 395 
{
671 dsic.upv.es!antodo 396 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
655 dsic.upv.es!jroman 397 
  PetscErrorCode ierr;
398 
 
399 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 400 
  switch (lanczos->reorthog) {
401 
    case EPSLANCZOS_ORTHOG_NONE:
891 dsic.upv.es!antodo 402 
      ierr = EPSSimpleLanczos(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 403 
      break;
772 dsic.upv.es!antodo 404 
    case EPSLANCZOS_ORTHOG_SELECTIVE:
891 dsic.upv.es!antodo 405 
      ierr = EPSSelectiveLanczos(eps,T,V,k,m,f,beta,breakdown,anorm);CHKERRQ(ierr);
671 dsic.upv.es!antodo 406 
      break;
891 dsic.upv.es!antodo 407 
    case EPSLANCZOS_ORTHOG_PARTIAL:
408 
      ierr = EPSPartialLanczos(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
409 
      break;
746 dsic.upv.es!antodo 410 
    case EPSLANCZOS_ORTHOG_PERIODIC:
891 dsic.upv.es!antodo 411 
      ierr = EPSPeriodicLanczos(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
746 dsic.upv.es!antodo 412 
      break;
772 dsic.upv.es!antodo 413 
    case EPSLANCZOS_ORTHOG_FULL:
891 dsic.upv.es!antodo 414 
      ierr = EPSFullLanczos(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
772 dsic.upv.es!antodo 415 
      break;
671 dsic.upv.es!antodo 416 
  }
655 dsic.upv.es!jroman 417 
  PetscFunctionReturn(0);
418 
}
419 
 
647 dsic.upv.es!antodo 420 
#undef __FUNCT__  
421 
#define __FUNCT__ "EPSSolve_LANCZOS"
422 
PetscErrorCode EPSSolve_LANCZOS(EPS eps)
423 
{
902 dsic.upv.es!antodo 424 
#if defined(SLEPC_MISSING_LAPACK_DSTEVR) || defined(SLEPC_MISSING_LAPACK_DLAMCH)
817 dsic.upv.es!antodo 425 
  PetscFunctionBegin;
902 dsic.upv.es!antodo 426 
  SETERRQ(PETSC_ERR_SUP,"DSTEVR - Lapack routine is unavailable.");
806 dsic.upv.es!antodo 427 
#else
891 dsic.upv.es!antodo 428 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
647 dsic.upv.es!antodo 429 
  PetscErrorCode ierr;
902 dsic.upv.es!antodo 430 
  int            nconv,i,j,k,n,m,*perm,restart,
895 dsic.upv.es!antodo 431 
                 *isuppz,*iwork,mout,lwork,liwork,vl,vu,info,
746 dsic.upv.es!antodo 432 
                 ncv=eps->ncv;
870 dsic.upv.es!antodo 433 
  Vec            f=eps->work[0];
902 dsic.upv.es!antodo 434 
  PetscScalar    *T=eps->T;
435 
  PetscReal      *ritz,*Y,*bnd,*D,*E,*work,anorm,beta,norm,abstol,restart_ritz;
897 dsic.upv.es!antodo 436 
  PetscTruth     breakdown;
437 
  char           *conv;
647 dsic.upv.es!antodo 438 
 
439 
  PetscFunctionBegin;
895 dsic.upv.es!antodo 440 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&D);CHKERRQ(ierr);
441 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&E);CHKERRQ(ierr);
666 dsic.upv.es!antodo 442 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&ritz);CHKERRQ(ierr);
902 dsic.upv.es!antodo 443 
  ierr = PetscMalloc(ncv*ncv*sizeof(PetscReal),&Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 444 
  ierr = PetscMalloc(2*ncv*sizeof(int),&isuppz);CHKERRQ(ierr);
445 
  lwork = 20*ncv;
446 
  ierr = PetscMalloc(lwork*sizeof(PetscReal),&work);CHKERRQ(ierr);
447 
  liwork = 10*ncv;
448 
  ierr = PetscMalloc(liwork*sizeof(int),&iwork);CHKERRQ(ierr);
449 
 
450 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&bnd);CHKERRQ(ierr);
687 dsic.upv.es!antodo 451 
  ierr = PetscMalloc(ncv*sizeof(int),&perm);CHKERRQ(ierr);
897 dsic.upv.es!antodo 452 
  ierr = PetscMalloc(ncv*sizeof(char),&conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 453 
 
655 dsic.upv.es!jroman 454 
  /* The first Lanczos vector is the normalized initial vector */
772 dsic.upv.es!antodo 455 
  ierr = EPSGetStartVector(eps,0,eps->V[0]);CHKERRQ(ierr);
647 dsic.upv.es!antodo 456 
 
895 dsic.upv.es!antodo 457 
  abstol = 2.0*LAPACKlamch_("S",1);
891 dsic.upv.es!antodo 458 
  anorm = 1.0;
655 dsic.upv.es!jroman 459 
  nconv = 0;
647 dsic.upv.es!antodo 460 
  eps->its = 0;
902 dsic.upv.es!antodo 461 
  restart_ritz = 0.0;
687 dsic.upv.es!antodo 462 
  for (i=0;i<eps->ncv;i++) eps->eigi[i]=0.0;
704 dsic.upv.es!antodo 463 
  EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,ncv);
464 
 
655 dsic.upv.es!jroman 465 
  /* Restart loop */
891 dsic.upv.es!antodo 466 
  eps->reason = EPS_CONVERGED_ITERATING;
467 
  while (eps->reason == EPS_CONVERGED_ITERATING) {
655 dsic.upv.es!jroman 468 
    /* Compute an ncv-step Lanczos factorization */
772 dsic.upv.es!antodo 469 
    m = ncv;
891 dsic.upv.es!antodo 470 
    ierr = EPSBasicLanczos(eps,T,eps->V,nconv,&m,f,&beta,&breakdown,anorm);CHKERRQ(ierr);
647 dsic.upv.es!antodo 471 
 
655 dsic.upv.es!jroman 472 
    /* At this point, T has the following structure
473 
 
474 
              | *     |               |
475 
              |     * |               |
476 
              | ------|-------------- |
477 
          T = |       | *   *         |
478 
              |       | *   *   *     |
479 
              |       |     *   *   * |
480 
              |       |         *   * |
481 
 
482 
       that is, a real symmetric tridiagonal matrix of order ncv whose
483 
       principal submatrix of order nconv is diagonal.  */
484 
 
485 
    /* Extract the tridiagonal block from T. Store the diagonal elements in D
486 
       and the off-diagonal elements in E  */
891 dsic.upv.es!antodo 487 
 
772 dsic.upv.es!antodo 488 
    n = m - nconv;
891 dsic.upv.es!antodo 489 
    for (i=0;i<n-1;i++) {
895 dsic.upv.es!antodo 490 
      D[i] = PetscRealPart(T[(i+nconv)*(ncv+1)]);
891 dsic.upv.es!antodo 491 
      E[i] = PetscRealPart(T[(i+nconv)*(ncv+1)+1]);
647 dsic.upv.es!antodo 492 
    }
895 dsic.upv.es!antodo 493 
    D[n-1] = PetscRealPart(T[(n-1+nconv)*(ncv+1)]);
655 dsic.upv.es!jroman 494 
 
668 dsic.upv.es!antodo 495 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
895 dsic.upv.es!antodo 496 
    LAPACKstevr_("V","A",&n,D,E,&vl,&vu,PETSC_NULL,PETSC_NULL,&abstol,&mout,ritz,Y,&n,isuppz,work,&lwork,iwork,&liwork,&info,1,1);
497 
    if (info) SETERRQ1(PETSC_ERR_LIB,"Error in Lapack xSTEVR %i",info);
891 dsic.upv.es!antodo 498 
 
499 
    /* Estimate ||A|| */
500 
    for (i=0;i<n;i++)
501 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
647 dsic.upv.es!antodo 502 
 
902 dsic.upv.es!antodo 503 
    /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
897 dsic.upv.es!antodo 504 
    for (i=0;i<n;i++)
902 dsic.upv.es!antodo 505 
      bnd[i] = beta*PetscAbsReal(Y[i*n+n-1]) + PETSC_MACHINE_EPSILON*anorm;
891 dsic.upv.es!antodo 506 
 
870 dsic.upv.es!antodo 507 
#ifdef PETSC_USE_COMPLEX
666 dsic.upv.es!antodo 508 
    for (i=0;i<n;i++) {
870 dsic.upv.es!antodo 509 
      eps->eigr[i+nconv] = ritz[i];
666 dsic.upv.es!antodo 510 
    }
870 dsic.upv.es!antodo 511 
    ierr = EPSSortEigenvalues(n,eps->eigr+nconv,eps->eigi,eps->which,n,perm);CHKERRQ(ierr);
512 
#else
513 
    ierr = EPSSortEigenvalues(n,ritz,eps->eigi,eps->which,n,perm);CHKERRQ(ierr);
514 
#endif
515 
 
649 dsic.upv.es!antodo 516 
    /* Look for converged eigenpairs */
870 dsic.upv.es!antodo 517 
    k = nconv;
518 
    for (i=0;i<n;i++) {
891 dsic.upv.es!antodo 519 
      eps->eigr[k] = ritz[perm[i]];
902 dsic.upv.es!antodo 520 
      eps->errest[k] = bnd[perm[i]] / PetscAbsScalar(eps->eigr[k]);    
891 dsic.upv.es!antodo 521 
      if (eps->errest[k] < eps->tol) {
902 dsic.upv.es!antodo 522 
 
523 
        if (lanczos->reorthog == EPSLANCZOS_ORTHOG_NONE) {
524 
          if (i>0 && PetscAbsScalar((eps->eigr[k]-ritz[perm[i-1]])/eps->eigr[k]) < eps->tol) {
525 
            /* Discard repeated eigenvalues */
526 
            conv[i] = 'R';
527 
            continue;
528 
          }
529 
        }
530 
 
870 dsic.upv.es!antodo 531 
        ierr = VecSet(eps->AV[k],0.0);CHKERRQ(ierr);
902 dsic.upv.es!antodo 532 
#ifndef PETSC_USE_COMPLEX
870 dsic.upv.es!antodo 533 
        ierr = VecMAXPY(eps->AV[k],n,Y+perm[i]*n,eps->V+nconv);CHKERRQ(ierr);
902 dsic.upv.es!antodo 534 
#else
535 
        for (j=0;j<n;j++) {
536 
          ierr = VecAXPY(eps->AV[k],Y[perm[i]*n+j],eps->V[nconv+j]);CHKERRQ(ierr);
537 
        }
538 
#endif   
897 dsic.upv.es!antodo 539 
 
540 
        if (lanczos->reorthog == EPSLANCZOS_ORTHOG_NONE) {
902 dsic.upv.es!antodo 541 
          ierr = VecNorm(eps->AV[k],NORM_2,&norm);CHKERRQ(ierr);
542 
          ierr = VecScale(eps->AV[k],1.0/norm);CHKERRQ(ierr);
543 
          eps->errest[k] = eps->errest[k] / norm;
544 
          if (i<n-1 && PetscAbsScalar((eps->eigr[k]-ritz[perm[i+1]])/eps->eigr[k]) >= eps->tol) {
545 
            PetscReal res;
546 
            ierr = STApply(eps->OP,eps->AV[k],f);CHKERRQ(ierr);
547 
            ierr = VecAXPY(f,-eps->eigr[k],eps->AV[k]);CHKERRQ(ierr);
548 
            ierr = VecNorm(f,NORM_2,&res);CHKERRQ(ierr);
549 
            eps->errest[k] = res / PetscAbsScalar(eps->eigr[k]);
550 
          }
551 
          if (eps->errest[k] >= eps->tol) {
552 
            conv[i] = 'S';
553 
            continue;
897 dsic.upv.es!antodo 554 
          }
895 dsic.upv.es!antodo 555 
        }
902 dsic.upv.es!antodo 556 
 
557 
        conv[i] = 'C';
558 
        k++;
559 
      } else conv[i] = 'N';
687 dsic.upv.es!antodo 560 
    }
647 dsic.upv.es!antodo 561 
 
891 dsic.upv.es!antodo 562 
    /* Look for unconverged eigenpairs */
563 
    j = k;
564 
    restart = -1;
565 
    for (i=0;i<n;i++) {
897 dsic.upv.es!antodo 566 
      if (conv[i] != 'C') {
567 
        if (restart == -1 && conv[i] == 'N') restart = i;
891 dsic.upv.es!antodo 568 
        eps->eigr[j] = ritz[perm[i]];
569 
        eps->errest[j] = bnd[perm[i]] / ritz[perm[i]];
570 
        j++;
571 
      }
572 
    }
902 dsic.upv.es!antodo 573 
 
574 
    if (breakdown) {
575 
      restart = -1;
576 
      PetscLogInfo((eps,"Breakdown in Lanczos method (norm=%g)\n",beta));
577 
    }
891 dsic.upv.es!antodo 578 
 
897 dsic.upv.es!antodo 579 
    if (k<eps->nev) {
580 
      if (lanczos->reorthog == EPSLANCZOS_ORTHOG_SELECTIVE && restart != -1) {
902 dsic.upv.es!antodo 581 
        /* avoid stagnation in selective reorthogonalization */
897 dsic.upv.es!antodo 582 
        if (PetscAbsScalar(restart_ritz - ritz[perm[restart]]) < PETSC_MACHINE_EPSILON) {
583 
          restart = -1;
584 
          restart_ritz = 0;
585 
        } else restart_ritz = ritz[perm[restart]];
586 
      }
587 
      if (restart != -1) {
588 
        /* use first not converged vector for restarting */
589 
        ierr = VecSet(eps->AV[k],0.0);CHKERRQ(ierr);
902 dsic.upv.es!antodo 590 
#ifndef PETSC_USE_COMPLEX
897 dsic.upv.es!antodo 591 
        ierr = VecMAXPY(eps->AV[k],n,Y+perm[restart]*n,eps->V+nconv);CHKERRQ(ierr);
902 dsic.upv.es!antodo 592 
#else
593 
        for (j=0;j<n;j++) {
594 
          ierr = VecAXPY(eps->AV[k],Y[perm[restart]*n+j],eps->V[nconv+j]);CHKERRQ(ierr);
595 
        }
596 
#endif
897 dsic.upv.es!antodo 597 
        ierr = VecCopy(eps->AV[k],eps->V[k]);CHKERRQ(ierr);
598 
      } else {
599 
        /* use random vector for restarting */
600 
        ierr = SlepcVecSetRandom(eps->V[k]);CHKERRQ(ierr);
601 
        PetscLogInfo((eps,"Using random vector for restart\n"));
602 
      }
870 dsic.upv.es!antodo 603 
    }
604 
 
605 
    /* copy converged vectors to V */
606 
    for (i=nconv;i<k;i++) {
607 
      ierr = VecCopy(eps->AV[i],eps->V[i]);CHKERRQ(ierr);
608 
    }
609 
 
897 dsic.upv.es!antodo 610 
    if (k<eps->nev) {
611 
      if (restart == -1 || lanczos->reorthog == EPSLANCZOS_ORTHOG_NONE) {
612 
        /* reorthonormalize restart vector */
895 dsic.upv.es!antodo 613 
        ierr = EPSOrthogonalize(eps,eps->nds+k,eps->DSV,eps->V[k],PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
897 dsic.upv.es!antodo 614 
        if (breakdown) {
615 
          eps->reason = EPS_DIVERGED_BREAKDOWN;
616 
          PetscLogInfo((eps,"Unable to generate more start vectors\n"));
617 
        } else {
618 
          ierr = VecScale(eps->V[k],1.0/norm);CHKERRQ(ierr);
619 
        }
891 dsic.upv.es!antodo 620 
      }
621 
    }
622 
 
870 dsic.upv.es!antodo 623 
    nconv = k;
891 dsic.upv.es!antodo 624 
    EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nconv+n);
625 
    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
626 
    if (nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
647 dsic.upv.es!antodo 627 
  }
628 
 
655 dsic.upv.es!jroman 629 
  eps->nconv = nconv;
647 dsic.upv.es!antodo 630 
 
895 dsic.upv.es!antodo 631 
  ierr = PetscFree(D);CHKERRQ(ierr);
632 
  ierr = PetscFree(E);CHKERRQ(ierr);
666 dsic.upv.es!antodo 633 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
647 dsic.upv.es!antodo 634 
  ierr = PetscFree(Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 635 
  ierr = PetscFree(isuppz);CHKERRQ(ierr);
636 
  ierr = PetscFree(work);CHKERRQ(ierr);
637 
  ierr = PetscFree(iwork);CHKERRQ(ierr);
638 
 
639 
  ierr = PetscFree(bnd);CHKERRQ(ierr);
687 dsic.upv.es!antodo 640 
  ierr = PetscFree(perm);CHKERRQ(ierr);
895 dsic.upv.es!antodo 641 
  ierr = PetscFree(conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 642 
  PetscFunctionReturn(0);
806 dsic.upv.es!antodo 643 
#endif
647 dsic.upv.es!antodo 644 
}
645 
 
671 dsic.upv.es!antodo 646 
#undef __FUNCT__  
647 
#define __FUNCT__ "EPSSetFromOptions_LANCZOS"
648 
PetscErrorCode EPSSetFromOptions_LANCZOS(EPS eps)
649 
{
650 
  PetscErrorCode ierr;
651 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
652 
  PetscTruth     flg;
906 dsic.upv.es!antodo 653 
  const char     *list[5] = { "none", "full", "selective", "periodic", "partial" };
671 dsic.upv.es!antodo 654 
 
655 
  PetscFunctionBegin;
656 
  ierr = PetscOptionsHead("LANCZOS options");CHKERRQ(ierr);
906 dsic.upv.es!antodo 657 
  ierr = PetscOptionsEList("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",list,5,list[lanczos->reorthog],(int*)&lanczos->reorthog,&flg);CHKERRQ(ierr);
671 dsic.upv.es!antodo 658 
  ierr = PetscOptionsTail();CHKERRQ(ierr);
659 
  PetscFunctionReturn(0);
660 
}
661 
 
647 dsic.upv.es!antodo 662 
EXTERN_C_BEGIN
663 
#undef __FUNCT__  
906 dsic.upv.es!antodo 664 
#define __FUNCT__ "EPSLanczosSetReorthog_LANCZOS"
665 
PetscErrorCode EPSLanczosSetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 666 
{
667 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
668 
 
669 
  PetscFunctionBegin;
670 
  switch (reorthog) {
671 
    case EPSLANCZOS_ORTHOG_NONE:
672 
    case EPSLANCZOS_ORTHOG_FULL:
772 dsic.upv.es!antodo 673 
    case EPSLANCZOS_ORTHOG_SELECTIVE:
746 dsic.upv.es!antodo 674 
    case EPSLANCZOS_ORTHOG_PERIODIC:
675 
    case EPSLANCZOS_ORTHOG_PARTIAL:
671 dsic.upv.es!antodo 676 
      lanczos->reorthog = reorthog;
677 
      break;
678 
    default:
679 
      SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
680 
  }
681 
  PetscFunctionReturn(0);
682 
}
683 
EXTERN_C_END
684 
 
685 
#undef __FUNCT__  
906 dsic.upv.es!antodo 686 
#define __FUNCT__ "EPSLanczosSetReorthog"
671 dsic.upv.es!antodo 687 
/*@
906 dsic.upv.es!antodo 688 
   EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 689 
   iteration.
690 
 
691 
   Collective on EPS
692 
 
693 
   Input Parameters:
694 
+  eps - the eigenproblem solver context
695 
-  reorthog - the type of reorthogonalization
696 
 
697 
   Options Database Key:
698 
.  -eps_lanczos_orthog - Sets the reorthogonalization type (either 'none' or
699 
                           'full')
700 
 
701 
   Level: advanced
702 
 
906 dsic.upv.es!antodo 703 
.seealso: EPSLanczosGetReorthog()
671 dsic.upv.es!antodo 704 
@*/
906 dsic.upv.es!antodo 705 
PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 706 
{
906 dsic.upv.es!antodo 707 
  PetscErrorCode ierr, (*f)(EPS,EPSLanczosReorthogType);
671 dsic.upv.es!antodo 708 
 
709 
  PetscFunctionBegin;
710 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
906 dsic.upv.es!antodo 711 
  ierr = PetscObjectQueryFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",(void (**)())&f);CHKERRQ(ierr);
671 dsic.upv.es!antodo 712 
  if (f) {
713 
    ierr = (*f)(eps,reorthog);CHKERRQ(ierr);
714 
  }
715 
  PetscFunctionReturn(0);
716 
}
717 
 
718 
EXTERN_C_BEGIN
719 
#undef __FUNCT__  
906 dsic.upv.es!antodo 720 
#define __FUNCT__ "EPSLanczosGetReorthog_LANCZOS"
721 
PetscErrorCode EPSLanczosGetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 722 
{
723 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
724 
  PetscFunctionBegin;
725 
  *reorthog = lanczos->reorthog;
726 
  PetscFunctionReturn(0);
727 
}
728 
EXTERN_C_END
729 
 
730 
#undef __FUNCT__  
906 dsic.upv.es!antodo 731 
#define __FUNCT__ "EPSLanczosGetReorthog"
707 dsic.upv.es!antodo 732 
/*@C
906 dsic.upv.es!antodo 733 
   EPSLanczosGetReorthog - Gets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 734 
   iteration.
735 
 
736 
   Collective on EPS
737 
 
738 
   Input Parameter:
739 
.  eps - the eigenproblem solver context
740 
 
741 
   Input Parameter:
742 
.  reorthog - the type of reorthogonalization
743 
 
744 
   Level: advanced
745 
 
906 dsic.upv.es!antodo 746 
.seealso: EPSLanczosSetReorthog()
671 dsic.upv.es!antodo 747 
@*/
906 dsic.upv.es!antodo 748 
PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 749 
{
906 dsic.upv.es!antodo 750 
  PetscErrorCode ierr, (*f)(EPS,EPSLanczosReorthogType*);
671 dsic.upv.es!antodo 751 
 
752 
  PetscFunctionBegin;
753 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
906 dsic.upv.es!antodo 754 
  ierr = PetscObjectQueryFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",(void (**)())&f);CHKERRQ(ierr);
671 dsic.upv.es!antodo 755 
  if (f) {
756 
    ierr = (*f)(eps,reorthog);CHKERRQ(ierr);
757 
  }
758 
  PetscFunctionReturn(0);
759 
}
760 
 
761 
#undef __FUNCT__  
762 
#define __FUNCT__ "EPSView_LANCZOS"
763 
PetscErrorCode EPSView_LANCZOS(EPS eps,PetscViewer viewer)
764 
{
765 
  PetscErrorCode ierr;
766 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
767 
  PetscTruth     isascii;
906 dsic.upv.es!antodo 768 
  const char     *list[5] = { "none", "full" , "selective", "periodic", "partial" };
671 dsic.upv.es!antodo 769 
 
770 
  PetscFunctionBegin;
771 
  ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
772 
  if (!isascii) {
773 
    SETERRQ1(1,"Viewer type %s not supported for EPSLANCZOS",((PetscObject)viewer)->type_name);
774 
  }  
775 
  ierr = PetscViewerASCIIPrintf(viewer,"reorthogonalization: %s\n",list[lanczos->reorthog]);CHKERRQ(ierr);
776 
  PetscFunctionReturn(0);
777 
}
778 
 
906 dsic.upv.es!antodo 779 
EXTERN PetscErrorCode EPSSolve_TS_LANCZOS(EPS);
780 
 
671 dsic.upv.es!antodo 781 
EXTERN_C_BEGIN
782 
#undef __FUNCT__  
647 dsic.upv.es!antodo 783 
#define __FUNCT__ "EPSCreate_LANCZOS"
784 
PetscErrorCode EPSCreate_LANCZOS(EPS eps)
785 
{
671 dsic.upv.es!antodo 786 
  PetscErrorCode ierr;
787 
  EPS_LANCZOS    *lanczos;
788 
 
647 dsic.upv.es!antodo 789 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 790 
  ierr = PetscNew(EPS_LANCZOS,&lanczos);CHKERRQ(ierr);
791 
  PetscMemzero(lanczos,sizeof(EPS_LANCZOS));
792 
  PetscLogObjectMemory(eps,sizeof(EPS_LANCZOS));
793 
  eps->data                      = (void *) lanczos;
647 dsic.upv.es!antodo 794 
  eps->ops->solve                = EPSSolve_LANCZOS;
885 dsic.upv.es!jroman 795 
  eps->ops->solvets              = EPSSolve_TS_LANCZOS;
647 dsic.upv.es!antodo 796 
  eps->ops->setup                = EPSSetUp_LANCZOS;
671 dsic.upv.es!antodo 797 
  eps->ops->setfromoptions       = EPSSetFromOptions_LANCZOS;
647 dsic.upv.es!antodo 798 
  eps->ops->destroy              = EPSDestroy_Default;
671 dsic.upv.es!antodo 799 
  eps->ops->view                 = EPSView_LANCZOS;
647 dsic.upv.es!antodo 800 
  eps->ops->backtransform        = EPSBackTransform_Default;
885 dsic.upv.es!jroman 801 
  if (eps->solverclass==EPS_TWO_SIDE)
802 
       eps->ops->computevectors       = EPSComputeVectors_Schur;
803 
  else eps->ops->computevectors       = EPSComputeVectors_Default;
671 dsic.upv.es!antodo 804 
  lanczos->reorthog              = EPSLANCZOS_ORTHOG_FULL;
906 dsic.upv.es!antodo 805 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","EPSLanczosSetReorthog_LANCZOS",EPSLanczosSetReorthog_LANCZOS);CHKERRQ(ierr);
806 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","EPSLanczosGetReorthog_LANCZOS",EPSLanczosGetReorthog_LANCZOS);CHKERRQ(ierr);
647 dsic.upv.es!antodo 807 
  PetscFunctionReturn(0);
808 
}
809 
EXTERN_C_END
810