Subversion Repositories slepc-dev

Rev

Details | Compare with Previous | Last modification | View Log | RSS feed

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