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 
 
950 dsic.upv.es!jroman 7 
   Algorithm:
647 dsic.upv.es!antodo 8 
 
950 dsic.upv.es!jroman 9 
       Lanczos method for symmetric (Hermitian) problems, with explicit
10 
       restart and deflation. Several reorthogonalization strategies can
655 dsic.upv.es!jroman 11 
       be selected.
12 
 
647 dsic.upv.es!antodo 13 
   References:
14 
 
950 dsic.upv.es!jroman 15 
       [1] "Lanczos Methods in SLEPc", SLEPc Technical Report STR-5,
16 
           available at http://www.grycap.upv.es/slepc.
647 dsic.upv.es!antodo 17 
 
1673 slepc 18 
   Last update: Feb 2009
655 dsic.upv.es!jroman 19 
 
1376 slepc 20 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1672 slepc 21 
   SLEPc - Scalable Library for Eigenvalue Problem Computations
2116 eromero 22 
   Copyright (c) 2002-2010, Universidad Politecnica de Valencia, Spain
1376 slepc 23 
 
1672 slepc 24 
   This file is part of SLEPc.
25 
 
26 
   SLEPc is free software: you can redistribute it and/or modify it under  the
27 
   terms of version 3 of the GNU Lesser General Public License as published by
28 
   the Free Software Foundation.
29 
 
30 
   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
31 
   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
32 
   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
33 
   more details.
34 
 
35 
   You  should have received a copy of the GNU Lesser General  Public  License
36 
   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
1376 slepc 37 
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
647 dsic.upv.es!antodo 38 
*/
1376 slepc 39 
 
2283 jroman 40 
#include <private/epsimpl.h>                /*I "slepceps.h" I*/
647 dsic.upv.es!antodo 41 
 
1947 jroman 42 
PetscErrorCode EPSSolve_LANCZOS(EPS);
43 
 
671 dsic.upv.es!antodo 44 
typedef struct {
906 dsic.upv.es!antodo 45 
  EPSLanczosReorthogType reorthog;
1953 jroman 46 
  Vec *AV;
671 dsic.upv.es!antodo 47 
} EPS_LANCZOS;
48 
 
647 dsic.upv.es!antodo 49 
#undef __FUNCT__  
50 
#define __FUNCT__ "EPSSetUp_LANCZOS"
51 
PetscErrorCode EPSSetUp_LANCZOS(EPS eps)
52 
{
1638 slepc 53 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
647 dsic.upv.es!antodo 54 
  PetscErrorCode ierr;
55 
 
56 
  PetscFunctionBegin;
1577 slepc 57 
  if (eps->ncv) { /* ncv set */
2214 jroman 58 
    if (eps->ncv<eps->nev) SETERRQ(((PetscObject)eps)->comm,1,"The value of ncv must be at least nev");
647 dsic.upv.es!antodo 59 
  }
1577 slepc 60 
  else if (eps->mpd) { /* mpd set */
1928 jroman 61 
    eps->ncv = PetscMin(eps->n,eps->nev+eps->mpd);
1577 slepc 62 
  }
63 
  else { /* neither set: defaults depend on nev being small or large */
1928 jroman 64 
    if (eps->nev<500) eps->ncv = PetscMin(eps->n,PetscMax(2*eps->nev,eps->nev+15));
65 
    else { eps->mpd = 500; eps->ncv = PetscMin(eps->n,eps->nev+eps->mpd); }
1577 slepc 66 
  }
67 
  if (!eps->mpd) eps->mpd = eps->ncv;
2214 jroman 68 
  if (eps->ncv>eps->nev+eps->mpd) SETERRQ(((PetscObject)eps)->comm,1,"The value of ncv must not be larger than nev+mpd");
1928 jroman 69 
  if (!eps->max_it) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
1220 slepc 70 
 
1942 jroman 71 
  if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
1947 jroman 72 
  switch (eps->which) {
73 
    case EPS_LARGEST_IMAGINARY:
74 
    case EPS_SMALLEST_IMAGINARY:
75 
    case EPS_TARGET_IMAGINARY:
2214 jroman 76 
      SETERRQ(((PetscObject)eps)->comm,1,"Wrong value of eps->which");
1947 jroman 77 
    default: ; /* default case to remove warning */
906 dsic.upv.es!antodo 78 
  }
1947 jroman 79 
  if (!eps->ishermitian)
2214 jroman 80 
    SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
1560 slepc 81 
  if (!eps->extraction) {
82 
    ierr = EPSSetExtraction(eps,EPS_RITZ);CHKERRQ(ierr);
83 
  } else if (eps->extraction!=EPS_RITZ) {
2214 jroman 84 
    SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Unsupported extraction type\n");
1426 slepc 85 
  }
86 
 
647 dsic.upv.es!antodo 87 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
1940 jroman 88 
  if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_SELECTIVE) {
1953 jroman 89 
    ierr = VecDuplicateVecs(eps->V[0],eps->ncv,&lanczos->AV);CHKERRQ(ierr);
1638 slepc 90 
  }
1947 jroman 91 
  if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
1755 antodo 92 
    ierr = EPSDefaultGetWork(eps,2);CHKERRQ(ierr);
93 
  } else {
94 
    ierr = EPSDefaultGetWork(eps,1);CHKERRQ(ierr);
95 
  }
1947 jroman 96 
 
97 
  /* dispatch solve method */
2214 jroman 98 
  if (eps->leftvecs) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Left vectors not supported in this solver");
1947 jroman 99 
  eps->ops->solve = EPSSolve_LANCZOS;
647 dsic.upv.es!antodo 100 
  PetscFunctionReturn(0);
101 
}
102 
 
