Subversion Repositories slepc-dev

Rev

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
2575 eromero 22 
   Copyright (c) 2002-2011, Universitat 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 
 
2729 jroman 40 
#include <slepc-private/epsimpl.h>                /*I "slepceps.h" I*/
647 dsic.upv.es!antodo 41 
 
2324 jroman 42 
PetscErrorCode EPSSolve_Lanczos(EPS);
1947 jroman 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__  
2324 jroman 50 
#define __FUNCT__ "EPSSetUp_Lanczos"
51 
PetscErrorCode EPSSetUp_Lanczos(EPS eps)
647 dsic.upv.es!antodo 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 
 
2613 jroman 71 
  if (!eps->which) { ierr = EPSDefaultSetWhich(eps);CHKERRQ(ierr); }
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 
  }
2762 jroman 79 
  if (!eps->ishermitian) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
1560 slepc 80 
  if (!eps->extraction) {
81 
    ierr = EPSSetExtraction(eps,EPS_RITZ);CHKERRQ(ierr);
2762 jroman 82 
  } else if (eps->extraction!=EPS_RITZ) SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_SUP,"Unsupported extraction type\n");
1426 slepc 83 
 
647 dsic.upv.es!antodo 84 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
1940 jroman 85 
  if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_SELECTIVE) {
2410 jroman 86 
    ierr = VecDuplicateVecs(eps->t,eps->ncv,&lanczos->AV);CHKERRQ(ierr);
1638 slepc 87 
  }
2796 jroman 88 
  ierr = PSSetType(eps->ps,PSHEP);CHKERRQ(ierr);
89 
  ierr = PSSetCompact(eps->ps,PETSC_TRUE);CHKERRQ(ierr);
2792 jroman 90 
  ierr = PSAllocate(eps->ps,eps->ncv);CHKERRQ(ierr);
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");
2324 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)
2316 jroman 218 
        nritzo++;
1154 slepc 219 
 
220 
    if (nritzo>nritz) {
221 
      nritz = 0;
222 
      for (i=0;i<n;i++) {
2316 jroman 223 
        if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
224 
          ierr = SlepcVecMAXPBY(lanczos->AV[nritz],0.0,1.0,n,Y+i*n,V+k);CHKERRQ(ierr);
1154 slepc 225 
          nritz++;
2316 jroman 226 
        }
1154 slepc 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) {
2316 jroman 233 
        *M = j+1;
234 
        break;
1154 slepc 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;
2316 jroman 328 
        else which[k] = PETSC_TRUE;
1128 slepc 329 
      }
330 
      /* find right interval */
331 
      for (k=i+1;k<j;k++) {
1178 slepc 332 
        if (PetscAbsReal(mu[k])<eta || which[k]) break;
2316 jroman 333 
        else which[k] = PETSC_TRUE;
1128 slepc 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],
2317 jroman 351 
                 reorth = PETSC_FALSE,force_reorth = PETSC_FALSE,
352 
                 fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
1755 antodo 353 
  PetscScalar    *hwork,lhwork[100];
746 dsic.upv.es!antodo 354 
 
355 
  PetscFunctionBegin;
1163 slepc 356 
  if (m>100) {
1178 slepc 357 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega);CHKERRQ(ierr);
358 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega_old);CHKERRQ(ierr);
1128 slepc 359 
  } else {
360 
    omega = lomega;
361 
    omega_old = lomega_old;
1629 slepc 362 
  }
1755 antodo 363 
  if (m > 100) {
2216 jroman 364 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which);CHKERRQ(ierr);
365 
    ierr = PetscMalloc(sizeof(PetscBool)*m,&which2);CHKERRQ(ierr);
1755 antodo 366 
    ierr = PetscMalloc(m*sizeof(PetscScalar),&hwork);CHKERRQ(ierr);
1629 slepc 367 
  } else {
1128 slepc 368 
    which = lwhich;
1163 slepc 369 
    which2 = lwhich2;
1629 slepc 370 
    hwork = lhwork;
1128 slepc 371 
  }
1629 slepc 372 
 
2473 jroman 373 
  eps1 = PetscSqrtReal((PetscReal)eps->n)*PETSC_MACHINE_EPSILON/2;
