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
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;
2792 jroman 529 
  PetscInt       nconv,i,j,k,l,x,n,m,*perm,restart,ncv=eps->ncv,r,ld;
1755 antodo 530 
  Vec            w=eps->work[1],f=eps->work[0];
1623 slepc 531 
  PetscScalar    *Y,stmp;
2070 jroman 532 
  PetscReal      *d,*e,*ritz,*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);
666 dsic.upv.es!antodo 538 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&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 */
1577 slepc 554 
    m = PetscMin(nconv+eps->mpd,ncv);
2792 jroman 555 
    ierr = PSGetArrayReal(eps->ps,PS_MAT_T,&d);CHKERRQ(ierr);
556 
    e = d + ld;
1616 slepc 557 
    ierr = EPSBasicLanczos(eps,d,e,eps->V,nconv,&m,f,&breakdown,anorm);CHKERRQ(ierr);
772 dsic.upv.es!antodo 558 
    n = m - nconv;
1616 slepc 559 
    beta = e[n-1];
2792 jroman 560 
    ierr = PSRestoreArrayReal(eps->ps,PS_MAT_T,&d);CHKERRQ(ierr);
561 
    ierr = PSSetDimensions(eps->ps,n,0,0);CHKERRQ(ierr);
562 
    ierr = PSSetState(eps->ps,PS_STATE_INTERMEDIATE);CHKERRQ(ierr);
563 
 
564 
    /* Solve projected problem */
565 
    ierr = PSSolve(eps->ps,ritz,PETSC_NULL);CHKERRQ(ierr);
566 
    ierr = PSSort(eps->ps,ritz,PETSC_NULL,eps->which_func,eps->which_ctx);CHKERRQ(ierr);
1638 slepc 567 
 
891 dsic.upv.es!antodo 568 
    /* Estimate ||A|| */
569 
    for (i=0;i<n;i++)
570 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
647 dsic.upv.es!antodo 571 
 
902 dsic.upv.es!antodo 572 
    /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
2792 jroman 573 
    ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 574 
    for (i=0;i<n;i++) {
2792 jroman 575 
      resnorm = beta*PetscAbsScalar(Y[n-1+i*ld]) + PETSC_MACHINE_EPSILON*anorm;
2070 jroman 576 
      ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],resnorm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
577 
      if (bnd[i]<eps->tol) {
902 dsic.upv.es!antodo 578 
        conv[i] = 'C';
1638 slepc 579 
      } else {
580 
        conv[i] = 'N';
581 
      }
687 dsic.upv.es!antodo 582 
    }
2792 jroman 583 
    ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
647 dsic.upv.es!antodo 584 
 
1638 slepc 585 
    /* purge repeated ritz values */
1940 jroman 586 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL)
1638 slepc 587 
      for (i=1;i<n;i++)
588 
        if (conv[i] == 'C')
589 
          if (PetscAbsScalar((ritz[i]-ritz[i-1])/ritz[i]) < eps->tol)
590 
            conv[i] = 'R';
902 dsic.upv.es!antodo 591 
 
1638 slepc 592 
    /* Compute restart vector */
902 dsic.upv.es!antodo 593 
    if (breakdown) {
2499 jroman 594 
      ierr = PetscInfo2(eps,"Breakdown in Lanczos method (it=%D norm=%G)\n",eps->its,beta);CHKERRQ(ierr);
1638 slepc 595 
    } else {
596 
      restart = 0;
597 
      while (restart<n && conv[restart] != 'N') restart++;
598 
      if (restart >= n) {
599 
        breakdown = PETSC_TRUE;
600 
      } else {
601 
        for (i=restart+1;i<n;i++)
602 
          if (conv[i] == 'N') {
2161 jroman 603 
            ierr = EPSCompareEigenvalues(eps,ritz[restart],0.0,ritz[i],0.0,&r);CHKERRQ(ierr);
604 
            if (r>0) restart = i;
1638 slepc 605 
          }
2792 jroman 606 
        ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
607 
        ierr = SlepcVecMAXPBY(f,0.0,1.0,n,Y+restart*ld,eps->V+nconv);CHKERRQ(ierr);
608 
        ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 609 
      }
902 dsic.upv.es!antodo 610 
    }