103 
#undef __FUNCT__  
1068 slepc 104 
#define __FUNCT__ "EPSLocalLanczos"
683 dsic.upv.es!jroman 105 
/*
1068 slepc 106 
   EPSLocalLanczos - Local reorthogonalization.
683 dsic.upv.es!jroman 107 
 
108 
   This is the simplest variant. At each Lanczos step, the corresponding Lanczos vector
109 
   is orthogonalized with respect to the two previous Lanczos vectors, according to
110 
   the three term Lanczos recurrence. WARNING: This variant does not track the loss of
111 
   orthogonality that occurs in finite-precision arithmetic and, therefore, the
112 
   generated vectors are not guaranteed to be (semi-)orthogonal.
113 
*/
2216 jroman 114 
static PetscErrorCode EPSLocalLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,Vec *V,PetscInt k,PetscInt *M,Vec f,PetscBool *breakdown)
671 dsic.upv.es!antodo 115 
{
116 
  PetscErrorCode ierr;
1509 slepc 117 
  PetscInt       i,j,m = *M;
1133 slepc 118 
  PetscReal      norm;
2216 jroman 119 
  PetscBool      *which,lwhich[100];
1755 antodo 120 
  PetscScalar    *hwork,lhwork[100];
1121 slepc 121 
 
1163 slepc 122 
  PetscFunctionBegin;  
1755 antodo 123 
  if (m > 100) {
2216 jroman 124 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which);CHKERRQ(ierr);
1755 antodo 125 
    ierr = PetscMalloc(m*sizeof(PetscScalar),&hwork);CHKERRQ(ierr);
1616 slepc 126 
  } else {
127 
    which = lwhich;
128 
    hwork = lhwork;
129 
  }
1755 antodo 130 
  for (i=0;i<k;i++)
1163 slepc 131 
    which[i] = PETSC_TRUE;
1121 slepc 132 
 
1616 slepc 133 
  for (j=k;j<m-1;j++) {
134 
    ierr = STApply(eps->OP,V[j],V[j+1]);CHKERRQ(ierr);
1755 antodo 135 
    which[j] = PETSC_TRUE;
136 
    if (j-2>=k) which[j-2] = PETSC_FALSE;
137 
    ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,which,V,V[j+1],hwork,&norm,breakdown);CHKERRQ(ierr);
138 
    alpha[j-k] = PetscRealPart(hwork[j]);
1616 slepc 139 
    beta[j-k] = norm;
1121 slepc 140 
    if (*breakdown) {
141 
      *M = j+1;
1755 antodo 142 
      if (m > 100) {
1616 slepc 143 
        ierr = PetscFree(which);CHKERRQ(ierr);
144 
        ierr = PetscFree(hwork);CHKERRQ(ierr);
145 
      }
146 
      PetscFunctionReturn(0);
147 
    } else {
148 
      ierr = VecScale(V[j+1],1.0/norm);CHKERRQ(ierr);
1121 slepc 149 
    }
150 
  }
1616 slepc 151 
  ierr = STApply(eps->OP,V[m-1],f);CHKERRQ(ierr);
1755 antodo 152 
  ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,m,PETSC_NULL,V,f,hwork,&norm,PETSC_NULL);CHKERRQ(ierr);
153 
  alpha[m-1-k] = PetscRealPart(hwork[m-1]);
1639 slepc 154 
  beta[m-1-k] = norm;
1133 slepc 155 
 
1755 antodo 156 
  if (m > 100) {
1616 slepc 157 
    ierr = PetscFree(which);CHKERRQ(ierr);
158 
    ierr = PetscFree(hwork);CHKERRQ(ierr);
159 
  }
1121 slepc 160 
  PetscFunctionReturn(0);
161 
}
162 
 
163 
#undef __FUNCT__  
891 dsic.upv.es!antodo 164 
#define __FUNCT__ "EPSSelectiveLanczos"
928 dsic.upv.es!antodo 165 
/*
166 
   EPSSelectiveLanczos - Selective reorthogonalization.
167 
*/
2216 jroman 168 
static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,Vec *V,PetscInt k,PetscInt *M,Vec f,PetscBool *breakdown,PetscReal anorm)
891 dsic.upv.es!antodo 169 
{
170 
  PetscErrorCode ierr;
1953 jroman 171 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
1509 slepc 172 
  PetscInt       i,j,m = *M,n,nritz=0,nritzo;
1616 slepc 173 
  PetscReal      *d,*e,*ritz,norm;
1755 antodo 174 
  PetscScalar    *Y,*hwork,lhwork[100];
2216 jroman 175 
  PetscBool      *which,lwhich[100];
891 dsic.upv.es!antodo 176 
 
177 
  PetscFunctionBegin;
1616 slepc 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);
1178 slepc 181 
  ierr = PetscMalloc(m*m*sizeof(PetscScalar),&Y);CHKERRQ(ierr);
1755 antodo 182 
  if (m > 100) {
2216 jroman 183 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which);CHKERRQ(ierr);
1755 antodo 184 
    ierr = PetscMalloc(m*sizeof(PetscScalar),&hwork);CHKERRQ(ierr);
1616 slepc 185 
  } else {
186 
    which = lwhich;
187 
    hwork = lhwork;
188 
  }
1755 antodo 189 
  for (i=0;i<k;i++)
1163 slepc 190 
    which[i] = PETSC_TRUE;
1124 slepc 191 
 
891 dsic.upv.es!antodo 192 
  for (j=k;j<m;j++) {
1154 slepc 193 
    /* Lanczos step */
891 dsic.upv.es!antodo 194 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1755 antodo 195 
    which[j] = PETSC_TRUE;
196 
    if (j-2>=k) which[j-2] = PETSC_FALSE;
197 
    ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,which,V,f,hwork,&norm,breakdown);CHKERRQ(ierr);
198 
    alpha[j-k] = PetscRealPart(hwork[j]);
1616 slepc 199 
    beta[j-k] = norm;
891 dsic.upv.es!antodo 200 
    if (*breakdown) {
201 
      *M = j+1;
202 
      break;
203 
    }
204 
 
1154 slepc 205 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
891 dsic.upv.es!antodo 206 
    n = j-k+1;
1616 slepc 207 
    for (i=0;i<n;i++) { d[i] = alpha[i]; e[i] = beta[i]; }
208 
    ierr = EPSDenseTridiagonal(n,d,e,ritz,Y);CHKERRQ(ierr);
891 dsic.upv.es!antodo 209 
 
210 
    /* Estimate ||A|| */
211 
    for (i=0;i<n;i++)
212 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
1154 slepc 213 
 
214 
    /* Compute nearly converged Ritz vectors */
215 
    nritzo = 0;