374 
  delta = PETSC_SQRT_MACHINE_EPSILON/PetscSqrtReal((PetscReal)eps->ncv);
375 
  eta = pow(PETSC_MACHINE_EPSILON,3.0/4.0)/PetscSqrtReal((PetscReal)eps->ncv);
1143 slepc 376 
  if (anorm < 0.0) {
377 
    anorm = 1.0;
378 
    estimate_anorm = PETSC_TRUE;
379 
  }
1163 slepc 380 
  for (i=0;i<m-k;i++)
1152 slepc 381 
    omega[i] = omega_old[i] = 0.0;
1755 antodo 382 
  for (i=0;i<k;i++)
1163 slepc 383 
    which[i] = PETSC_TRUE;  
1126 slepc 384 
 
891 dsic.upv.es!antodo 385 
  for (j=k;j<m;j++) {
386 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1128 slepc 387 
    if (fro) {
388 
      /* Lanczos step with full reorthogonalization */
1755 antodo 389 
      ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,PETSC_NULL,V,f,hwork,&norm,breakdown);CHKERRQ(ierr);      
390 
      alpha[j-k] = PetscRealPart(hwork[j]);
1128 slepc 391 
    } else {
392 
      /* Lanczos step */
1755 antodo 393 
      which[j] = PETSC_TRUE;
394 
      if (j-2>=k) which[j-2] = PETSC_FALSE;
395 
      ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,j+1,which,V,f,hwork,&norm,breakdown);CHKERRQ(ierr);
396 
      alpha[j-k] = PetscRealPart(hwork[j]);
1629 slepc 397 
      beta[j-k] = norm;
1149 slepc 398 
 
1154 slepc 399 
      /* Estimate ||A|| if needed */
1152 slepc 400 
      if (estimate_anorm) {
1629 slepc 401 
        if (j>k) anorm = PetscMax(anorm,PetscAbsReal(alpha[j-k])+norm+beta[j-k-1]);
2316 jroman 402 
        else anorm = PetscMax(anorm,PetscAbsReal(alpha[j-k])+norm);
1152 slepc 403 
      }
1154 slepc 404 
 
405 
      /* Check if reorthogonalization is needed */
1128 slepc 406 
      reorth = PETSC_FALSE;
1154 slepc 407 
      if (j>k) {      
2316 jroman 408 
        update_omega(omega,omega_old,j-k,alpha,beta-1,eps1,anorm);
409 
        for (i=0;i<j-k;i++)
410 
          if (PetscAbsScalar(omega[i]) > delta) reorth = PETSC_TRUE;
1154 slepc 411 
      }
1128 slepc 412 
 
413 
      if (reorth || force_reorth) {
2316 jroman 414 
        if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_PERIODIC) {
415 
          /* Periodic reorthogonalization */
416 
          if (force_reorth) force_reorth = PETSC_FALSE;
417 
          else force_reorth = PETSC_TRUE;
418 
          ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,j-k,PETSC_NULL,V+k,f,hwork,&norm,breakdown);CHKERRQ(ierr);
419 
          for (i=0;i<j-k;i++)
1128 slepc 420 
            omega[i] = eps1;
2316 jroman 421 
        } else {
422 
          /* Partial reorthogonalization */
423 
          if (force_reorth) force_reorth = PETSC_FALSE;
424 
          else {
425 
            force_reorth = PETSC_TRUE;
426 
            compute_int(which2,omega,j-k,delta,eta);
427 
            for (i=0;i<j-k;i++)
428 
              if (which2[i]) omega[i] = eps1;
429 
          }
430 
          ierr = IPOrthogonalize(eps->ip,0,PETSC_NULL,j-k,which2,V+k,f,hwork,&norm,breakdown);CHKERRQ(ierr);
431 
        }
1126 slepc 432 
      }
746 dsic.upv.es!antodo 433 
    }
1150 slepc 434 
 
1929 jroman 435 
    if (*breakdown || norm < eps->n*anorm*PETSC_MACHINE_EPSILON) {
891 dsic.upv.es!antodo 436 
      *M = j+1;
437 
      break;
438 
    }
