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 
 
647 dsic.upv.es!antodo 20 
*/
707 dsic.upv.es!antodo 21 
#include "src/eps/epsimpl.h"                /*I "slepceps.h" I*/
647 dsic.upv.es!antodo 22 
#include "slepcblaslapack.h"
23 
 
671 dsic.upv.es!antodo 24 
typedef struct {
906 dsic.upv.es!antodo 25 
  EPSLanczosReorthogType reorthog;
671 dsic.upv.es!antodo 26 
} EPS_LANCZOS;
27 
 
647 dsic.upv.es!antodo 28 
#undef __FUNCT__  
29 
#define __FUNCT__ "EPSSetUp_LANCZOS"
30 
PetscErrorCode EPSSetUp_LANCZOS(EPS eps)
31 
{
32 
  PetscErrorCode ierr;
982 slepc 33 
  PetscInt       N;
647 dsic.upv.es!antodo 34 
 
35 
  PetscFunctionBegin;
36 
  ierr = VecGetSize(eps->vec_initial,&N);CHKERRQ(ierr);
37 
  if (eps->ncv) {
38 
    if (eps->ncv<eps->nev) SETERRQ(1,"The value of ncv must be at least nev");
39 
  }
651 dsic.upv.es!antodo 40 
  else eps->ncv = PetscMin(N,PetscMax(2*eps->nev,eps->nev+15));
1220 slepc 41 
  if (!eps->max_it) eps->max_it = PetscMax(100,2*N/eps->ncv);
42 
 
906 dsic.upv.es!antodo 43 
  if (eps->solverclass==EPS_ONE_SIDE) {
44 
    if (eps->which == EPS_LARGEST_IMAGINARY || eps->which == EPS_SMALLEST_IMAGINARY)
45 
      SETERRQ(1,"Wrong value of eps->which");
46 
    if (!eps->ishermitian)
47 
      SETERRQ(PETSC_ERR_SUP,"Requested method is only available for Hermitian problems");
48 
  } else {
49 
    if (eps->which != EPS_LARGEST_MAGNITUDE)
50 
      SETERRQ(1,"Wrong value of eps->which");
51 
  }
647 dsic.upv.es!antodo 52 
  ierr = EPSAllocateSolution(eps);CHKERRQ(ierr);
1040 slepc 53 
  ierr = PetscFree(eps->T);CHKERRQ(ierr);
647 dsic.upv.es!antodo 54 
  ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&eps->T);CHKERRQ(ierr);
885 dsic.upv.es!jroman 55 
  if (eps->solverclass==EPS_TWO_SIDE) {
1040 slepc 56 
    ierr = PetscFree(eps->Tl);CHKERRQ(ierr);
885 dsic.upv.es!jroman 57 
    ierr = PetscMalloc(eps->ncv*eps->ncv*sizeof(PetscScalar),&eps->Tl);CHKERRQ(ierr);
58 
    ierr = EPSDefaultGetWork(eps,2);CHKERRQ(ierr);
59 
  }
60 
  else { ierr = EPSDefaultGetWork(eps,1);CHKERRQ(ierr); }
647 dsic.upv.es!antodo 61 
  PetscFunctionReturn(0);
62 
}
63 
 