216 
    for (i=0;i<n;i++)
891 dsic.upv.es!antodo 217 
      if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm)
1154 slepc 218 
        nritzo++;
219 
 
220 
    if (nritzo>nritz) {
221 
      nritz = 0;
222 
      for (i=0;i<n;i++) {
223 
        if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
1953 jroman 224 
          ierr = SlepcVecMAXPBY(lanczos->AV[nritz],0.0,1.0,n,Y+i*n,V+k);CHKERRQ(ierr);
1154 slepc 225 
          nritz++;
226 
        }
227 
      }
891 dsic.upv.es!antodo 228 
    }
229 
 
1154 slepc 230 
    if (nritz > 0) {
1953 jroman 231 
      ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,nritz,PETSC_NULL,lanczos->AV,f,hwork,&norm,breakdown);CHKERRQ(ierr);
1154 slepc 232 
      if (*breakdown) {
233 
        *M = j+1;
234 
        break;
235 
      }
236 
    }
237 
 
891 dsic.upv.es!antodo 238 
    if (j<m-1) {
239 
      ierr = VecScale(f,1.0 / norm);CHKERRQ(ierr);
240 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
241 
    }
671 dsic.upv.es!antodo 242 
  }
1133 slepc 243 
 
1616 slepc 244 
  ierr = PetscFree(d);CHKERRQ(ierr);
245 
  ierr = PetscFree(e);CHKERRQ(ierr);
891 dsic.upv.es!antodo 246 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 247 
  ierr = PetscFree(Y);CHKERRQ(ierr);
1755 antodo 248 
  if (m > 100) {
1616 slepc 249 
    ierr = PetscFree(which);CHKERRQ(ierr);
250 
    ierr = PetscFree(hwork);CHKERRQ(ierr);
251 
  }
671 dsic.upv.es!antodo 252 
  PetscFunctionReturn(0);
253 
}
254 
 
255 
#undef __FUNCT__  
1126 slepc 256 
#define __FUNCT__ "update_omega"
1509 slepc 257 
static void update_omega(PetscReal *omega,PetscReal *omega_old,PetscInt j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
772 dsic.upv.es!antodo 258 
{
1509 slepc 259 
  PetscInt       k;
2139 jroman 260 
  PetscReal      T,binv;
1121 slepc 261 
 
262 
  PetscFunctionBegin;
1126 slepc 263 
  /* Estimate of contribution to roundoff errors from A*v
264 
       fl(A*v) = A*v + f,
265 
     where ||f|| \approx eps1*||A||.
266 
     For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps. */
267 
  T = eps1*anorm;
268 
  binv = 1.0/beta[j+1];
269 
 
270 
  /* Update omega(1) using omega(0)==0. */
1178 slepc 271 
  omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] -
1126 slepc 272 
                beta[j]*omega_old[0];
273 
  if (omega_old[0] > 0)
274 
    omega_old[0] = binv*(omega_old[0] + T);
275 
  else
276 
    omega_old[0] = binv*(omega_old[0] - T);  
277 
 
278 
  /* Update remaining components. */
279 
  for (k=1;k<j-1;k++) {
1178 slepc 280 
    omega_old[k] = beta[k+1]*omega[k+1] + (alpha[k]-alpha[j])*omega[k] +
1126 slepc 281 
                   beta[k]*omega[k-1] - beta[j]*omega_old[k];
282 
    if (omega_old[k] > 0)
283 
      omega_old[k] = binv*(omega_old[k] + T);      
284 
    else
285 
      omega_old[k] = binv*(omega_old[k] - T);      
1121 slepc 286 
  }
1126 slepc 287 
  omega_old[j-1] = binv*T;
1121 slepc 288 
 
1126 slepc 289 
  /* Swap omega and omega_old. */
290 
  for (k=0;k<j;k++) {
291 
    omega[k] = omega_old[k];
292 
    omega_old[k] = omega[k];
1121 slepc 293 
  }
1126 slepc 294 
  omega[j] = eps1;  
1128 slepc 295 
  PetscFunctionReturnVoid();
1121 slepc 296 
}
297 
 
298 
#undef __FUNCT__  
1128 slepc 299 
#define __FUNCT__ "compute_int"
2216 jroman 300 
static void compute_int(PetscBool *which,PetscReal *mu,PetscInt j,PetscReal delta,PetscReal eta)
1128 slepc 301 
{
2216 jroman 302 
  PetscInt  i,k,maxpos;
303 
  PetscReal max;
304 
  PetscBool found;
1128 slepc 305 
 
1129 slepc 306 
  PetscFunctionBegin;  
1128 slepc 307 
  /* initialize which */
308 
  found = PETSC_FALSE;
1297 slepc 309 
  maxpos = 0;
1128 slepc 310 
  max = 0.0;
311 
  for (i=0;i<j;i++) {
1178 slepc 312 
    if (PetscAbsReal(mu[i]) >= delta) {
1128 slepc 313 
      which[i] = PETSC_TRUE;
314 
      found = PETSC_TRUE;
315 
    } else which[i] = PETSC_FALSE;
1178 slepc 316 
    if (PetscAbsReal(mu[i]) > max) {
1128 slepc 317 
      maxpos = i;
1178 slepc 318 
      max = PetscAbsReal(mu[i]);
1128 slepc 319 
    }
320 
  }
321 
  if (!found) which[maxpos] = PETSC_TRUE;    
322 
 
323 
  for (i=0;i<j;i++)
324 
    if (which[i]) {
325 
      /* find left interval */
326 
      for (k=i;k>=0;k--) {
1178 slepc 327 
        if (PetscAbsReal(mu[k])<eta || which[k]) break;
1128 slepc 328 
        else which[k] = PETSC_TRUE;
329 
      }
330 
      /* find right interval */
331 
      for (k=i+1;k<j;k++) {
1178 slepc 332 
        if (PetscAbsReal(mu[k])<eta || which[k]) break;
1128 slepc 333 
        else which[k] = PETSC_TRUE;
334 
      }
335 
    }
336 
  PetscFunctionReturnVoid();
337 
}
338 
 