1128 slepc 439 
    if (!fro && norm*delta < anorm*eps1) {
440 
      fro = PETSC_TRUE;
2499 jroman 441 
      ierr = PetscInfo1(eps,"Switching to full reorthogonalization at iteration %D\n",eps->its);CHKERRQ(ierr);
1128 slepc 442 
    }
1639 slepc 443 
 
444 
    beta[j-k] = norm;
891 dsic.upv.es!antodo 445 
    if (j<m-1) {
446 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
447 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
746 dsic.upv.es!antodo 448 
    }
449 
  }
1126 slepc 450 
 
1163 slepc 451 
  if (m>100) {
1128 slepc 452 
    ierr = PetscFree(omega);CHKERRQ(ierr);
453 
    ierr = PetscFree(omega_old);CHKERRQ(ierr);
454 
    ierr = PetscFree(which);CHKERRQ(ierr);
1163 slepc 455 
    ierr = PetscFree(which2);CHKERRQ(ierr);
1629 slepc 456 
    ierr = PetscFree(hwork);CHKERRQ(ierr);
1128 slepc 457 
  }
746 dsic.upv.es!antodo 458 
  PetscFunctionReturn(0);
459 
}
460 
 
461 
#undef __FUNCT__  
655 dsic.upv.es!jroman 462 
#define __FUNCT__ "EPSBasicLanczos"
463 
/*
464 
   EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
465 
   columns are assumed to be locked and therefore they are not modified. On
466 
   exit, the following relation is satisfied:
647 dsic.upv.es!antodo 467 
 
655 dsic.upv.es!jroman 468 
                    OP * V - V * T = f * e_m^T
469 
 
470 
   where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
471 
   f is the residual vector and e_m is the m-th vector of the canonical basis.
472 
   The Lanczos vectors (together with vector f) are B-orthogonal (to working
473 
   accuracy) if full reorthogonalization is being used, otherwise they are
474 
   (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
475 
   Lanczos vector can be computed as v_{m+1} = f / beta.
476 
 
477 
   This function simply calls another function which depends on the selected
478 
   reorthogonalization strategy.
479 
*/
2216 jroman 480 
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 481 
{
2370 jroman 482 
  EPS_LANCZOS        *lanczos = (EPS_LANCZOS *)eps->data;
483 
  PetscScalar        *T;
484 
  PetscInt           i,n=*m;
485 
  PetscReal          betam;
486 
  PetscErrorCode     ierr;
487 
  IPOrthogRefineType orthog_ref;
655 dsic.upv.es!jroman 488 
 
489 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 490 
  switch (lanczos->reorthog) {
1940 jroman 491 
    case EPS_LANCZOS_REORTHOG_LOCAL:
1616 slepc 492 
      ierr = EPSLocalLanczos(eps,alpha,beta,V,k,m,f,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 493 
      break;
1940 jroman 494 
    case EPS_LANCZOS_REORTHOG_SELECTIVE:
1616 slepc 495 
      ierr = EPSSelectiveLanczos(eps,alpha,beta,V,k,m,f,breakdown,anorm);CHKERRQ(ierr);
671 dsic.upv.es!antodo 496 
      break;
1940 jroman 497 
    case EPS_LANCZOS_REORTHOG_FULL:
1616 slepc 498 
      ierr = EPSFullLanczos(eps,alpha,beta,V,k,m,f,breakdown);CHKERRQ(ierr);
1124 slepc 499 
      break;
1940 jroman 500 
    case EPS_LANCZOS_REORTHOG_PARTIAL:
501 
    case EPS_LANCZOS_REORTHOG_PERIODIC:
1629 slepc 502 
      ierr = EPSPartialLanczos(eps,alpha,beta,V,k,m,f,breakdown,anorm);CHKERRQ(ierr);
503 
      break;
1940 jroman 504 
    case EPS_LANCZOS_REORTHOG_DELAYED:
1616 slepc 505 
      ierr = PetscMalloc(n*n*sizeof(PetscScalar),&T);CHKERRQ(ierr);
1345 slepc 506 
      ierr = IPGetOrthogonalization(eps->ip,PETSC_NULL,&orthog_ref,PETSC_NULL);CHKERRQ(ierr);
2370 jroman 507 
      if (orthog_ref == IP_ORTHOG_REFINE_NEVER) {
1616 slepc 508 
        ierr = EPSDelayedArnoldi1(eps,T,n,V,k,m,f,&betam,breakdown);CHKERRQ(ierr);      
1124 slepc 509 
      } else {
1616 slepc 510 
        ierr = EPSDelayedArnoldi(eps,T,n,V,k,m,f,&betam,breakdown);CHKERRQ(ierr);
1121 slepc 511 
      }
1682 slepc 512 
      for (i=k;i<n-1;i++) { alpha[i-k] = PetscRealPart(T[n*i+i]); beta[i-k] = PetscRealPart(T[n*i+i+1]); }
513 
      alpha[n-1] = PetscRealPart(T[n*(n-1)+n-1]);
1616 slepc 514 
      beta[n-1] = betam;
515 
      ierr = PetscFree(T);CHKERRQ(ierr);
772 dsic.upv.es!antodo 516 
      break;
1124 slepc 517 
    default:
2214 jroman 518 
      SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 519 
  }
655 dsic.upv.es!jroman 520 
  PetscFunctionReturn(0);
521 
}
522 
 
647 dsic.upv.es!antodo 523 
#undef __FUNCT__  
2324 jroman 524 
#define __FUNCT__ "EPSSolve_Lanczos"
525 
PetscErrorCode EPSSolve_Lanczos(EPS eps)
647 dsic.upv.es!antodo 526 
{
2317 jroman 527 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
647 dsic.upv.es!antodo 528 
  PetscErrorCode ierr;
2818 jroman 529 
  PetscInt       nconv,i,j,k,l,x,n,*perm,restart,ncv=eps->ncv,r,ld;
1755 antodo 530 
  Vec            w=eps->work[1],f=eps->work[0];
2798 jroman 531 
  PetscScalar    *Y,*ritz,stmp;
532 
  PetscReal      *d,*e,*bnd,anorm,beta,norm,rtmp,resnorm;
2216 jroman 533 
  PetscBool      breakdown;
1638 slepc 534 
  char           *conv,ctmp;
647 dsic.upv.es!antodo 535 
 
536 
  PetscFunctionBegin;
2792 jroman 537 
  ierr = PSGetLeadingDimension(eps->ps,&ld);CHKERRQ(ierr);
2798 jroman 538 
  ierr = PetscMalloc(ncv*sizeof(PetscScalar),&ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 539 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&bnd);CHKERRQ(ierr);
1638 slepc 540 
  ierr = PetscMalloc(ncv*sizeof(PetscInt),&perm);CHKERRQ(ierr);
897 dsic.upv.es!antodo 541 
  ierr = PetscMalloc(ncv*sizeof(char),&conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 542 
 
655 dsic.upv.es!jroman 543 
  /* The first Lanczos vector is the normalized initial vector */
1057 slepc 544 
  ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);
647 dsic.upv.es!antodo 545 
 
1143 slepc 546 
  anorm = -1.0;
655 dsic.upv.es!jroman 547 
  nconv = 0;
704 dsic.upv.es!antodo 548 
 
655 dsic.upv.es!jroman 549 
  /* Restart loop */
891 dsic.upv.es!antodo 550 
  while (eps->reason == EPS_CONVERGED_ITERATING) {
1220 slepc 551 
    eps->its++;
2792 jroman 552 
 
655 dsic.upv.es!jroman 553 
    /* Compute an ncv-step Lanczos factorization */
2818 jroman 554 
    n = PetscMin(nconv+eps->mpd,ncv);
2792 jroman 555 
    ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&d);CHKERRQ(ierr);
556 
    e = d + ld;
2818 jroman 557 
    ierr = EPSBasicLanczos(eps,d,e,eps->V,nconv,&n,f,&breakdown,anorm);CHKERRQ(ierr);
1616 slepc 558 
    beta = e[n-1];
2792 jroman 559 
    ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&d);CHKERRQ(ierr);
560 
    ierr = PSSetDimensions(eps->ps,n,0,0);CHKERRQ(ierr);
561 
    ierr = PSSetState(eps->ps,PS_STATE_INTERMEDIATE);CHKERRQ(ierr);
562 
 
563 
    /* Solve projected problem */
564 
    ierr = PSSolve(eps->ps,ritz,PETSC_NULL);CHKERRQ(ierr);
565 
    ierr = PSSort(eps->ps,ritz,PETSC_NULL,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
1638 slepc 566 
 
891 dsic.upv.es!antodo 567 
    /* Estimate ||A|| */
2818 jroman 568 
    for (i=nconv;i<n;i++)
2798 jroman 569 
      anorm = PetscMax(anorm,PetscAbsReal(PetscRealPart(ritz[i])));
647 dsic.upv.es!antodo 570 
 
902 dsic.upv.es!antodo 571 
    /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
2792 jroman 572 
    ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
2818 jroman 573 
    for (i=nconv;i<n;i++) {
2792 jroman 574 
      resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
2070 jroman 575 
      ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
576 
      if (bnd[i]<eps->tol) {
902 dsic.upv.es!antodo 577 
        conv[i] = 'C';
1638 slepc 578 
      } else {
579 
        conv[i] = 'N';
580 
      }
687 dsic.upv.es!antodo 581 
    }
2792 jroman 582 
    ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
647 dsic.upv.es!antodo 583 
 
1638 slepc 584 
    /* purge repeated ritz values */
1940 jroman 585 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL)
2818 jroman 586 
      for (i=nconv+1;i<n;i++)
1638 slepc 587 
        if (conv[i] == 'C')
588 
          if (PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol)
589 
            conv[i] = 'R';
902 dsic.upv.es!antodo 590 
 
1638 slepc 591 
    /* Compute restart vector */
902 dsic.upv.es!antodo 592 
    if (breakdown) {
2499 jroman 593 
      ierr = PetscInfo2(eps,"Breakdown in Lanczos method (it=%D norm=%G)\n",eps->its,beta);CHKERRQ(ierr);
1638 slepc 594 
    } else {
2818 jroman 595 
      restart = nconv;
1638 slepc 596 
      while (restart<n && conv[restart] != 'N') restart++;
597 
      if (restart >= n) {
598 
        breakdown = PETSC_TRUE;
599 
      } else {
600 
        for (i=restart+1;i<n;i++)
601 
          if (conv[i] == 'N') {
2798 jroman 602 
            ierr = (*eps->which_func)(ritz[restart],0.0,ritz[i],0.0,&r,eps->which_ctx);CHKERRQ(ierr);
2161 jroman 603 
            if (r>0) restart = i;
1638 slepc 604 
          }
2792 jroman 605 
        ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
606 
        ierr = SlepcVecMAXPBY(f,0.0,1.0,n,Y+restart*ld,eps->V+nconv);CHKERRQ(ierr);
607 
        ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 608 
      }
902 dsic.upv.es!antodo 609 
    }
1638 slepc 610 
 
611 
    /* Count and put converged eigenvalues first */
2818 jroman 612 
    for (i=nconv;i<n;i++) perm[i] = i;
613 
    for (k=nconv;k<n;k++)
1638 slepc 614 
      if (conv[perm[k]] != 'C') {
615 
        j = k + 1;
616 
        while (j<n && conv[perm[j]] != 'C') j++;
617 
        if (j>=n) break;
618 
        l = perm[k]; perm[k] = perm[j]; perm[j] = l;
619 
      }
620 
 
621 
    /* Sort eigenvectors according to permutation */
2792 jroman 622 
    ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
2818 jroman 623 
    for (i=nconv;i<k;i++) {
1638 slepc 624 
      x = perm[i];
625 
      if (x != i) {
626 
        j = i + 1;
627 
        while (perm[j] != i) j++;
628 
        /* swap eigenvalues i and j */
2798 jroman 629 
        stmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = stmp;
1638 slepc 630 
        rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
631 
        ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
632 
        perm[j] = x; perm[i] = i;
633 
        /* swap eigenvectors i and j */
634 
        for (l=0;l<n;l++) {
2792 jroman 635 
          stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
1638 slepc 636 
        }
637 
      }
638 
    }
891 dsic.upv.es!antodo 639 
 
1638 slepc 640 
    /* compute converged eigenvectors */
2818 jroman 641 
    ierr = SlepcUpdateVectors(n,eps->V,nconv,k,Y,n,PETSC_FALSE);CHKERRQ(ierr);
2792 jroman 642 
    ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 643 
 
644 
    /* purge spurious ritz values */
1940 jroman 645 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
2818 jroman 646 
      for (i=nconv;i<k;i++) {
647 
        ierr = VecNorm(eps->V[i],NORM_2,&norm);CHKERRQ(ierr);
648 
        ierr = VecScale(eps->V[i],1.0/norm);CHKERRQ(ierr);
649 
        ierr = STApply(eps->OP,eps->V[i],w);CHKERRQ(ierr);
650 
        ierr = VecAXPY(w,-ritz[i],eps->V[i]);CHKERRQ(ierr);
2316 jroman 651 
        ierr = VecNorm(w,NORM_2,&norm);CHKERRQ(ierr);
2070 jroman 652 
        ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
653 
        if (bnd[i]>=eps->tol) conv[i] = 'S';
897 dsic.upv.es!antodo 654 
      }
2818 jroman 655 
      for (i=nconv;i<k;i++)
1638 slepc 656 
        if (conv[i] != 'C') {
657 
          j = i + 1;
658 
          while (j<k && conv[j] != 'C') j++;
659 
          if (j>=k) break;
660 
          /* swap eigenvalues i and j */
2798 jroman 661 
          stmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = stmp;
1638 slepc 662 
          rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
663 
          ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
664 
          /* swap eigenvectors i and j */
2818 jroman 665 
          ierr = VecSwap(eps->V[i],eps->V[j]);CHKERRQ(ierr);
1638 slepc 666 
        }
667 
      k = i;
870 dsic.upv.es!antodo 668 
    }