64 
#undef __FUNCT__  
1068 slepc 65 
#define __FUNCT__ "EPSLocalLanczos"
683 dsic.upv.es!jroman 66 
/*
1068 slepc 67 
   EPSLocalLanczos - Local reorthogonalization.
683 dsic.upv.es!jroman 68 
 
69 
   This is the simplest variant. At each Lanczos step, the corresponding Lanczos vector
70 
   is orthogonalized with respect to the two previous Lanczos vectors, according to
71 
   the three term Lanczos recurrence. WARNING: This variant does not track the loss of
72 
   orthogonality that occurs in finite-precision arithmetic and, therefore, the
73 
   generated vectors are not guaranteed to be (semi-)orthogonal.
74 
*/
1220 slepc 75 
static PetscErrorCode EPSLocalLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta,PetscTruth *breakdown)
671 dsic.upv.es!antodo 76 
{
77 
  PetscErrorCode ierr;
1163 slepc 78 
  int            i,j,m = *M;
1133 slepc 79 
  PetscReal      norm;
1163 slepc 80 
  PetscTruth     *which,lwhich[100];
1121 slepc 81 
 
1163 slepc 82 
  PetscFunctionBegin;  
83 
  if (m>100) {
84 
    ierr = PetscMalloc(sizeof(PetscTruth)*m,&which);CHKERRQ(ierr);
85 
  } else which = lwhich;
86 
  for (i=0;i<k;i++)
87 
    which[i] = PETSC_TRUE;
1121 slepc 88 
 
89 
  for (j=k;j<m;j++) {
90 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1179 slepc 91 
    ierr = EPSOrthogonalize(eps,eps->nds,PETSC_NULL,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
1163 slepc 92 
    which[j] = PETSC_TRUE;
93 
    if (j-2>=k) which[j-2] = PETSC_FALSE;
94 
    ierr = EPSOrthogonalize(eps,j+1,which,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);
1121 slepc 95 
    if (*breakdown) {
96 
      *M = j+1;
97 
      break;
98 
    }
99 
    if (j<m-1) {
100 
      T[m*j+j+1] = norm;
101 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
102 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
103 
    }
104 
  }
1133 slepc 105 
  *beta = norm;
106 
 
1163 slepc 107 
  if (m>100) { ierr = PetscFree(which);CHKERRQ(ierr); }
1121 slepc 108 
  PetscFunctionReturn(0);
109 
}
110 
 
111 
#undef __FUNCT__  
891 dsic.upv.es!antodo 112 
#define __FUNCT__ "EPSSelectiveLanczos"
928 dsic.upv.es!antodo 113 
/*
114 
   EPSSelectiveLanczos - Selective reorthogonalization.
115 
*/
891 dsic.upv.es!antodo 116 
static PetscErrorCode EPSSelectiveLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta,PetscTruth *breakdown,PetscReal anorm)
117 
{
118 
  PetscErrorCode ierr;
1178 slepc 119 
  int            i,j,m = *M,n,nritz=0,nritzo;
120 
  PetscReal      *ritz,norm;
121 
  PetscScalar    *Y;
1163 slepc 122 
  PetscTruth     *which,lwhich[100];
891 dsic.upv.es!antodo 123 
 
124 
  PetscFunctionBegin;
125 
  ierr = PetscMalloc(m*sizeof(PetscReal),&ritz);CHKERRQ(ierr);
1178 slepc 126 
  ierr = PetscMalloc(m*m*sizeof(PetscScalar),&Y);CHKERRQ(ierr);
891 dsic.upv.es!antodo 127 
 
1163 slepc 128 
  if (m>100) {
129 
    ierr = PetscMalloc(sizeof(PetscTruth)*m,&which);CHKERRQ(ierr);
130 
  } else which = lwhich;
131 
  for (i=0;i<k;i++)
132 
    which[i] = PETSC_TRUE;
1124 slepc 133 
 
891 dsic.upv.es!antodo 134 
  for (j=k;j<m;j++) {
1154 slepc 135 
    /* Lanczos step */
891 dsic.upv.es!antodo 136 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1179 slepc 137 
    ierr = EPSOrthogonalize(eps,eps->nds,PETSC_NULL,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
1163 slepc 138 
    which[j] = PETSC_TRUE;
139 
    if (j-2>=k) which[j-2] = PETSC_FALSE;
140 
    ierr = EPSOrthogonalize(eps,j+1,which,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);
891 dsic.upv.es!antodo 141 
    if (*breakdown) {
142 
      *M = j+1;
143 
      break;
144 
    }
145 
 
1154 slepc 146 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
891 dsic.upv.es!antodo 147 
    n = j-k+1;
1178 slepc 148 
    ierr = EPSDenseTridiagonal(n,T+k*(m+1),m,ritz,Y);CHKERRQ(ierr);
891 dsic.upv.es!antodo 149 
 
150 
    /* Estimate ||A|| */
151 
    for (i=0;i<n;i++)
152 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
1154 slepc 153 
 
154 
    /* Compute nearly converged Ritz vectors */
155 
    nritzo = 0;
156 
    for (i=0;i<n;i++)
891 dsic.upv.es!antodo 157 
      if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm)
1154 slepc 158 
        nritzo++;
159 
 
160 
    if (nritzo>nritz) {
161 
      nritz = 0;
162 
      for (i=0;i<n;i++) {
163 
        if (norm*PetscAbsScalar(Y[i*n+n-1]) < PETSC_SQRT_MACHINE_EPSILON*anorm) {
164 
          ierr = VecSet(eps->AV[nritz],0.0);CHKERRQ(ierr);
165 
          ierr = VecMAXPY(eps->AV[nritz],n,Y+i*n,V+k);CHKERRQ(ierr);
166 
          nritz++;
167 
        }
168 
      }
891 dsic.upv.es!antodo 169 
    }
170 
 
1154 slepc 171 
    if (nritz > 0) {
1163 slepc 172 
      ierr = EPSOrthogonalize(eps,nritz,PETSC_NULL,eps->AV,f,PETSC_NULL,&norm,breakdown);CHKERRQ(ierr);
1154 slepc 173 
      if (*breakdown) {
174 
        *M = j+1;
175 
        break;
176 
      }
177 
    }
178 
 
891 dsic.upv.es!antodo 179 
    if (j<m-1) {
180 
      T[m*j+j+1] = norm;
181 
      ierr = VecScale(f,1.0 / norm);CHKERRQ(ierr);
182 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
183 
    }
671 dsic.upv.es!antodo 184 
  }
891 dsic.upv.es!antodo 185 
  *beta = norm;
1133 slepc 186 
 
891 dsic.upv.es!antodo 187 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
895 dsic.upv.es!antodo 188 
  ierr = PetscFree(Y);CHKERRQ(ierr);
1163 slepc 189 
  if (m>100) { ierr = PetscFree(which);CHKERRQ(ierr); }
671 dsic.upv.es!antodo 190 
  PetscFunctionReturn(0);
191 
}
192 
 
193 
#undef __FUNCT__  
1126 slepc 194 
#define __FUNCT__ "update_omega"
1178 slepc 195 
static void update_omega(PetscReal *omega,PetscReal *omega_old,int j,PetscReal *alpha,PetscReal *beta,PetscReal eps1,PetscReal anorm)
772 dsic.upv.es!antodo 196 
{
1126 slepc 197 
  int            k;
198 
  PetscReal      T,binv,temp;
1121 slepc 199 
 
200 
  PetscFunctionBegin;
1126 slepc 201 
  /* Estimate of contribution to roundoff errors from A*v
202 
       fl(A*v) = A*v + f,
203 
     where ||f|| \approx eps1*||A||.
204 
     For a full matrix A, a rule-of-thumb estimate is eps1 = sqrt(n)*eps. */
205 
  T = eps1*anorm;
206 
  binv = 1.0/beta[j+1];
207 
 
208 
  /* Update omega(1) using omega(0)==0. */
1178 slepc 209 
  omega_old[0]= beta[1]*omega[1] + (alpha[0]-alpha[j])*omega[0] -
1126 slepc 210 
                beta[j]*omega_old[0];
211 
  if (omega_old[0] > 0)
212 
    omega_old[0] = binv*(omega_old[0] + T);
213 
  else
214 
    omega_old[0] = binv*(omega_old[0] - T);  
215 
 
216 
  /* Update remaining components. */
217 
  for (k=1;k<j-1;k++) {
1178 slepc 218 
    omega_old[k] = beta[k+1]*omega[k+1] + (alpha[k]-alpha[j])*omega[k] +
1126 slepc 219 
                   beta[k]*omega[k-1] - beta[j]*omega_old[k];
220 
    if (omega_old[k] > 0)
221 
      omega_old[k] = binv*(omega_old[k] + T);      
222 
    else
223 
      omega_old[k] = binv*(omega_old[k] - T);      
1121 slepc 224 
  }
1126 slepc 225 
  omega_old[j-1] = binv*T;
1121 slepc 226 
 
1126 slepc 227 
  /* Swap omega and omega_old. */
228 
  for (k=0;k<j;k++) {
229 
    temp = omega[k];
230 
    omega[k] = omega_old[k];
231 
    omega_old[k] = omega[k];
1121 slepc 232 
  }
1126 slepc 233 
  omega[j] = eps1;  
1128 slepc 234 
  PetscFunctionReturnVoid();
1121 slepc 235 
}
236 
 