339 
#undef __FUNCT__  
891 dsic.upv.es!antodo 340 
#define __FUNCT__ "EPSPartialLanczos"
928 dsic.upv.es!antodo 341 
/*
342 
   EPSPartialLanczos - Partial reorthogonalization.
343 
*/
2216 jroman 344 
static PetscErrorCode EPSPartialLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,Vec *V,PetscInt k,PetscInt *M,Vec f,PetscBool *breakdown,PetscReal anorm)
746 dsic.upv.es!antodo 345 
{
1929 jroman 346 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
746 dsic.upv.es!antodo 347 
  PetscErrorCode ierr;
1509 slepc 348 
  PetscInt       i,j,m = *M;
1629 slepc 349 
  PetscReal      norm,*omega,lomega[100],*omega_old,lomega_old[100],eps1,delta,eta;
2216 jroman 350 
  PetscBool      *which,lwhich[100],*which2,lwhich2[100],
1163 slepc 351 
                 reorth = PETSC_FALSE,force_reorth = PETSC_FALSE,fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
1755 antodo 352 
  PetscScalar    *hwork,lhwork[100];
746 dsic.upv.es!antodo 353 
 
354 
  PetscFunctionBegin;
1163 slepc 355 
  if (m>100) {
1178 slepc 356 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega);CHKERRQ(ierr);
357 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega_old);CHKERRQ(ierr);
1128 slepc 358 
  } else {
359 
    omega = lomega;
360 
    omega_old = lomega_old;
1629 slepc 361 
  }
1755 antodo 362 
  if (m > 100) {
2216 jroman 363 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which);CHKERRQ(ierr);
364 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which2);CHKERRQ(ierr);
1755 antodo 365 
    ierr = PetscMalloc(m*sizeof(PetscScalar),&hwork);CHKERRQ(ierr);
1629 slepc 366 
  } else {
1128 slepc 367 
    which = lwhich;
1163 slepc 368 
    which2 = lwhich2;
1629 slepc 369 
    hwork = lhwork;
1128 slepc 370 
  }
1629 slepc 371 
 
1929 jroman 372 
  eps1 = sqrt((PetscReal)eps->n)*PETSC_MACHINE_EPSILON/2;
1126 slepc 373 
  delta = PETSC_SQRT_MACHINE_EPSILON/sqrt((PetscReal)eps->ncv);
1128 slepc 374 
  eta = pow(PETSC_MACHINE_EPSILON,3.0/4.0)/sqrt((PetscReal)eps->ncv);
1143 slepc 375 
  if (anorm < 0.0) {
376 
    anorm = 1.0;
377 
    estimate_anorm = PETSC_TRUE;
378 
  }
1163 slepc 379 
  for (i=0;i<m-k;i++)
1152 slepc 380 
    omega[i] = omega_old[i] = 0.0;
1755 antodo 381 
  for (i=0;i<k;i++)
1163 slepc 382 
    which[i] = PETSC_TRUE;  
1126 slepc 383 
 
891 dsic.upv.es!antodo 384 
  for (j=k;j<m;j++) {
385 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1128 slepc 386 
    if (fro) {
387 
      /* Lanczos step with full reorthogonalization */
1755 antodo 388 
      ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,PETSC_NULL,V,f,hwork,&norm,breakdown);CHKERRQ(ierr);      
389 
      alpha[j-k] = PetscRealPart(hwork[j]);
1128 slepc 390 
    } else {
391 
      /* Lanczos step */
1755 antodo 392 
      which[j] = PETSC_TRUE;
393 
      if (j-2>=k) which[j-2] = PETSC_FALSE;
394 
      ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,which,V,f,hwork,&norm,breakdown);CHKERRQ(ierr);
395 
      alpha[j-k] = PetscRealPart(hwork[j]);
1629 slepc 396 
      beta[j-k] = norm;
1149 slepc 397 
 
1154 slepc 398 
      /* Estimate ||A|| if needed */
1152 slepc 399 
      if (estimate_anorm) {
1629 slepc 400 
        if (j>k) anorm = PetscMax(anorm,PetscAbsReal(alpha[j-k])+norm+beta[j-k-1]);
401 
        else anorm = PetscMax(anorm,PetscAbsReal(alpha[j-k])+norm);
1152 slepc 402 
      }
1154 slepc 403 
 
404 
      /* Check if reorthogonalization is needed */
1128 slepc 405 
      reorth = PETSC_FALSE;
1154 slepc 406 
      if (j>k) {      
1629 slepc 407 
        update_omega(omega,omega_old,j-k,alpha,beta-1,eps1,anorm);
1154 slepc 408 
        for (i=0;i<j-k;i++)
409 
          if (PetscAbsScalar(omega[i]) > delta) reorth = PETSC_TRUE;
410 
      }
1128 slepc 411 
 
412 
      if (reorth || force_reorth) {
1940 jroman 413 
        if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
1128 slepc 414 
          /* Periodic reorthogonalization */
1129 slepc 415 
          if (force_reorth) force_reorth = PETSC_FALSE;
416 
          else force_reorth = PETSC_TRUE;
1755 antodo 417 
          ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,j-k,PETSC_NULL,V+k,f,hwork,&norm,breakdown);CHKERRQ(ierr);
1128 slepc 418 
          for (i=0;i<j-k;i++)
419 
            omega[i] = eps1;
420 
        } else {
421 
          /* Partial reorthogonalization */
1129 slepc 422 
          if (force_reorth) force_reorth = PETSC_FALSE;
423 
          else {
424 
            force_reorth = PETSC_TRUE;
1163 slepc 425 
            compute_int(which2,omega,j-k,delta,eta);
1129 slepc 426 
            for (i=0;i<j-k;i++)
1163 slepc 427 
              if (which2[i]) omega[i] = eps1;
1150 slepc 428 
          }
1755 antodo 429 
          ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,j-k,which2,V+k,f,hwork,&norm,breakdown);CHKERRQ(ierr);   
1129 slepc 430 
        }      
1126 slepc 431 
      }
746 dsic.upv.es!antodo 432 
    }
1150 slepc 433 
 
1929 jroman 434 
    if (*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON) {
891 dsic.upv.es!antodo 435 
      *M = j+1;
436 
      break;
437 
    }