669 
 
1638 slepc 670 
    /* store ritz values and estimated errors */
2818 jroman 671 
    for (i=nconv;i<n;i++) {
672 
      eps->eigr[i] = ritz[i];
673 
      eps->errest[i] = bnd[i];
870 dsic.upv.es!antodo 674 
    }
2818 jroman 675 
    ierr = EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,n);CHKERRQ(ierr);
676 
    nconv = k;
1638 slepc 677 
    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
678 
    if (nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
679 
 
680 
     if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
1940 jroman 681 
      if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
1638 slepc 682 
        /* Reorthonormalize restart vector */
2316 jroman 683 
        ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,nconv,PETSC_NULL,eps->V,f,PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
684 
        ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
1638 slepc 685 
      }
686 
      if (breakdown) {
2316 jroman 687 
        /* Use random vector for restarting */
2499 jroman 688 
        ierr = PetscInfo(eps,"Using random vector for restart\n");CHKERRQ(ierr);
2316 jroman 689 
        ierr = EPSGetStartVector(eps,nconv,f,&breakdown);CHKERRQ(ierr);
1057 slepc 690 
      }
1638 slepc 691 
      if (breakdown) { /* give up */
692 
        eps->reason = EPS_DIVERGED_BREAKDOWN;
2499 jroman 693 
        ierr = PetscInfo(eps,"Unable to generate more start vectors\n");CHKERRQ(ierr);
1638 slepc 694 
      } else {
695 
        ierr = VecCopy(f,eps->V[nconv]);CHKERRQ(ierr);
696 
      }