237 
#undef __FUNCT__  
1128 slepc 238 
#define __FUNCT__ "compute_int"
239 
static void compute_int(PetscTruth *which,PetscReal *mu,int j,PetscReal delta,PetscReal eta)
240 
{
241 
  int        i,k,maxpos;
242 
  PetscReal  max;
243 
  PetscTruth found;
244 
 
1129 slepc 245 
  PetscFunctionBegin;  
1128 slepc 246 
  /* initialize which */
247 
  found = PETSC_FALSE;
248 
  max = 0.0;
249 
  for (i=0;i<j;i++) {
1178 slepc 250 
    if (PetscAbsReal(mu[i]) >= delta) {
1128 slepc 251 
      which[i] = PETSC_TRUE;
252 
      found = PETSC_TRUE;
253 
    } else which[i] = PETSC_FALSE;
1178 slepc 254 
    if (PetscAbsReal(mu[i]) > max) {
1128 slepc 255 
      maxpos = i;
1178 slepc 256 
      max = PetscAbsReal(mu[i]);
1128 slepc 257 
    }
258 
  }
259 
  if (!found) which[maxpos] = PETSC_TRUE;    
260 
 
261 
  for (i=0;i<j;i++)
262 
    if (which[i]) {
263 
      /* find left interval */
264 
      for (k=i;k>=0;k--) {
1178 slepc 265 
        if (PetscAbsReal(mu[k])<eta || which[k]) break;
1128 slepc 266 
        else which[k] = PETSC_TRUE;
267 
      }
268 
      /* find right interval */
269 
      for (k=i+1;k<j;k++) {
1178 slepc 270 
        if (PetscAbsReal(mu[k])<eta || which[k]) break;
1128 slepc 271 
        else which[k] = PETSC_TRUE;
272 
      }
273 
    }
274 
  PetscFunctionReturnVoid();
275 
}
276 
 
277 
#undef __FUNCT__  
891 dsic.upv.es!antodo 278 
#define __FUNCT__ "EPSPartialLanczos"
928 dsic.upv.es!antodo 279 
/*
280 
   EPSPartialLanczos - Partial reorthogonalization.
281 
*/
1126 slepc 282 
static PetscErrorCode EPSPartialLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *M,Vec f,PetscReal *beta, PetscTruth *breakdown,PetscReal anorm)
746 dsic.upv.es!antodo 283 
{
1126 slepc 284 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
746 dsic.upv.es!antodo 285 
  PetscErrorCode ierr;
1126 slepc 286 
  Mat            A;
1204 slepc 287 
  int            i,j,m = *M;
288 
  PetscInt       n;
1178 slepc 289 
  PetscReal      norm,*omega,lomega[100],*omega_old,lomega_old[100],eps1,delta,eta,*b,lb[101],*a,la[100];
1163 slepc 290 
  PetscTruth     *which,lwhich[100],*which2,lwhich2[100],
291 
                 reorth = PETSC_FALSE,force_reorth = PETSC_FALSE,fro = PETSC_FALSE,estimate_anorm = PETSC_FALSE;
746 dsic.upv.es!antodo 292 
 
293 
  PetscFunctionBegin;
1163 slepc 294 
  if (m>100) {
1178 slepc 295 
    ierr = PetscMalloc(m*sizeof(PetscReal),&a);CHKERRQ(ierr);
1173 slepc 296 
    ierr = PetscMalloc((m+1)*sizeof(PetscReal),&b);CHKERRQ(ierr);
1178 slepc 297 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega);CHKERRQ(ierr);
298 
    ierr = PetscMalloc(m*sizeof(PetscReal),&omega_old);CHKERRQ(ierr);
1163 slepc 299 
    ierr = PetscMalloc(sizeof(PetscTruth)*m,&which);CHKERRQ(ierr);
300 
    ierr = PetscMalloc(sizeof(PetscTruth)*m,&which2);CHKERRQ(ierr);
1128 slepc 301 
  } else {
302 
    a = la;
303 
    b = lb;
304 
    omega = lomega;
305 
    omega_old = lomega_old;
306 
    which = lwhich;
1163 slepc 307 
    which2 = lwhich2;
1128 slepc 308 
  }
1126 slepc 309 
 
310 
  ierr = STGetOperators(eps->OP,&A,PETSC_NULL);CHKERRQ(ierr);
311 
  ierr = MatGetSize(A,&n,PETSC_NULL);CHKERRQ(ierr);
312 
  eps1 = sqrt((PetscReal)n)*PETSC_MACHINE_EPSILON/2;
313 
  delta = PETSC_SQRT_MACHINE_EPSILON/sqrt((PetscReal)eps->ncv);
1128 slepc 314 
  eta = pow(PETSC_MACHINE_EPSILON,3.0/4.0)/sqrt((PetscReal)eps->ncv);
1143 slepc 315 
  if (anorm < 0.0) {
316 
    anorm = 1.0;
317 
    estimate_anorm = PETSC_TRUE;
318 
  }
1163 slepc 319 
  for (i=0;i<m-k;i++)
1152 slepc 320 
    omega[i] = omega_old[i] = 0.0;
1163 slepc 321 
  for (i=0;i<k;i++)
322 
    which[i] = PETSC_TRUE;  
1126 slepc 323 
 
