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