697 
    }    
647 dsic.upv.es!antodo 698 
  }
699 
 
655 dsic.upv.es!jroman 700 
  eps->nconv = nconv;
647 dsic.upv.es!antodo 701 
 
666 dsic.upv.es!antodo 702 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 703 
  ierr = PetscFree(bnd);CHKERRQ(ierr);
1638 slepc 704 
  ierr = PetscFree(perm);CHKERRQ(ierr);
895 dsic.upv.es!antodo 705 
  ierr = PetscFree(conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 706 
  PetscFunctionReturn(0);
707 
}
708 
 
671 dsic.upv.es!antodo 709 
#undef __FUNCT__  
2324 jroman 710 
#define __FUNCT__ "EPSSetFromOptions_Lanczos"
711 
PetscErrorCode EPSSetFromOptions_Lanczos(EPS eps)
671 dsic.upv.es!antodo 712 
{
2322 jroman 713 
  PetscErrorCode         ierr;
714 
  EPS_LANCZOS            *lanczos = (EPS_LANCZOS *)eps->data;
715 
  PetscBool              flg;
716 
  EPSLanczosReorthogType reorthog;
671 dsic.upv.es!antodo 717 
 
718 
  PetscFunctionBegin;
2384 jroman 719 
  ierr = PetscOptionsHead("EPS Lanczos Options");CHKERRQ(ierr);
2322 jroman 720 
  ierr = PetscOptionsEnum("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",EPSLanczosReorthogTypes,(PetscEnum)lanczos->reorthog,(PetscEnum*)&reorthog,&flg);CHKERRQ(ierr);
721 
  if (flg) { ierr = EPSLanczosSetReorthog(eps,reorthog);CHKERRQ(ierr); }
2384 jroman 722 
  ierr = PetscOptionsTail();CHKERRQ(ierr);
671 dsic.upv.es!antodo 723 
  PetscFunctionReturn(0);
724 
}
725 
 