1638 slepc 611 
 
612 
    /* Count and put converged eigenvalues first */
613 
    for (i=0;i<n;i++) perm[i] = i;
614 
    for (k=0;k<n;k++)
615 
      if (conv[perm[k]] != 'C') {
616 
        j = k + 1;
617 
        while (j<n && conv[perm[j]] != 'C') j++;
618 
        if (j>=n) break;
619 
        l = perm[k]; perm[k] = perm[j]; perm[j] = l;
620 
      }
621 
 
622 
    /* Sort eigenvectors according to permutation */
2792 jroman 623 
    ierr = PSGetArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 624 
    for (i=0;i<k;i++) {
625 
      x = perm[i];
626 
      if (x != i) {
627 
        j = i + 1;
628 
        while (perm[j] != i) j++;
629 
        /* swap eigenvalues i and j */
630 
        rtmp = ritz[x]; ritz[x] = ritz[i]; ritz[i] = rtmp;
631 
        rtmp = bnd[x]; bnd[x] = bnd[i]; bnd[i] = rtmp;
632 
        ctmp = conv[x]; conv[x] = conv[i]; conv[i] = ctmp;
633 
        perm[j] = x; perm[i] = i;
634 
        /* swap eigenvectors i and j */
635 
        for (l=0;l<n;l++) {
2792 jroman 636 
          stmp = Y[l+x*ld]; Y[l+x*ld] = Y[l+i*ld]; Y[l+i*ld] = stmp;
1638 slepc 637 
        }
638 
      }
639 
    }
891 dsic.upv.es!antodo 640 
 
1638 slepc 641 
    /* compute converged eigenvectors */
642 
    ierr = SlepcUpdateVectors(n,eps->V+nconv,0,k,Y,n,PETSC_FALSE);CHKERRQ(ierr);
2792 jroman 643 
    ierr = PSRestoreArray(eps->ps,PS_MAT_Q,&Y);CHKERRQ(ierr);
1638 slepc 644 
 
645 
    /* purge spurious ritz values */
1940 jroman 646 
    if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL) {
1638 slepc 647 
      for (i=0;i<k;i++) {
2316 jroman 648 
        ierr = VecNorm(eps->V[nconv+i],NORM_2,&norm);CHKERRQ(ierr);
1638 slepc 649 
        ierr = VecScale(eps->V[nconv+i],1.0/norm);CHKERRQ(ierr);
650 
        ierr = STApply(eps->OP,eps->V[nconv+i],w);CHKERRQ(ierr);
2316 jroman 651 
        ierr = VecAXPY(w,-ritz[i],eps->V[nconv+i]);CHKERRQ(ierr);
652 
        ierr = VecNorm(w,NORM_2,&norm);CHKERRQ(ierr);
2070 jroman 653 
        ierr = (*eps->conv_func)(eps,ritz[i],eps->eigi[i],norm,&bnd[i],eps->conv_ctx);CHKERRQ(ierr);
654 
        if (bnd[i]>=eps->tol) conv[i] = 'S';
897 dsic.upv.es!antodo 655 
      }
1638 slepc 656 
      for (i=0;i<k;i++)
657 
        if (conv[i] != 'C') {
658 
          j = i + 1;
659 
          while (j<k && conv[j] != 'C') j++;
660 
          if (j>=k) break;
661 
          /* swap eigenvalues i and j */
662 
          rtmp = ritz[j]; ritz[j] = ritz[i]; ritz[i] = rtmp;
663 
          rtmp = bnd[j]; bnd[j] = bnd[i]; bnd[i] = rtmp;
664 
          ctmp = conv[j]; conv[j] = conv[i]; conv[i] = ctmp;
665 
          /* swap eigenvectors i and j */
666 
          ierr = VecSwap(eps->V[nconv+i],eps->V[nconv+j]);CHKERRQ(ierr);
667 
        }
668 
      k = i;
870 dsic.upv.es!antodo 669 
    }
670 
 
1638 slepc 671 
    /* store ritz values and estimated errors */