891 dsic.upv.es!antodo 324 
  for (j=k;j<m;j++) {
325 
    ierr = STApply(eps->OP,V[j],f);CHKERRQ(ierr);
1163 slepc 326 
    ierr = EPSOrthogonalize(eps,eps->nds,PETSC_NULL,eps->DS,f,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
1128 slepc 327 
    if (fro) {
328 
      /* Lanczos step with full reorthogonalization */
1163 slepc 329 
      ierr = EPSOrthogonalize(eps,j+1,PETSC_NULL,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);      
1128 slepc 330 
    } else {
331 
      /* Lanczos step */
1163 slepc 332 
      which[j] = PETSC_TRUE;
333 
      if (j-2>=k) which[j-2] = PETSC_FALSE;
334 
      ierr = EPSOrthogonalize(eps,j+1,which,V,f,T+m*j,&norm,breakdown);CHKERRQ(ierr);
1178 slepc 335 
      a[j-k] = PetscRealPart(T[m*j+j]);
1151 slepc 336 
      b[j-k+1] = norm;
1149 slepc 337 
 
1154 slepc 338 
      /* Estimate ||A|| if needed */
1152 slepc 339 
      if (estimate_anorm) {
1178 slepc 340 
        if (j>k) anorm = PetscMax(anorm,PetscAbsReal(a[j-k])+norm+b[j-k]);
341 
        else anorm = PetscMax(anorm,PetscAbsReal(a[j-k])+norm);
1152 slepc 342 
      }
1154 slepc 343 
 
344 
      /* Check if reorthogonalization is needed */
1128 slepc 345 
      reorth = PETSC_FALSE;
1154 slepc 346 
      if (j>k) {      
347 
        update_omega(omega,omega_old,j-k,a,b,eps1,anorm);
348 
        for (i=0;i<j-k;i++)
349 
          if (PetscAbsScalar(omega[i]) > delta) reorth = PETSC_TRUE;
350 
      }
1128 slepc 351 
 
352 
      if (reorth || force_reorth) {
353 
        if (lanczos->reorthog == EPSLANCZOS_REORTHOG_PERIODIC) {
354 
          /* Periodic reorthogonalization */
1129 slepc 355 
          if (force_reorth) force_reorth = PETSC_FALSE;
356 
          else force_reorth = PETSC_TRUE;
1163 slepc 357 
          ierr = EPSOrthogonalize(eps,j-k,PETSC_NULL,V+k,f,PETSC_NULL,&norm,breakdown);CHKERRQ(ierr);
1128 slepc 358 
          for (i=0;i<j-k;i++)
359 
            omega[i] = eps1;
360 
        } else {
361 
          /* Partial reorthogonalization */
1129 slepc 362 
          if (force_reorth) force_reorth = PETSC_FALSE;
363 
          else {
364 
            force_reorth = PETSC_TRUE;
1163 slepc 365 
            compute_int(which2,omega,j-k,delta,eta);
1129 slepc 366 
            for (i=0;i<j-k;i++)
1163 slepc 367 
              if (which2[i]) omega[i] = eps1;
1150 slepc 368 
          }
1163 slepc 369 
          ierr = EPSOrthogonalize(eps,j-k,which2,V+k,f,PETSC_NULL,&norm,breakdown);CHKERRQ(ierr);      
1129 slepc 370 
        }      
1126 slepc 371 
      }
746 dsic.upv.es!antodo 372 
    }
1150 slepc 373 
 
1128 slepc 374 
    if (*breakdown || norm < n*anorm*PETSC_MACHINE_EPSILON) {
891 dsic.upv.es!antodo 375 
      *M = j+1;
376 
      break;
377 
    }
1128 slepc 378 
    if (!fro && norm*delta < anorm*eps1) {
379 
      fro = PETSC_TRUE;
380 
      PetscInfo1(eps,"Switching to full reorthogonalization at iteration %i\n",eps->its);      
381 
    }
891 dsic.upv.es!antodo 382 
    if (j<m-1) {
1126 slepc 383 
      T[m*j+j+1] = b[j-k+1] = norm;
891 dsic.upv.es!antodo 384 
      ierr = VecScale(f,1.0/norm);CHKERRQ(ierr);
385 
      ierr = VecCopy(f,V[j+1]);CHKERRQ(ierr);
746 dsic.upv.es!antodo 386 
    }
387 
  }
891 dsic.upv.es!antodo 388 
  *beta = norm;
1126 slepc 389 
 
1163 slepc 390 
  if (m>100) {
1128 slepc 391 
    ierr = PetscFree(a);CHKERRQ(ierr);
392 
    ierr = PetscFree(b);CHKERRQ(ierr);
393 
    ierr = PetscFree(omega);CHKERRQ(ierr);
394 
    ierr = PetscFree(omega_old);CHKERRQ(ierr);
395 
    ierr = PetscFree(which);CHKERRQ(ierr);
1163 slepc 396 
    ierr = PetscFree(which2);CHKERRQ(ierr);
1128 slepc 397 
  }
746 dsic.upv.es!antodo 398 
  PetscFunctionReturn(0);
399 
}
400 
 