647 dsic.upv.es!antodo 726 
EXTERN_C_BEGIN
727 
#undef __FUNCT__  
2324 jroman 728 
#define __FUNCT__ "EPSLanczosSetReorthog_Lanczos"
729 
PetscErrorCode EPSLanczosSetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 730 
{
731 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
732 
 
733 
  PetscFunctionBegin;
734 
  switch (reorthog) {
1940 jroman 735 
    case EPS_LANCZOS_REORTHOG_LOCAL:
736 
    case EPS_LANCZOS_REORTHOG_FULL:
737 
    case EPS_LANCZOS_REORTHOG_DELAYED:
738 
    case EPS_LANCZOS_REORTHOG_SELECTIVE:
739 
    case EPS_LANCZOS_REORTHOG_PERIODIC:
740 
    case EPS_LANCZOS_REORTHOG_PARTIAL:
671 dsic.upv.es!antodo 741 
      lanczos->reorthog = reorthog;
742 
      break;
743 
    default:
2214 jroman 744 
      SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 745 
  }
746 
  PetscFunctionReturn(0);
747 
}
748 
EXTERN_C_END
749 
 
750 
#undef __FUNCT__  
906 dsic.upv.es!antodo 751 
#define __FUNCT__ "EPSLanczosSetReorthog"
671 dsic.upv.es!antodo 752 
/*@
906 dsic.upv.es!antodo 753 
   EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 754 
   iteration.
755 
 
2328 jroman 756 
   Logically Collective on EPS
671 dsic.upv.es!antodo 757 
 
758 
   Input Parameters:
759 
+  eps - the eigenproblem solver context
760 
-  reorthog - the type of reorthogonalization
761 
 
762 
   Options Database Key:
1486 slepc 763 
.  -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
764 
                         'periodic', 'partial', 'full' or 'delayed')
671 dsic.upv.es!antodo 765 
 
766 
   Level: advanced
767 
 
1364 slepc 768 
.seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 769 
@*/
906 dsic.upv.es!antodo 770 
PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 771 
{
2221 jroman 772 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 773 
 
774 
  PetscFunctionBegin;
2213 jroman 775 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 776 
  PetscValidLogicalCollectiveEnum(eps,reorthog,2);
2221 jroman 777 
  ierr = PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 778 
  PetscFunctionReturn(0);
779 
}
780 
 