672 
    for (i=0;i<n;i++) {
673 
      eps->eigr[nconv+i] = ritz[i];
674 
      eps->errest[nconv+i] = bnd[i];
870 dsic.upv.es!antodo 675 
    }
2313 jroman 676 
    ierr = EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nconv+n);CHKERRQ(ierr);
1638 slepc 677 
    nconv = nconv+k;
678 
    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
679 
    if (nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
680 
 
681 
     if (eps->reason == EPS_CONVERGED_ITERATING) { /* copy restart vector */
1940 jroman 682 
      if (lanczos->reorthog == EPS_LANCZOS_REORTHOG_LOCAL && !breakdown) {
1638 slepc 683 
        /* Reorthonormalize restart vector */
2316 jroman 684 
        ierr = IPOrthogonalize(eps->ip,eps->nds,eps->DS,nconv,PETSC_NULL,eps->V,f,PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
685 
        ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
1638 slepc 686 
      }
687 
      if (breakdown) {
2316 jroman 688 
        /* Use random vector for restarting */
2499 jroman 689 
        ierr = PetscInfo(eps,"Using random vector for restart\n");CHKERRQ(ierr);
2316 jroman 690 
        ierr = EPSGetStartVector(eps,nconv,f,&breakdown);CHKERRQ(ierr);
1057 slepc 691 
      }
1638 slepc 692 
      if (breakdown) { /* give up */
693 
        eps->reason = EPS_DIVERGED_BREAKDOWN;
2499 jroman 694 
        ierr = PetscInfo(eps,"Unable to generate more start vectors\n");CHKERRQ(ierr);
1638 slepc 695 
      } else {
696 
        ierr = VecCopy(f,eps->V[nconv]);CHKERRQ(ierr);
697 
      }
698 
    }    
647 dsic.upv.es!antodo 699 
  }
700 
 
655 dsic.upv.es!jroman 701 
  eps->nconv = nconv;
647 dsic.upv.es!antodo 702 
 
666 dsic.upv.es!antodo 703 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 704 
  ierr = PetscFree(bnd);CHKERRQ(ierr);
1638 slepc 705 
  ierr = PetscFree(perm);CHKERRQ(ierr);
895 dsic.upv.es!antodo 706 
  ierr = PetscFree(conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 707 
  PetscFunctionReturn(0);
708 
}
709 
 
671 dsic.upv.es!antodo 710 
#undef __FUNCT__  
2324 jroman 711 
#define __FUNCT__ "EPSSetFromOptions_Lanczos"
712 
PetscErrorCode EPSSetFromOptions_Lanczos(EPS eps)
671 dsic.upv.es!antodo 713 
{
2322 jroman 714 
  PetscErrorCode         ierr;
715 
  EPS_LANCZOS            *lanczos = (EPS_LANCZOS *)eps->data;
716 
  PetscBool              flg;
717 
  EPSLanczosReorthogType reorthog;
671 dsic.upv.es!antodo 718 
 
719 
  PetscFunctionBegin;
2384 jroman 720 
  ierr = PetscOptionsHead("EPS Lanczos Options");CHKERRQ(ierr);
2322 jroman 721 
  ierr = PetscOptionsEnum("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",EPSLanczosReorthogTypes,(PetscEnum)lanczos->reorthog,(PetscEnum*)&reorthog,&flg);CHKERRQ(ierr);
722 
  if (flg) { ierr = EPSLanczosSetReorthog(eps,reorthog);CHKERRQ(ierr); }
2384 jroman 723 
  ierr = PetscOptionsTail();CHKERRQ(ierr);
671 dsic.upv.es!antodo 724 
  PetscFunctionReturn(0);
725 
}
726 
 
647 dsic.upv.es!antodo 727 
EXTERN_C_BEGIN
728 
#undef __FUNCT__  
2324 jroman 729 
#define __FUNCT__ "EPSLanczosSetReorthog_Lanczos"
730 
PetscErrorCode EPSLanczosSetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 731 
{
732 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
733 
 
734 
  PetscFunctionBegin;
735 
  switch (reorthog) {
1940 jroman 736 
    case EPS_LANCZOS_REORTHOG_LOCAL:
737 
    case EPS_LANCZOS_REORTHOG_FULL:
738 
    case EPS_LANCZOS_REORTHOG_DELAYED:
739 
    case EPS_LANCZOS_REORTHOG_SELECTIVE:
740 
    case EPS_LANCZOS_REORTHOG_PERIODIC:
741 
    case EPS_LANCZOS_REORTHOG_PARTIAL:
671 dsic.upv.es!antodo 742 
      lanczos->reorthog = reorthog;
743 
      break;
744 
    default:
2214 jroman 745 
      SETERRQ(((PetscObject)eps)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 746 
  }
747 
  PetscFunctionReturn(0);
748 
}
749 
EXTERN_C_END
750 
 
751 
#undef __FUNCT__  
906 dsic.upv.es!antodo 752 
#define __FUNCT__ "EPSLanczosSetReorthog"
671 dsic.upv.es!antodo 753 
/*@
906 dsic.upv.es!antodo 754 
   EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 755 
   iteration.
756 
 
2328 jroman 757 
   Logically Collective on EPS
671 dsic.upv.es!antodo 758 
 
759 
   Input Parameters:
760 
+  eps - the eigenproblem solver context
761 
-  reorthog - the type of reorthogonalization
762 
 
763 
   Options Database Key:
1486 slepc 764 
.  -eps_lanczos_reorthog - Sets the reorthogonalization type (either 'local', 'selective',
765 
                         'periodic', 'partial', 'full' or 'delayed')
671 dsic.upv.es!antodo 766 
 
767 
   Level: advanced
768 
 
1364 slepc 769 
.seealso: EPSLanczosGetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 770 
@*/
906 dsic.upv.es!antodo 771 
PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 772 
{
2221 jroman 773 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 774 
 
775 
  PetscFunctionBegin;
2213 jroman 776 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 777 
  PetscValidLogicalCollectiveEnum(eps,reorthog,2);
2221 jroman 778 
  ierr = PetscTryMethod(eps,"EPSLanczosSetReorthog_C",(EPS,EPSLanczosReorthogType),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 779 
  PetscFunctionReturn(0);
780 
}
781 
 
782 
EXTERN_C_BEGIN
783 
#undef __FUNCT__  
2324 jroman 784 
#define __FUNCT__ "EPSLanczosGetReorthog_Lanczos"
785 
PetscErrorCode EPSLanczosGetReorthog_Lanczos(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 786 
{
787 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
2317 jroman 788 
 
671 dsic.upv.es!antodo 789 
  PetscFunctionBegin;
790 
  *reorthog = lanczos->reorthog;
791 
  PetscFunctionReturn(0);
792 
}
793 
EXTERN_C_END
794 
 
795 
#undef __FUNCT__  
906 dsic.upv.es!antodo 796 
#define __FUNCT__ "EPSLanczosGetReorthog"
707 dsic.upv.es!antodo 797 
/*@C
906 dsic.upv.es!antodo 798 
   EPSLanczosGetReorthog - Gets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 799 
   iteration.
800 
 
2328 jroman 801 
   Not Collective
671 dsic.upv.es!antodo 802 
 
803 
   Input Parameter:
804 
.  eps - the eigenproblem solver context
805 
 
806 
   Input Parameter:
807 
.  reorthog - the type of reorthogonalization
808 
 
809 
   Level: advanced
810 
 
1364 slepc 811 
.seealso: EPSLanczosSetReorthog(), EPSLanczosReorthogType
671 dsic.upv.es!antodo 812 
@*/
906 dsic.upv.es!antodo 813 
PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 814 
{
2221 jroman 815 
  PetscErrorCode ierr;
671 dsic.upv.es!antodo 816 
 
817 
  PetscFunctionBegin;
2213 jroman 818 
  PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
2326 jroman 819 
  PetscValidPointer(reorthog,2);
2221 jroman 820 
  ierr = PetscTryMethod(eps,"EPSLanczosGetReorthog_C",(EPS,EPSLanczosReorthogType*),(eps,reorthog));CHKERRQ(ierr);
671 dsic.upv.es!antodo 821 
  PetscFunctionReturn(0);
822 
}
823 
 
824 
#undef __FUNCT__  
2348 jroman 825 
#define __FUNCT__ "EPSReset_Lanczos"
826 
PetscErrorCode EPSReset_Lanczos(EPS eps)
827 
{
828 
  PetscErrorCode ierr;
829 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
830 
 
831 
  PetscFunctionBegin;
2410 jroman 832 
  ierr = VecDestroyVecs(eps->ncv,&lanczos->AV);CHKERRQ(ierr);
2348 jroman 833 
  ierr = EPSReset_Default(eps);CHKERRQ(ierr);
834 
  PetscFunctionReturn(0);
835 
}
836 
 
837 
#undef __FUNCT__  
2324 jroman 838 
#define __FUNCT__ "EPSDestroy_Lanczos"
839 
PetscErrorCode EPSDestroy_Lanczos(EPS eps)
1925 jroman 840 
{
841 
  PetscErrorCode ierr;
842 
 
843 
  PetscFunctionBegin;
2348 jroman 844 
  ierr = PetscFree(eps->data);CHKERRQ(ierr);
1925 jroman 845 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
846 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","",PETSC_NULL);CHKERRQ(ierr);
847 
  PetscFunctionReturn(0);
848 
}
849 
 
850 
#undef __FUNCT__  
2324 jroman 851 
#define __FUNCT__ "EPSView_Lanczos"
852 
PetscErrorCode EPSView_Lanczos(EPS eps,PetscViewer viewer)
671 dsic.upv.es!antodo 853 
{
854 
  PetscErrorCode ierr;
855 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
2216 jroman 856 
  PetscBool      isascii;
671 dsic.upv.es!antodo 857 
 
858 
  PetscFunctionBegin;
2215 jroman 859 
  ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr);
2762 jroman 860 
  if (!isascii) SETERRQ1(((PetscObject)eps)->comm,1,"Viewer type %s not supported for EPS Lanczos",((PetscObject)viewer)->type_name);
2365 jroman 861 
  ierr = PetscViewerASCIIPrintf(viewer,"  Lanczos: %s reorthogonalization\n",EPSLanczosReorthogTypes[lanczos->reorthog]);CHKERRQ(ierr);
671 dsic.upv.es!antodo 862 
  PetscFunctionReturn(0);
863 
}
864 
 
865 
EXTERN_C_BEGIN
866 
#undef __FUNCT__  
2324 jroman 867 
#define __FUNCT__ "EPSCreate_Lanczos"
868 
PetscErrorCode EPSCreate_Lanczos(EPS eps)
647 dsic.upv.es!antodo 869 
{
671 dsic.upv.es!antodo 870 
  PetscErrorCode ierr;
871 
 
647 dsic.upv.es!antodo 872 
  PetscFunctionBegin;
2329 jroman 873 
  ierr = PetscNewLog(eps,EPS_LANCZOS,&eps->data);CHKERRQ(ierr);
2324 jroman 874 
  eps->ops->setup                = EPSSetUp_Lanczos;
875 
  eps->ops->setfromoptions       = EPSSetFromOptions_Lanczos;
876 
  eps->ops->destroy              = EPSDestroy_Lanczos;
2348 jroman 877 
  eps->ops->reset                = EPSReset_Lanczos;
2324 jroman 878 
  eps->ops->view                 = EPSView_Lanczos;
647 dsic.upv.es!antodo 879 
  eps->ops->backtransform        = EPSBackTransform_Default;
1928 jroman 880 
  eps->ops->computevectors       = EPSComputeVectors_Hermitian;
2324 jroman 881 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","EPSLanczosSetReorthog_Lanczos",EPSLanczosSetReorthog_Lanczos);CHKERRQ(ierr);
882 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","EPSLanczosGetReorthog_Lanczos",EPSLanczosGetReorthog_Lanczos);CHKERRQ(ierr);
647 dsic.upv.es!antodo 883 
  PetscFunctionReturn(0);
884 
}
885 
EXTERN_C_END
886