1128 slepc 438 
    if (!fro && norm*delta < anorm*eps1) {
439 
      fro = PETSC_TRUE;
440 
      PetscInfo1(eps,"Switching to full reorthogonalization at iteration %i\n",eps->its);      
441 
    }
1639 slepc 442 
 
443 
    beta[j-k] = norm;
891 dsic.upv.es!antodo 444 
    if (j<m-1) {
445 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
446 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
746 dsic.upv.es!antodo 447 
    }
448 
  }
1126 slepc 449 
 
1163 slepc 450 
  if (m>100) {
1128 slepc 451 
    ierr = PetscFree(omega);CHKERRQ(ierr);
452 
    ierr = PetscFree(omega_old);CHKERRQ(ierr);
1629 slepc 453 
  }
1755 antodo 454 
  if (m > 100) {
1128 slepc 455 
    ierr = PetscFree(which);CHKERRQ(ierr);
1163 slepc 456 
    ierr = PetscFree(which2);CHKERRQ(ierr);
1629 slepc 457 
    ierr = PetscFree(hwork);CHKERRQ(ierr);
1128 slepc 458 
  }
746 dsic.upv.es!antodo 459 
  PetscFunctionReturn(0);
460 
}
461 
 
462 
#undef __FUNCT__  
655 dsic.upv.es!jroman 463 
#define __FUNCT__ "EPSBasicLanczos"
464 
/*
465 
   EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
466 
   columns are assumed to be locked and therefore they are not modified. On
467 
   exit, the following relation is satisfied:
647 dsic.upv.es!antodo 468 
 
655 dsic.upv.es!jroman 469 
                    OP * V - V * T = f * e_m^T
470 
 
471 
   where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
472 
   f is the residual vector and e_m is the m-th vector of the canonical basis.
473 
   The Lanczos vectors (together with vector f) are B-orthogonal (to working
474 
   accuracy) if full reorthogonalization is being used, otherwise they are
475 
   (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
476 
   Lanczos vector can be computed as v_{m+1} = f / beta.
477 
 
478 
   This function simply calls another function which depends on the selected
479 
   reorthogonalization strategy.
480 
*/
2216 jroman 481 
static PetscErrorCode EPSBasicLanczos(EPS eps,PetscReal *alpha,PetscReal *beta,Vec *V,PetscInt k,PetscInt *m,Vec f,PetscBool *breakdown,PetscReal anorm)
655 dsic.upv.es!jroman 482 
{
671 dsic.upv.es!antodo 483 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
655 dsic.upv.es!jroman 484 
  PetscErrorCode ierr;
1345 slepc 485 
  IPOrthogonalizationRefinementType orthog_ref;
1616 slepc 486 
  PetscScalar *T;
487 
  PetscInt i,n=*m;
1662 slepc 488 
  PetscReal betam;
655 dsic.upv.es!jroman 489 
 
490 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 491 
  switch (lanczos->reorthog) {
1940 jroman 492 
    case EPS_LANCZOS_REORTHOG_LOCAL:
1616 slepc 493 
      ierr = EPSLocalLanczos(eps,alpha,beta,V,k,m,f,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 494 
      break;
1940 jroman 495 
    case EPS_LANCZOS_REORTHOG_SELECTIVE:
1616 slepc 496 
      ierr = EPSSelectiveLanczos(eps,alpha,beta,V,k,m,f,breakdown,anorm);CHKERRQ(ierr);
671 dsic.upv.es!antodo 497 
      break;
1940 jroman 498 
    case EPS_LANCZOS_REORTHOG_FULL:
1616 slepc 499 
      ierr = EPSFullLanczos(eps,alpha,beta,V,k,m,f,breakdown);CHKERRQ(ierr);
1124 slepc 500 
      break;
1940 jroman 501 
    case EPS_LANCZOS_REORTHOG_PARTIAL:
502 
    case EPS_LANCZOS_REORTHOG_PERIODIC:
1629 slepc 503 
      ierr = EPSPartialLanczos(eps,alpha,beta,V,k,m,f,breakdown,anorm);CHKERRQ(ierr);
504 
      break;
1940 jroman 505 
    case EPS_LANCZOS_REORTHOG_DELAYED:
1616 slepc 506 
      ierr = PetscMalloc(n*n*sizeof(PetscScalar),&T);CHKERRQ(ierr);
1345 slepc 507 
      ierr = IPGetOrthogonalization(eps->ip,PETSC_NULL,&orthog_ref,PETSC_NULL);CHKERRQ(ierr);
1629 slepc 508 
      if (orthog_ref == IP_ORTH_REFINE_NEVER) {
1616 slepc 509 
        ierr = EPSDelayedArnoldi1(eps,T,n,V,k,m,f,&betam,breakdown);CHKERRQ(ierr);      
1124 slepc 510 
      } else {
1616 slepc 511 
        ierr = EPSDelayedArnoldi(eps,T,n,V,k,m,f,&betam,breakdown);CHKERRQ(ierr);
1121 slepc 512 
      }
1682 slepc 513 
      for (i=k;i<n-1;i++) { alpha[i-k] = PetscRealPart(T[n*i+i]); beta[i-k] = PetscRealPart(T[n*i+i+1]); }
514 
      alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
1616 slepc 515 
      beta[n-1] = betam;
516 
      ierr = PetscFree(T);CHKERRQ(ierr);
772 dsic.upv.es!antodo 517 
      break;
1124 slepc 518 
    default:
2214 jroman 519 
      SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 520 
  }
655 dsic.upv.es!jroman 521 
  PetscFunctionReturn(0);
522 
}
523 
 
647 dsic.upv.es!antodo 524 
#undef __FUNCT__  
525 
#define __FUNCT__ "EPSSolve_LANCZOS"
526 
PetscErrorCode EPSSolve_LANCZOS(EPS eps)
527 
{
891 dsic.upv.es!antodo 528 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
647 dsic.upv.es!antodo 529 
  PetscErrorCode ierr;
2161 jroman 530 
  PetscInt       nconv,i,j,k,l,x,n,m,*perm,restart,ncv=eps->ncv,r;
1755 antodo 531 
  Vec            w=eps->work[1],f=eps->work[0];
1623 slepc 532 
  PetscScalar    *Y,stmp;
2070 jroman 533 
  PetscReal      *d,*e,*ritz,*bnd,anorm,beta,norm,rtmp,resnorm;
2216 jroman 534 
  PetscBool      breakdown;
1638 slepc 535 
  char           *conv,ctmp;
647 dsic.upv.es!antodo 536 
 
537 
  PetscFunctionBegin;
1616 slepc 538 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&d);CHKERRQ(ierr);
539 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&e);CHKERRQ(ierr);
666 dsic.upv.es!antodo 540 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&ritz);CHKERRQ(ierr);
1178 slepc 541 
  ierr = PetscMalloc(ncv*ncv*sizeof(PetscScalar),&Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 542 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&bnd);CHKERRQ(ierr);