401 
#undef __FUNCT__  
655 dsic.upv.es!jroman 402 
#define __FUNCT__ "EPSBasicLanczos"
403 
/*
404 
   EPSBasicLanczos - Computes an m-step Lanczos factorization. The first k
405 
   columns are assumed to be locked and therefore they are not modified. On
406 
   exit, the following relation is satisfied:
647 dsic.upv.es!antodo 407 
 
655 dsic.upv.es!jroman 408 
                    OP * V - V * T = f * e_m^T
409 
 
410 
   where the columns of V are the Lanczos vectors, T is a tridiagonal matrix,
411 
   f is the residual vector and e_m is the m-th vector of the canonical basis.
412 
   The Lanczos vectors (together with vector f) are B-orthogonal (to working
413 
   accuracy) if full reorthogonalization is being used, otherwise they are
414 
   (B-)semi-orthogonal. On exit, beta contains the B-norm of f and the next
415 
   Lanczos vector can be computed as v_{m+1} = f / beta.
416 
 
417 
   This function simply calls another function which depends on the selected
418 
   reorthogonalization strategy.
419 
*/
891 dsic.upv.es!antodo 420 
static PetscErrorCode EPSBasicLanczos(EPS eps,PetscScalar *T,Vec *V,int k,int *m,Vec f,PetscReal *beta,PetscTruth *breakdown,PetscReal anorm)
655 dsic.upv.es!jroman 421 
{
671 dsic.upv.es!antodo 422 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
655 dsic.upv.es!jroman 423 
  PetscErrorCode ierr;
424 
 
425 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 426 
  switch (lanczos->reorthog) {
1068 slepc 427 
    case EPSLANCZOS_REORTHOG_LOCAL:
1160 slepc 428 
      ierr = EPSLocalLanczos(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
671 dsic.upv.es!antodo 429 
      break;
961 dsic.upv.es!jroman 430 
    case EPSLANCZOS_REORTHOG_SELECTIVE:
891 dsic.upv.es!antodo 431 
      ierr = EPSSelectiveLanczos(eps,T,V,k,m,f,beta,breakdown,anorm);CHKERRQ(ierr);
671 dsic.upv.es!antodo 432 
      break;
961 dsic.upv.es!jroman 433 
    case EPSLANCZOS_REORTHOG_PARTIAL:
434 
    case EPSLANCZOS_REORTHOG_PERIODIC:
1126 slepc 435 
      ierr = EPSPartialLanczos(eps,T,V,k,m,f,beta,breakdown,anorm);CHKERRQ(ierr);
746 dsic.upv.es!antodo 436 
      break;
961 dsic.upv.es!jroman 437 
    case EPSLANCZOS_REORTHOG_FULL:
1154 slepc 438 
      ierr = EPSBasicArnoldi(eps,PETSC_FALSE,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
1124 slepc 439 
      break;
440 
    case EPSLANCZOS_REORTHOG_DELAYED:
441 
      if (eps->orthog_ref == EPS_ORTH_REFINE_NEVER) {
1121 slepc 442 
        ierr = EPSDelayedArnoldi1(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);      
1124 slepc 443 
      } else {
1142 slepc 444 
        ierr = EPSDelayedArnoldi(eps,T,V,k,m,f,beta,breakdown);CHKERRQ(ierr);
1121 slepc 445 
      }
772 dsic.upv.es!antodo 446 
      break;
1124 slepc 447 
    default:
448 
      SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
671 dsic.upv.es!antodo 449 
  }
655 dsic.upv.es!jroman 450 
  PetscFunctionReturn(0);
451 
}
452 
 
647 dsic.upv.es!antodo 453 
#undef __FUNCT__  
454 
#define __FUNCT__ "EPSSolve_LANCZOS"
455 
PetscErrorCode EPSSolve_LANCZOS(EPS eps)
456 
{
891 dsic.upv.es!antodo 457 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
647 dsic.upv.es!antodo 458 
  PetscErrorCode ierr;
1178 slepc 459 
  int            nconv,i,j,k,n,m,*perm,restart,ncv=eps->ncv;
870 dsic.upv.es!antodo 460 
  Vec            f=eps->work[0];
1178 slepc 461 
  PetscScalar    *T=eps->T,*Y;
462 
  PetscReal      *ritz,*bnd,anorm,beta,norm;
1173 slepc 463 
  PetscTruth     breakdown;
897 dsic.upv.es!antodo 464 
  char           *conv;
647 dsic.upv.es!antodo 465 
 
466 
  PetscFunctionBegin;
666 dsic.upv.es!antodo 467 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&ritz);CHKERRQ(ierr);
1178 slepc 468 
  ierr = PetscMalloc(ncv*ncv*sizeof(PetscScalar),&Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 469 
  ierr = PetscMalloc(ncv*sizeof(PetscReal),&bnd);CHKERRQ(ierr);
687 dsic.upv.es!antodo 470 
  ierr = PetscMalloc(ncv*sizeof(int),&perm);CHKERRQ(ierr);
897 dsic.upv.es!antodo 471 
  ierr = PetscMalloc(ncv*sizeof(char),&conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 472 
 
655 dsic.upv.es!jroman 473 
  /* The first Lanczos vector is the normalized initial vector */
1057 slepc 474 
  ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);
647 dsic.upv.es!antodo 475 
 
1143 slepc 476 
  anorm = -1.0;
655 dsic.upv.es!jroman 477 
  nconv = 0;
704 dsic.upv.es!antodo 478 
 
655 dsic.upv.es!jroman 479 
  /* Restart loop */
891 dsic.upv.es!antodo 480 
  while (eps->reason == EPS_CONVERGED_ITERATING) {
1220 slepc 481 
    eps->its++;
655 dsic.upv.es!jroman 482 
    /* Compute an ncv-step Lanczos factorization */
772 dsic.upv.es!antodo 483 
    m = ncv;
891 dsic.upv.es!antodo 484 
    ierr = EPSBasicLanczos(eps,T,eps->V,nconv,&m,f,&beta,&breakdown,anorm);CHKERRQ(ierr);
647 dsic.upv.es!antodo 485 
 
1178 slepc 486 
    /* Compute eigenvalues and eigenvectors Y of the tridiagonal block */
772 dsic.upv.es!antodo 487 
    n = m - nconv;
1178 slepc 488 
    ierr = EPSDenseTridiagonal(n,T+nconv*(ncv+1),ncv,ritz,Y);CHKERRQ(ierr);
655 dsic.upv.es!jroman 489 
 
891 dsic.upv.es!antodo 490 
    /* Estimate ||A|| */
491 
    for (i=0;i<n;i++)
492 
      if (PetscAbsReal(ritz[i]) > anorm) anorm = PetscAbsReal(ritz[i]);
647 dsic.upv.es!antodo 493 
 
902 dsic.upv.es!antodo 494 
    /* Compute residual norm estimates as beta*abs(Y(m,:)) + eps*||A|| */
897 dsic.upv.es!antodo 495 
    for (i=0;i<n;i++)
1178 slepc 496 
      bnd[i] = beta*PetscAbsScalar(Y[i*n+n-1]) + PETSC_MACHINE_EPSILON*anorm;
891 dsic.upv.es!antodo 497 
 
1178 slepc 498 
    /* Sort eigenvalues according to eps->which */
499 
    if (eps->which == EPS_SMALLEST_REAL) {
500 
      /* LAPACK function has already ordered the eigenvalues and eigenvectors */
501 
      for (i=0;i<n;i++)
502 
        perm[i] = i;
503 
    } else {
870 dsic.upv.es!antodo 504 
#ifdef PETSC_USE_COMPLEX
1178 slepc 505 
      for (i=0;i<n;i++)
506 
        eps->eigr[i+nconv] = ritz[i];
507 
      ierr = EPSSortEigenvalues(n,eps->eigr+nconv,eps->eigi,eps->which,n,perm);CHKERRQ(ierr);
870 dsic.upv.es!antodo 508 
#else
1178 slepc 509 
      ierr = EPSSortEigenvalues(n,ritz,eps->eigi,eps->which,n,perm);CHKERRQ(ierr);
870 dsic.upv.es!antodo 510 
#endif
1178 slepc 511 
    }
870 dsic.upv.es!antodo 512 
 
649 dsic.upv.es!antodo 513 
    /* Look for converged eigenpairs */
870 dsic.upv.es!antodo 514 
    k = nconv;
515 
    for (i=0;i<n;i++) {
891 dsic.upv.es!antodo 516 
      eps->eigr[k] = ritz[perm[i]];
902 dsic.upv.es!antodo 517 
      eps->errest[k] = bnd[perm[i]] / PetscAbsScalar(eps->eigr[k]);    
891 dsic.upv.es!antodo 518 
      if (eps->errest[k] < eps->tol) {
902 dsic.upv.es!antodo 519 
 
1173 slepc 520 
        if (lanczos->reorthog == EPSLANCZOS_REORTHOG_LOCAL) {
902 dsic.upv.es!antodo 521 
          if (i>0 && PetscAbsScalar((eps->eigr[k]-ritz[perm[i-1]])/eps->eigr[k]) < eps->tol) {
522 
            /* Discard repeated eigenvalues */
523 
            conv[i] = 'R';
524 
            continue;
525 
          }
526 
        }
527 
 
870 dsic.upv.es!antodo 528 
        ierr = VecSet(eps->AV[k],0.0);CHKERRQ(ierr);
529 
        ierr = VecMAXPY(eps->AV[k],n,Y+perm[i]*n,eps->V+nconv);CHKERRQ(ierr);
1124 slepc 530 
 
1173 slepc 531 
        if (lanczos->reorthog == EPSLANCZOS_REORTHOG_LOCAL) {
1121 slepc 532 
          /* normalize locked vector and compute residual norm */
902 dsic.upv.es!antodo 533 
          ierr = VecNorm(eps->AV[k],NORM_2,&norm);CHKERRQ(ierr);
534 
          ierr = VecScale(eps->AV[k],1.0/norm);CHKERRQ(ierr);
1121 slepc 535 
          ierr = STApply(eps->OP,eps->AV[k],f);CHKERRQ(ierr);
536 
          ierr = VecAXPY(f,-eps->eigr[k],eps->AV[k]);CHKERRQ(ierr);
537 
          ierr = VecNorm(f,NORM_2,&norm);CHKERRQ(ierr);
538 
          eps->errest[k] = norm / PetscAbsScalar(eps->eigr[k]);
902 dsic.upv.es!antodo 539 
          if (eps->errest[k] >= eps->tol) {
540 
            conv[i] = 'S';
541 
            continue;
897 dsic.upv.es!antodo 542 
          }
895 dsic.upv.es!antodo 543 
        }
902 dsic.upv.es!antodo 544 
 
545 
        conv[i] = 'C';
546 
        k++;
547 
      } else conv[i] = 'N';
687 dsic.upv.es!antodo 548 
    }
647 dsic.upv.es!antodo 549 
 
950 dsic.upv.es!jroman 550 
    /* Look for non-converged eigenpairs */
891 dsic.upv.es!antodo 551 
    j = k;
552 
    restart = -1;
553 
    for (i=0;i<n;i++) {
897 dsic.upv.es!antodo 554 
      if (conv[i] != 'C') {
555 
        if (restart == -1 && conv[i] == 'N') restart = i;
891 dsic.upv.es!antodo 556 
        eps->eigr[j] = ritz[perm[i]];
557 
        eps->errest[j] = bnd[perm[i]] / ritz[perm[i]];
558 
        j++;
559 
      }
560 
    }
902 dsic.upv.es!antodo 561 
 
562 
    if (breakdown) {
563 
      restart = -1;
1121 slepc 564 
      PetscInfo2(eps,"Breakdown in Lanczos method (it=%i norm=%g)\n",eps->its,beta);
902 dsic.upv.es!antodo 565 
    }
891 dsic.upv.es!antodo 566 
 
897 dsic.upv.es!antodo 567 
    if (k<eps->nev) {
568 
      if (restart != -1) {
950 dsic.upv.es!jroman 569 
        /* Use first non-converged vector for restarting */
897 dsic.upv.es!antodo 570 
        ierr = VecSet(eps->AV[k],0.0);CHKERRQ(ierr);
571 
        ierr = VecMAXPY(eps->AV[k],n,Y+perm[restart]*n,eps->V+nconv);CHKERRQ(ierr);
572 
        ierr = VecCopy(eps->AV[k],eps->V[k]);CHKERRQ(ierr);
573 
      }
870 dsic.upv.es!antodo 574 
    }
575 
 
950 dsic.upv.es!jroman 576 
    /* Copy converged vectors to V */
870 dsic.upv.es!antodo 577 
    for (i=nconv;i<k;i++) {
578 
      ierr = VecCopy(eps->AV[i],eps->V[i]);CHKERRQ(ierr);
579 
    }
580 
 
897 dsic.upv.es!antodo 581 
    if (k<eps->nev) {
1057 slepc 582 
      if (restart == -1) {
583 
        /* Use random vector for restarting */
584 
        PetscInfo(eps,"Using random vector for restart\n");
585 
        ierr = EPSGetStartVector(eps,k,eps->V[k],&breakdown);CHKERRQ(ierr);
1173 slepc 586 
      } else if (lanczos->reorthog == EPSLANCZOS_REORTHOG_LOCAL) {
950 dsic.upv.es!jroman 587 
        /* Reorthonormalize restart vector */
1163 slepc 588 
        ierr = EPSOrthogonalize(eps,eps->nds+k,PETSC_NULL,eps->DSV,eps->V[k],PETSC_NULL,&norm,&breakdown);CHKERRQ(ierr);
1057 slepc 589 
        ierr = VecScale(eps->V[k],1.0/norm);CHKERRQ(ierr);
1145 slepc 590 
      } else breakdown = PETSC_FALSE;
1057 slepc 591 
      if (breakdown) {
592 
        eps->reason = EPS_DIVERGED_BREAKDOWN;
593 
        PetscInfo(eps,"Unable to generate more start vectors\n");
594 
      }
891 dsic.upv.es!antodo 595 
    }
596 
 
870 dsic.upv.es!antodo 597 
    nconv = k;
891 dsic.upv.es!antodo 598 
    EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nconv+n);
599 
    if (eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
600 
    if (nconv >= eps->nev) eps->reason = EPS_CONVERGED_TOL;
647 dsic.upv.es!antodo 601 
  }
602 
 
655 dsic.upv.es!jroman 603 
  eps->nconv = nconv;
647 dsic.upv.es!antodo 604 
 
666 dsic.upv.es!antodo 605 
  ierr = PetscFree(ritz);CHKERRQ(ierr);
647 dsic.upv.es!antodo 606 
  ierr = PetscFree(Y);CHKERRQ(ierr);
895 dsic.upv.es!antodo 607 
  ierr = PetscFree(bnd);CHKERRQ(ierr);
687 dsic.upv.es!antodo 608 
  ierr = PetscFree(perm);CHKERRQ(ierr);
895 dsic.upv.es!antodo 609 
  ierr = PetscFree(conv);CHKERRQ(ierr);
647 dsic.upv.es!antodo 610 
  PetscFunctionReturn(0);
611 
}
612 
 
1173 slepc 613 
static const char *lanczoslist[6] = { "local", "full", "selective", "periodic", "partial" , "delayed" };
1068 slepc 614 
 
671 dsic.upv.es!antodo 615 
#undef __FUNCT__  
616 
#define __FUNCT__ "EPSSetFromOptions_LANCZOS"
617 
PetscErrorCode EPSSetFromOptions_LANCZOS(EPS eps)
618 
{
619 
  PetscErrorCode ierr;
620 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
621 
  PetscTruth     flg;
982 slepc 622 
  PetscInt       i;
671 dsic.upv.es!antodo 623 
 
624 
  PetscFunctionBegin;
625 
  ierr = PetscOptionsHead("LANCZOS options");CHKERRQ(ierr);
1173 slepc 626 
  ierr = PetscOptionsEList("-eps_lanczos_reorthog","Lanczos reorthogonalization","EPSLanczosSetReorthog",lanczoslist,6,lanczoslist[lanczos->reorthog],&i,&flg);CHKERRQ(ierr);
1002 slepc 627 
  if (flg) lanczos->reorthog = (EPSLanczosReorthogType)i;
671 dsic.upv.es!antodo 628 
  ierr = PetscOptionsTail();CHKERRQ(ierr);
629 
  PetscFunctionReturn(0);
630 
}
631 
 
647 dsic.upv.es!antodo 632 
EXTERN_C_BEGIN
633 
#undef __FUNCT__  
906 dsic.upv.es!antodo 634 
#define __FUNCT__ "EPSLanczosSetReorthog_LANCZOS"
635 
PetscErrorCode EPSLanczosSetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 636 
{
637 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
638 
 
639 
  PetscFunctionBegin;
640 
  switch (reorthog) {
1068 slepc 641 
    case EPSLANCZOS_REORTHOG_LOCAL:
961 dsic.upv.es!jroman 642 
    case EPSLANCZOS_REORTHOG_FULL:
1124 slepc 643 
    case EPSLANCZOS_REORTHOG_DELAYED:
961 dsic.upv.es!jroman 644 
    case EPSLANCZOS_REORTHOG_SELECTIVE:
645 
    case EPSLANCZOS_REORTHOG_PERIODIC:
646 
    case EPSLANCZOS_REORTHOG_PARTIAL:
671 dsic.upv.es!antodo 647 
      lanczos->reorthog = reorthog;
648 
      break;
649 
    default:
650 
      SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid reorthogonalization type");
651 
  }
652 
  PetscFunctionReturn(0);
653 
}
654 
EXTERN_C_END
655 
 
656 
#undef __FUNCT__  
906 dsic.upv.es!antodo 657 
#define __FUNCT__ "EPSLanczosSetReorthog"
671 dsic.upv.es!antodo 658 
/*@
906 dsic.upv.es!antodo 659 
   EPSLanczosSetReorthog - Sets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 660 
   iteration.
661 
 
662 
   Collective on EPS
663 
 
664 
   Input Parameters:
665 
+  eps - the eigenproblem solver context
666 
-  reorthog - the type of reorthogonalization
667 
 
668 
   Options Database Key:
1068 slepc 669 
.  -eps_lanczos_orthog - Sets the reorthogonalization type (either 'local', 'selective',
950 dsic.upv.es!jroman 670 
                         'periodic', 'partial' or 'full')
671 dsic.upv.es!antodo 671 
 
672 
   Level: advanced
673 
 
906 dsic.upv.es!antodo 674 
.seealso: EPSLanczosGetReorthog()
671 dsic.upv.es!antodo 675 
@*/
906 dsic.upv.es!antodo 676 
PetscErrorCode EPSLanczosSetReorthog(EPS eps,EPSLanczosReorthogType reorthog)
671 dsic.upv.es!antodo 677 
{
906 dsic.upv.es!antodo 678 
  PetscErrorCode ierr, (*f)(EPS,EPSLanczosReorthogType);
671 dsic.upv.es!antodo 679 
 
680 
  PetscFunctionBegin;
681 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
906 dsic.upv.es!antodo 682 
  ierr = PetscObjectQueryFunction((PetscObject)eps,"EPSLanczosSetReorthog_C",(void (**)())&f);CHKERRQ(ierr);
671 dsic.upv.es!antodo 683 
  if (f) {
684 
    ierr = (*f)(eps,reorthog);CHKERRQ(ierr);
685 
  }
686 
  PetscFunctionReturn(0);
687 
}
688 
 
689 
EXTERN_C_BEGIN
690 
#undef __FUNCT__  
906 dsic.upv.es!antodo 691 
#define __FUNCT__ "EPSLanczosGetReorthog_LANCZOS"
692 
PetscErrorCode EPSLanczosGetReorthog_LANCZOS(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 693 
{
694 
  EPS_LANCZOS *lanczos = (EPS_LANCZOS *)eps->data;
695 
  PetscFunctionBegin;
696 
  *reorthog = lanczos->reorthog;
697 
  PetscFunctionReturn(0);
698 
}
699 
EXTERN_C_END
700 
 
701 
#undef __FUNCT__  
906 dsic.upv.es!antodo 702 
#define __FUNCT__ "EPSLanczosGetReorthog"
707 dsic.upv.es!antodo 703 
/*@C
906 dsic.upv.es!antodo 704 
   EPSLanczosGetReorthog - Gets the type of reorthogonalization used during the Lanczos
671 dsic.upv.es!antodo 705 
   iteration.
706 
 
707 
   Collective on EPS
708 
 
709 
   Input Parameter:
710 
.  eps - the eigenproblem solver context
711 
 
712 
   Input Parameter:
713 
.  reorthog - the type of reorthogonalization
714 
 
715 
   Level: advanced
716 
 
906 dsic.upv.es!antodo 717 
.seealso: EPSLanczosSetReorthog()
671 dsic.upv.es!antodo 718 
@*/
906 dsic.upv.es!antodo 719 
PetscErrorCode EPSLanczosGetReorthog(EPS eps,EPSLanczosReorthogType *reorthog)
671 dsic.upv.es!antodo 720 
{
906 dsic.upv.es!antodo 721 
  PetscErrorCode ierr, (*f)(EPS,EPSLanczosReorthogType*);
671 dsic.upv.es!antodo 722 
 
723 
  PetscFunctionBegin;
724 
  PetscValidHeaderSpecific(eps,EPS_COOKIE,1);
906 dsic.upv.es!antodo 725 
  ierr = PetscObjectQueryFunction((PetscObject)eps,"EPSLanczosGetReorthog_C",(void (**)())&f);CHKERRQ(ierr);
671 dsic.upv.es!antodo 726 
  if (f) {
727 
    ierr = (*f)(eps,reorthog);CHKERRQ(ierr);
728 
  }
729 
  PetscFunctionReturn(0);
730 
}
731 
 
732 
#undef __FUNCT__  
733 
#define __FUNCT__ "EPSView_LANCZOS"
734 
PetscErrorCode EPSView_LANCZOS(EPS eps,PetscViewer viewer)
735 
{
736 
  PetscErrorCode ierr;
737 
  EPS_LANCZOS    *lanczos = (EPS_LANCZOS *)eps->data;
738 
  PetscTruth     isascii;
739 
 
740 
  PetscFunctionBegin;
741 
  ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr);
742 
  if (!isascii) {
743 
    SETERRQ1(1,"Viewer type %s not supported for EPSLANCZOS",((PetscObject)viewer)->type_name);
744 
  }  
1068 slepc 745 
  ierr = PetscViewerASCIIPrintf(viewer,"reorthogonalization: %s\n",lanczoslist[lanczos->reorthog]);CHKERRQ(ierr);
671 dsic.upv.es!antodo 746 
  PetscFunctionReturn(0);
747 
}
748 
 
961 dsic.upv.es!jroman 749 
/*
750 
EXTERN PetscErrorCode EPSSolve_TS_LANCZOS(EPS);
751 
*/
906 dsic.upv.es!antodo 752 
 
671 dsic.upv.es!antodo 753 
EXTERN_C_BEGIN
754 
#undef __FUNCT__  
647 dsic.upv.es!antodo 755 
#define __FUNCT__ "EPSCreate_LANCZOS"
756 
PetscErrorCode EPSCreate_LANCZOS(EPS eps)
757 
{
671 dsic.upv.es!antodo 758 
  PetscErrorCode ierr;
759 
  EPS_LANCZOS    *lanczos;
760 
 
647 dsic.upv.es!antodo 761 
  PetscFunctionBegin;
671 dsic.upv.es!antodo 762 
  ierr = PetscNew(EPS_LANCZOS,&lanczos);CHKERRQ(ierr);
763 
  PetscLogObjectMemory(eps,sizeof(EPS_LANCZOS));
764 
  eps->data                      = (void *) lanczos;
647 dsic.upv.es!antodo 765 
  eps->ops->solve                = EPSSolve_LANCZOS;
950 dsic.upv.es!jroman 766 
/*  eps->ops->solvets              = EPSSolve_TS_LANCZOS;*/
647 dsic.upv.es!antodo 767 
  eps->ops->setup                = EPSSetUp_LANCZOS;
671 dsic.upv.es!antodo 768 
  eps->ops->setfromoptions       = EPSSetFromOptions_LANCZOS;
647 dsic.upv.es!antodo 769 
  eps->ops->destroy              = EPSDestroy_Default;
671 dsic.upv.es!antodo 770 
  eps->ops->view                 = EPSView_LANCZOS;
647 dsic.upv.es!antodo 771 
  eps->ops->backtransform        = EPSBackTransform_Default;
950 dsic.upv.es!jroman 772 
  /*if (eps->solverclass==EPS_TWO_SIDE)
885 dsic.upv.es!jroman 773 
       eps->ops->computevectors       = EPSComputeVectors_Schur;
950 dsic.upv.es!jroman 774 
  else*/ eps->ops->computevectors       = EPSComputeVectors_Default;
1068 slepc 775 
  lanczos->reorthog              = EPSLANCZOS_REORTHOG_LOCAL;
906 dsic.upv.es!antodo 776 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosSetReorthog_C","EPSLanczosSetReorthog_LANCZOS",EPSLanczosSetReorthog_LANCZOS);CHKERRQ(ierr);
777 
  ierr = PetscObjectComposeFunctionDynamic((PetscObject)eps,"EPSLanczosGetReorthog_C","EPSLanczosGetReorthog_LANCZOS",EPSLanczosGetReorthog_LANCZOS);CHKERRQ(ierr);
647 dsic.upv.es!antodo 778 
  PetscFunctionReturn(0);
779 
}
780 
EXTERN_C_END
781