781 
EXTERN_C_BEGIN
782 
#undef __FUNCT__  
2324 jroman 783 
#define __FUNCT__ "EPSLanczosGetReorthog_Lanczos"
784 
PetscErrorCode EPSLanczosGetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 785 
{
786 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
2317 jroman 787 
 
671 dsic.upv.es!antodo 788 
  PetscFunctionBegin;
789 
  *reorthog = lanczos->reorthog;
790 
  PetscFunctionReturn(0);
791 
}
792 
EXTERN_C_END
793 
 
794 
#undef __FUNCT__  
906 dsic.upv.es!antodo 795 
#define __FUNCT__ "EPSLanczosGetReorthog"
707 dsic.upv.es!antodo 796 
/*@C
906 dsic.upv.es!antodo 797 
   EPSLanczosGetReorthog - Gets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 798 
   iteration.
799 
 
2328 jroman 800 
   Not Collective
671 dsic.upv.es!antodo 801 
 
802 
   Input Parameter:
803 
.  eps - the eigenproblem solver context
804 
 
805 
   Input Parameter:
806 
.  reorthog - the type of reorthogonalization
807 
 
808 
   Level: advanced
809 
 
1364 slepc 810 
.seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 811 
@*/
906 dsic.upv.es!antodo 812 
PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 813 
{
2221 jroman 814 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 815 
 
816 
  PetscFunctionBegin;
2213 jroman 817 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 818 
  PetscValidPointer(reorthog,2);
2221 jroman 819 
  ierr = PetscTryMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 820 
  PetscFunctionReturn(0);
821 
}
822 
 