1638 slepc 543 
  ierr = PetscMalloc(ncv*sizeof(PetscInt),&perm);CHKERRQ(ierr);
897 dsic.upv.es!antodo 544 
  ierr = PetscMalloc(ncv*sizeof(char),&conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 545 
 
655 dsic.upv.es!jroman 546 
  /* The first Lanczos vector is the normalized initial vector */
1057 slepc 547 
  ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);
647 dsic.upv.es!antodo 548 
 
1143 slepc 549 
  anorm = -1.0;
655 dsic.upv.es!jroman 550 
  nconv = 0;
704 dsic.upv.es!antodo 551 
 
655 dsic.upv.es!jroman 552 
  /* Restart loop */
891 dsic.upv.es!antodo 553 
  while (eps->reason == EPS_CONVERGED_ITERATING) {
1220 slepc 554 
    eps->its++;
655 dsic.upv.es!jroman 555 
    /* Compute an ncv-step Lanczos factorization */
1577 slepc 556 
    m = PetscMin(nconv+eps->mpd,ncv);
1616 slepc 557 
    ierr = EPSBasicLanczos(eps,d,e,eps->V,nconv,&m,f,&breakdown,anorm);CHKERRQ(ierr);
647 dsic.upv.es!antodo 558 
 
1178 slepc 559 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
772 dsic.upv.es!antodo 560 
    n = m - nconv;
1616 slepc 561 
    beta = e[n-1];
562 
    ierr = EPSDenseTridiagonal(n,d,e,ritz,Y);CHKERRQ(ierr);
1638 slepc 563 
 
891 dsic.upv.es!antodo 564 
    /* Estimate ||A|| */
565 
    for (i=0;i<n;i++)
566 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
647 dsic.upv.es!antodo 567 
 
902 dsic.upv.es!antodo 568 
    /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
1638 slepc 569 
    for (i=0;i<n;i++) {
2070 jroman 570 
      resnorm = beta*PetscAbsScalar(Y[i*n+n-1]) + PETSC_MACHINE_EPSILON*anorm;
571 
      ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
572 
      if (bnd[i]<eps->tol) {
902 dsic.upv.es!antodo 573 
        conv[i] = 'C';
1638 slepc 574 
      } else {
575 
        conv[i] = 'N';
576 
      }
687 dsic.upv.es!antodo 577 
    }
647 dsic.upv.es!antodo 578 
 
1638 slepc 579 
    /* purge repeated ritz values */
1940 jroman 580 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL)
1638 slepc 581 
      for (i=1;i<n;i++)
582 
        if (conv[i] == 'C')
583 
          if (PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol)
584 
            conv[i] = 'R';
902 dsic.upv.es!antodo 585 
 
1638 slepc 586 
    /* Compute restart vector */
902 dsic.upv.es!antodo 587 
    if (breakdown) {
1121 slepc 588 
      PetscInfo2(eps,"Breakdown in Lanczos method (it=%i norm=%g)\n",eps->its,beta);
1638 slepc 589 
    } else {
590 
      restart = 0;
591 
      while (restart<n && conv[restart] != 'N') restart++;
592 
      if (restart >= n) {
593 
        breakdown = PETSC_TRUE;
594 
      } else {
595 
        for (i=restart+1;i<n;i++)
596 
          if (conv[i] == 'N') {
2161 jroman 597 
            ierr = EPSCompareEigenvalues(eps,ritz[restart],0.0,ritz[i],0.0,&r);CHKERRQ(ierr);
598 
            if (r>0) restart = i;
1638 slepc 599 
          }
1755 antodo 600 
        ierr = SlepcVecMAXPBY(f,0.0,1.0,n,Y+restart*n,eps->V+nconv);CHKERRQ(ierr);
1638 slepc 601 
      }
902 dsic.upv.es!antodo 602 
    }
1638 slepc 603 
 
604 
    /* Count and put converged eigenvalues first */
605 
    for (i=0;i<n;i++) perm[i] = i;
606 
    for (k=0;k<n;k++)
607 
      if (conv[perm[k]] != 'C') {
608 
        j = k + 1;
609 
        while (j<n && conv[perm[j]] != 'C') j++;
610 
        if (j>=n) break;
611 
        l = perm[k]; perm[k] = perm[j]; perm[j] = l;
612 
      }
613 
 
614 
    /* Sort eigenvectors according to permutation */
615 
    for (i=0;i<k;i++) {
616 
      x = perm[i];
617 
      if (x != i) {
618 
        j = i + 1;
619 
        while (perm[j] != i) j++;
620 
        /* swap eigenvalues i and j */
621 
        rtmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = rtmp;
622 
        rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
623 
        ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
624 
        perm[j] = x; perm[i] = i;
625 
        /* swap eigenvectors i and j */
626 
        for (l=0;l<n;l++) {
627 
          stmp = Y[x*n+l]; Y[x*n+l] = Y[i*n+l]; Y[i*n+l] = stmp;
628 
        }
629 
      }
630 
    }
891 dsic.upv.es!antodo 631 
 
1638 slepc 632 
    /* compute converged eigenvectors */
633 
    ierr = SlepcUpdateVectors(n,eps->V+nconv,0,k,Y,n,PETSC_FALSE);CHKERRQ(ierr);
634 
 
635 
    /* purge spurious ritz values */
