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