823 
#undef __FUNCT__  
2348 jroman 824 
#define __FUNCT__ "EPSReset_Lanczos"
825 
PetscErrorCode EPSReset_Lanczos(EPS eps)
826 
{
827 
  PetscErrorCode ierr;
828 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
829 
 
830 
  PetscFunctionBegin;
2410 jroman 831 
  ierr = VecDestroyVecs(eps->ncv,&lanczos->AV);CHKERRQ(ierr);
2348 jroman 832 
  ierr = EPSReset_Default(eps);CHKERRQ(ierr);
833 
  PetscFunctionReturn(0);
834 
}
835 
 
836 
#undef __FUNCT__  
2324 jroman 837 
#define __FUNCT__ "EPSDestroy_Lanczos"
838 
PetscErrorCode EPSDestroy_Lanczos(EPS eps)
1925 jroman 839 
{
840 
  PetscErrorCode ierr;
841 
 
842 
  PetscFunctionBegin;
2348 jroman 843 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
1925 jroman 844 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
845 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
846 
  PetscFunctionReturn(0);
847 
}
848 
 
849 
#undef __FUNCT__  
2324 jroman 850 
#define __FUNCT__ "EPSView_Lanczos"
851 
PetscErrorCode EPSView_Lanczos(EPS eps,PetscViewer viewer)
671 dsic.upv.es!antodo 852 
{
853 
  PetscErrorCode ierr;
854 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
2216 jroman 855 
  PetscBool      isascii;
671 dsic.upv.es!antodo 856 
 
857 
  PetscFunctionBegin;
2823 jroman 858 
  ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
2762 jroman 859 
  if (!isascii) SETERRQ1(((PetscObject)eps)->comm,1,"Viewer type %s not supported for EPS Lanczos",((PetscObject)viewer)->type_name);
2365 jroman 860 
  ierr = PetscViewerASCIIPrintf(viewer,"  Lanczos: %s reorthogonalization\n",EPSLanczosReorthogTypes[lanczos->reorthog]);CHKERRQ(ierr);
671 dsic.upv.es!antodo 861 
  PetscFunctionReturn(0);
862 
}
863 
 
864 
EXTERN_C_BEGIN
865 
#undef __FUNCT__  
2324 jroman 866 
#define __FUNCT__ "EPSCreate_Lanczos"
867 
PetscErrorCode EPSCreate_Lanczos(EPS eps)
647 dsic.upv.es!antodo 868 
{
671 dsic.upv.es!antodo 869 
  PetscErrorCode ierr;
870 
 
647 dsic.upv.es!antodo 871 
  PetscFunctionBegin;
2329 jroman 872 
  ierr = PetscNewLog(eps,EPS_LANCZOS,&eps->data);CHKERRQ(ierr);
2324 jroman 873 
  eps->ops->setup                = EPSSetUp_Lanczos;
874 
  eps->ops->setfromoptions       = EPSSetFromOptions_Lanczos;
875 
  eps->ops->destroy              = EPSDestroy_Lanczos;
2348 jroman 876 
  eps->ops->reset                = EPSReset_Lanczos;
2324 jroman 877 
  eps->ops->view                 = EPSView_Lanczos;
647 dsic.upv.es!antodo 878 
  eps->ops->backtransform        = EPSBackTransform_Default;
1928 jroman 879 
  eps->ops->computevectors       = EPSComputeVectors_Hermitian;
2324 jroman 880 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","EPSLanczosSetReorthog_Lanczos",EPSLanczosSetReorthog_Lanczos);CHKERRQ(ierr);
881 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","EPSLanczosGetReorthog_Lanczos",EPSLanczosGetReorthog_Lanczos);CHKERRQ(ierr);
647 dsic.upv.es!antodo 882 
  PetscFunctionReturn(0);
883 
}
884 
EXTERN_C_END
885