1940 jroman 636 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
1638 slepc 637 
      for (i=0;i<k;i++) {
638 
        ierr = VecNorm(eps->V[nconv+i],NORM_2,&norm);CHKERRQ(ierr);
639 
        ierr = VecScale(eps->V[nconv+i],1.0/norm);CHKERRQ(ierr);
640 
        ierr = STApply(eps->OP,eps->V[nconv+i],w);CHKERRQ(ierr);
641 
        ierr = VecAXPY(w,-ritz[i],eps->V[nconv+i]);CHKERRQ(ierr);
642 
        ierr = VecNorm(w,NORM_2,&norm);CHKERRQ(ierr);
2070 jroman 643 
        ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
644 
        if (bnd[i]>=eps->tol) conv[i] = 'S';
897 dsic.upv.es!antodo 645 
      }
1638 slepc 646 
      for (i=0;i<k;i++)
647 
        if (conv[i] != 'C') {
648 
          j = i + 1;
649 
          while (j<k && conv[j] != 'C') j++;
650 
          if (j>=k) break;
651 
          /* swap eigenvalues i and j */
652 
          rtmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = rtmp;
653 
          rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
654 
          ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
655 
          /* swap eigenvectors i and j */
656 
          ierr = VecSwap(eps->V[nconv+i],eps->V[nconv+j]);CHKERRQ(ierr);
657 
        }
658 
      k = i;
870 dsic.upv.es!antodo 659 
    }
660 
 
1638 slepc 661 
    /* store ritz values and estimated errors */
662 
    for (i=0;i<n;i++) {
663 
      eps->eigr[nconv+i] = ritz[i];
664 
      eps->errest[nconv+i] = bnd[i];
870 dsic.upv.es!antodo 665 
    }
1638 slepc 666 
    EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nconv+n);
667 
    nconv = nconv+k;
668 
    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
669 
    if (nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
670 
 
671 
     if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
1940 jroman 672 
      if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
1638 slepc 673 
        /* Reorthonormalize restart vector */
1755 antodo 674 
        ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,nconv,PETSC_NULL,eps->V,f,PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
1638 slepc 675 
        ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
676 
      }
677 
      if (breakdown) {
1057 slepc 678 
        /* Use random vector for restarting */
679 
        PetscInfo(eps,"Using random vector for restart\n");
1638 slepc 680 
        ierr = EPSGetStartVector(eps,nconv,f,&breakdown);CHKERRQ(ierr);
1057 slepc 681 
      }
1638 slepc 682 
      if (breakdown) { /* give up */
683 
        eps->reason = EPS_DIVERGED_BREAKDOWN;
684 
        PetscInfo(eps,"Unable to generate more start vectors\n");
685 
      } else {
686 
        ierr = VecCopy(f,eps->V[nconv]);CHKERRQ(ierr);
687 
      }
688 
    }    
647 dsic.upv.es!antodo 689 
  }
690 
 
655 dsic.upv.es!jroman 691 
  eps->nconv = nconv;
647 dsic.upv.es!antodo 692 
 
1638 slepc 693 
  ierr = PetscFree(d);CHKERRQ(ierr);
694 
  ierr = PetscFree(e);CHKERRQ(ierr);
666 dsic.upv.es!antodo 695 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
647 dsic.upv.es!antodo 696 
  ierr = PetscFree(Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 697 
  ierr = PetscFree(bnd);CHKERRQ(ierr);
1638 slepc 698 
  ierr = PetscFree(perm);CHKERRQ(ierr);
895 dsic.upv.es!antodo 699 
  ierr = PetscFree(conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 700 
  PetscFunctionReturn(0);
701 
}
702 
 
1173 slepc 703 
static const char *lanczoslist[6] = { "local", "full", "selective", "periodic", "partial" , "delayed" };
1068 slepc 704 
 
671 dsic.upv.es!antodo 705 
#undef __FUNCT__  
706 
#define __FUNCT__ "EPSSetFromOptions_LANCZOS"
707 
PetscErrorCode EPSSetFromOptions_LANCZOS(EPS eps)
708 
{
709 
  PetscErrorCode ierr;
710 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
2216 jroman 711 
  PetscBool      flg;
982 slepc 712 
  PetscInt       i;
671 dsic.upv.es!antodo 713 
 
714 
  PetscFunctionBegin;
2102 eromero 715 
  ierr = PetscOptionsBegin(((PetscObject)eps)->comm,((PetscObject)eps)->prefix,"LANCZOS Options","EPS");CHKERRQ(ierr);
1173 slepc 716 
  ierr = PetscOptionsEList("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",lanczoslist,6,lanczoslist[lanczos->reorthog],&i,&flg);CHKERRQ(ierr);
1002 slepc 717 
  if (flg) lanczos->reorthog = (EPSLanczosReorthogType)i;
2102 eromero 718 
  ierr = PetscOptionsEnd();CHKERRQ(ierr);
671 dsic.upv.es!antodo 719 
  PetscFunctionReturn(0);
720 
}
721 
 
647 dsic.upv.es!antodo 722 
EXTERN_C_BEGIN
723 
#undef __FUNCT__  
906 dsic.upv.es!antodo 724 
#define __FUNCT__ "EPSLanczosSetReorthog_LANCZOS"
725 
PetscErrorCode EPSLanczosSetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 726 
{
727 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
728 
 
729 
  PetscFunctionBegin;
730 
  switch (reorthog) {
1940 jroman 731 
    case EPS_LANCZOS_REORTHOG_LOCAL:
732 
    case EPS_LANCZOS_REORTHOG_FULL:
733 
    case EPS_LANCZOS_REORTHOG_DELAYED:
734 
    case EPS_LANCZOS_REORTHOG_SELECTIVE:
735 
    case EPS_LANCZOS_REORTHOG_PERIODIC:
736 
    case EPS_LANCZOS_REORTHOG_PARTIAL:
671 dsic.upv.es!antodo 737 
      lanczos->reorthog = reorthog;
738 
      break;
739 
    default:
2214 jroman 740 
      SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 741 
  }
742 
  PetscFunctionReturn(0);
743 
}
744 
EXTERN_C_END
745 
 
746 
#undef __FUNCT__  
906 dsic.upv.es!antodo 747 
#define __FUNCT__ "EPSLanczosSetReorthog"
671 dsic.upv.es!antodo 748 
/*@
906 dsic.upv.es!antodo 749 
   EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 750 
   iteration.
751 
 
752 
   Collective on EPS
753 
 
754 
   Input Parameters:
755 
+  eps - the eigenproblem solver context
756 
-  reorthog - the type of reorthogonalization
757 
 
758 
   Options Database Key:
1486 slepc 759 
.  -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
760 
                         'periodic', 'partial', 'full' or 'delayed')
671 dsic.upv.es!antodo 761 
 
762 
   Level: advanced
763 
 
1364 slepc 764 
.seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 765 
@*/
906 dsic.upv.es!antodo 766 
PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 767 
{
2221 jroman 768 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 769 
 
770 
  PetscFunctionBegin;
2213 jroman 771 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2221 jroman 772 
  ierr = PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 773 
  PetscFunctionReturn(0);
774 
}
775 
 
776 
EXTERN_C_BEGIN
777 
#undef __FUNCT__  
906 dsic.upv.es!antodo 778 
#define __FUNCT__ "EPSLanczosGetReorthog_LANCZOS"
779 
PetscErrorCode EPSLanczosGetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 780 
{
781 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
782 
  PetscFunctionBegin;
783 
  *reorthog = lanczos->reorthog;
784 
  PetscFunctionReturn(0);
785 
}
786 
EXTERN_C_END
787 
 
788 
#undef __FUNCT__  
906 dsic.upv.es!antodo 789 
#define __FUNCT__ "EPSLanczosGetReorthog"
707 dsic.upv.es!antodo 790 
/*@C
906 dsic.upv.es!antodo 791 
   EPSLanczosGetReorthog - Gets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 792 
   iteration.
793 
 
794 
   Collective on EPS
795 
 
796 
   Input Parameter:
797 
.  eps - the eigenproblem solver context
798 
 
799 
   Input Parameter:
800 
.  reorthog - the type of reorthogonalization
801 
 
802 
   Level: advanced
803 
 
1364 slepc 804 
.seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 805 
@*/
906 dsic.upv.es!antodo 806 
PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 807 
{
2221 jroman 808 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 809 
 
810 
  PetscFunctionBegin;
2213 jroman 811 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2221 jroman 812 
  ierr = PetscTryMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 813 
  PetscFunctionReturn(0);
814 
}
815 
 
816 
#undef __FUNCT__  
1925 jroman 817 
#define __FUNCT__ "EPSDestroy_LANCZOS"
818 
PetscErrorCode EPSDestroy_LANCZOS(EPS eps)
819 
{
820 
  PetscErrorCode ierr;
1953 jroman 821 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
1925 jroman 822 
 
823 
  PetscFunctionBegin;
2213 jroman 824 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2271 jroman 825 
  if (lanczos->AV) { ierr = VecDestroyVecs(eps->ncv,&lanczos->AV);CHKERRQ(ierr); }
1925 jroman 826 
  ierr = EPSDestroy_Default(eps);CHKERRQ(ierr);
827 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
828 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
829 
  PetscFunctionReturn(0);
830 
}
831 
 
832 
#undef __FUNCT__  
671 dsic.upv.es!antodo 833 
#define __FUNCT__ "EPSView_LANCZOS"
834 
PetscErrorCode EPSView_LANCZOS(EPS eps,PetscViewer viewer)
835 
{
836 
  PetscErrorCode ierr;
837 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
2216 jroman 838 
  PetscBool      isascii;
671 dsic.upv.es!antodo 839 
 
840 
  PetscFunctionBegin;
2215 jroman 841 
  ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
671 dsic.upv.es!antodo 842 
  if (!isascii) {
2214 jroman 843 
    SETERRQ1(((PetscObject)eps)->comm,1,"Viewer type %s not supported for EPSLANCZOS",((PetscObject)viewer)->type_name);
671 dsic.upv.es!antodo 844 
  }  
1068 slepc 845 
  ierr = PetscViewerASCIIPrintf(viewer,"reorthogonalization: %s\n",lanczoslist[lanczos->reorthog]);CHKERRQ(ierr);
671 dsic.upv.es!antodo 846 
  PetscFunctionReturn(0);
847 
}
848 
 
849 
EXTERN_C_BEGIN
850 
#undef __FUNCT__  
647 dsic.upv.es!antodo 851 
#define __FUNCT__ "EPSCreate_LANCZOS"
852 
PetscErrorCode EPSCreate_LANCZOS(EPS eps)
853 
{
671 dsic.upv.es!antodo 854 
  PetscErrorCode ierr;
855 
  EPS_LANCZOS    *lanczos;
856 
 
647 dsic.upv.es!antodo 857 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 858 
  ierr = PetscNew(EPS_LANCZOS,&lanczos);CHKERRQ(ierr);
859 
  PetscLogObjectMemory(eps,sizeof(EPS_LANCZOS));
860 
  eps->data                      = (void *) lanczos;
647 dsic.upv.es!antodo 861 
  eps->ops->setup                = EPSSetUp_LANCZOS;
671 dsic.upv.es!antodo 862 
  eps->ops->setfromoptions       = EPSSetFromOptions_LANCZOS;
1925 jroman 863 
  eps->ops->destroy              = EPSDestroy_LANCZOS;
671 dsic.upv.es!antodo 864 
  eps->ops->view                 = EPSView_LANCZOS;
647 dsic.upv.es!antodo 865 
  eps->ops->backtransform        = EPSBackTransform_Default;
1928 jroman 866 
  eps->ops->computevectors       = EPSComputeVectors_Hermitian;
1940 jroman 867 
  lanczos->reorthog              = EPS_LANCZOS_REORTHOG_LOCAL;
906 dsic.upv.es!antodo 868 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","EPSLanczosSetReorthog_LANCZOS",EPSLanczosSetReorthog_LANCZOS);CHKERRQ(ierr);
869 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","EPSLanczosGetReorthog_LANCZOS",EPSLanczosGetReorthog_LANCZOS);CHKERRQ(ierr);
647 dsic.upv.es!antodo 870 
  PetscFunctionReturn(0);
871 
}
872 
EXTERN_C_END
873