Actual source code: epsdefault.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    This file contains some simple default routines for common operations
     12: */

     14: #include <slepc/private/epsimpl.h>
     15: #include <slepcvec.h>

     17: PetscErrorCode EPSBackTransform_Default(EPS eps)
     18: {
     19:   PetscFunctionBegin;
     20:   PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
     21:   PetscFunctionReturn(PETSC_SUCCESS);
     22: }

     24: /*
     25:   EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
     26:   using purification for generalized eigenproblems.
     27:  */
     28: PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
     29: {
     30:   PetscBool      iscayley,indef;
     31:   Mat            B,C;

     33:   PetscFunctionBegin;
     34:   if (eps->purify) {
     35:     PetscCall(EPS_Purify(eps,eps->nconv));
     36:     PetscCall(BVNormalize(eps->V,NULL));
     37:   } else {
     38:     /* In the case of Cayley transform, eigenvectors need to be B-normalized */
     39:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
     40:     if (iscayley && eps->isgeneralized) {
     41:       PetscCall(STGetMatrix(eps->st,1,&B));
     42:       PetscCall(BVGetMatrix(eps->V,&C,&indef));
     43:       PetscCheck(!indef,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
     44:       PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
     45:       PetscCall(BVNormalize(eps->V,NULL));
     46:       PetscCall(BVSetMatrix(eps->V,C,PETSC_FALSE));  /* restore original matrix */
     47:     }
     48:   }
     49:   PetscFunctionReturn(PETSC_SUCCESS);
     50: }

     52: /*
     53:   EPSComputeVectors_Indefinite - similar to the Schur version but
     54:   for indefinite problems
     55:  */
     56: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
     57: {
     58:   PetscInt       n;
     59:   Mat            X;

     61:   PetscFunctionBegin;
     62:   PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
     63:   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
     64:   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
     65:   PetscCall(BVMultInPlace(eps->V,X,0,n));
     66:   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));

     68:   /* purification */
     69:   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));

     71:   /* normalization */
     72:   PetscCall(BVNormalize(eps->V,eps->eigi));
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: /*
     77:   EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
     78:  */
     79: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
     80: {
     81:   PetscInt       i;
     82:   Vec            w,y;

     84:   PetscFunctionBegin;
     85:   if (!eps->twosided || !eps->isgeneralized) PetscFunctionReturn(PETSC_SUCCESS);
     86:   PetscCall(EPSSetWorkVecs(eps,1));
     87:   w = eps->work[0];
     88:   for (i=0;i<eps->nconv;i++) {
     89:     PetscCall(BVCopyVec(eps->W,i,w));
     90:     PetscCall(BVGetColumn(eps->W,i,&y));
     91:     PetscCall(STMatSolveHermitianTranspose(eps->st,w,y));
     92:     PetscCall(BVRestoreColumn(eps->W,i,&y));
     93:   }
     94:   PetscFunctionReturn(PETSC_SUCCESS);
     95: }

     97: /*
     98:   EPSComputeVectors_Schur - Compute eigenvectors from the vectors
     99:   provided by the eigensolver. This version is intended for solvers
    100:   that provide Schur vectors. Given the partial Schur decomposition
    101:   OP*V=V*T, the following steps are performed:
    102:       1) compute eigenvectors of T: T*Z=Z*D
    103:       2) compute eigenvectors of OP: X=V*Z
    104:  */
    105: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
    106: {
    107:   PetscInt       i;
    108:   Mat            Z;
    109:   Vec            z;

    111:   PetscFunctionBegin;
    112:   if (eps->ishermitian) {
    113:     if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
    114:     else PetscCall(EPSComputeVectors_Hermitian(eps));
    115:     PetscFunctionReturn(PETSC_SUCCESS);
    116:   }

    118:   /* right eigenvectors */
    119:   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));

    121:   /* V = V * Z */
    122:   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
    123:   PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
    124:   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));

    126:   /* Purify eigenvectors */
    127:   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));

    129:   /* Fix eigenvectors if balancing was used */
    130:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
    131:     for (i=0;i<eps->nconv;i++) {
    132:       PetscCall(BVGetColumn(eps->V,i,&z));
    133:       PetscCall(VecPointwiseDivide(z,z,eps->D));
    134:       PetscCall(BVRestoreColumn(eps->V,i,&z));
    135:     }
    136:   }

    138:   /* normalize eigenvectors (when using purification or balancing) */
    139:   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));

    141:   /* left eigenvectors */
    142:   if (eps->twosided) {
    143:     PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
    144:     /* W = W * Z */
    145:     PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
    146:     PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
    147:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
    148:     /* Fix left eigenvectors if balancing was used */
    149:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
    150:       for (i=0;i<eps->nconv;i++) {
    151:         PetscCall(BVGetColumn(eps->W,i,&z));
    152:         PetscCall(VecPointwiseMult(z,z,eps->D));
    153:         PetscCall(BVRestoreColumn(eps->W,i,&z));
    154:       }
    155:     }
    156:     PetscCall(EPSComputeVectors_Twosided(eps));
    157:     /* normalize */
    158:     PetscCall(BVNormalize(eps->W,eps->eigi));
    159: #if !defined(PETSC_USE_COMPLEX)
    160:     for (i=0;i<eps->nconv-1;i++) {
    161:       if (eps->eigi[i] != 0.0) {
    162:         if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
    163:         i++;
    164:       }
    165:     }
    166: #endif
    167:   }
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: /*@
    172:    EPSSetWorkVecs - Sets a number of work vectors into an `EPS` object.

    174:    Collective

    176:    Input Parameters:
    177: +  eps - the linear eigensolver context
    178: -  nw  - number of work vectors to allocate

    180:    Developer Note:
    181:    This is `SLEPC_EXTERN` because it may be required by user plugin `EPS`
    182:    implementations.

    184:    Level: developer

    186: .seealso: [](ch:eps), `EPSSetUp()`
    187: @*/
    188: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
    189: {
    190:   Vec            t;

    192:   PetscFunctionBegin;
    195:   PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
    196:   if (eps->nwork < nw) {
    197:     PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
    198:     eps->nwork = nw;
    199:     PetscCall(BVGetColumn(eps->V,0,&t));
    200:     PetscCall(VecDuplicateVecs(t,nw,&eps->work));
    201:     PetscCall(BVRestoreColumn(eps->V,0,&t));
    202:   }
    203:   PetscFunctionReturn(PETSC_SUCCESS);
    204: }

    206: /*
    207:   EPSSetWhichEigenpairs_Default - Sets the default value for which,
    208:   depending on the ST.
    209:  */
    210: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
    211: {
    212:   PetscBool      target;

    214:   PetscFunctionBegin;
    215:   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
    216:   if (target) eps->which = EPS_TARGET_MAGNITUDE;
    217:   else eps->which = EPS_LARGEST_MAGNITUDE;
    218:   PetscFunctionReturn(PETSC_SUCCESS);
    219: }

    221: /*
    222:   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
    223: */
    224: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
    225: {
    226:   PetscReal w;

    228:   PetscFunctionBegin;
    229:   w = SlepcAbsEigenvalue(eigr,eigi);
    230:   *errest = (w!=0.0)? res/w: PETSC_MAX_REAL;
    231:   PetscFunctionReturn(PETSC_SUCCESS);
    232: }

    234: /*
    235:   EPSConvergedAbsolute - Checks convergence absolutely.
    236: */
    237: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
    238: {
    239:   PetscFunctionBegin;
    240:   *errest = res;
    241:   PetscFunctionReturn(PETSC_SUCCESS);
    242: }

    244: /*
    245:   EPSConvergedNorm - Checks convergence relative to the eigenvalue and
    246:   the matrix norms.
    247: */
    248: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
    249: {
    250:   PetscReal w;

    252:   PetscFunctionBegin;
    253:   w = SlepcAbsEigenvalue(eigr,eigi);
    254:   *errest = res / (eps->nrma + w*eps->nrmb);
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: /*@C
    259:    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
    260:    iteration must be stopped.

    262:    Collective

    264:    Input Parameters:
    265: +  eps    - the linear eigensolver context
    266: .  its    - current number of iterations
    267: .  max_it - maximum number of iterations
    268: .  nconv  - number of currently converged eigenpairs
    269: .  nev    - number of requested eigenpairs
    270: -  ctx    - context (not used here)

    272:    Output Parameter:
    273: .  reason - result of the stopping test

    275:    Notes:
    276:    `EPSStoppingBasic()` will stop if all requested eigenvalues are converged, or if
    277:    the maximum number of iterations has been reached.

    279:    This is the default stopping test.
    280:    Use `EPSSetStoppingTest()` to provide your own test instead of using this one.

    282:    Level: advanced

    284: .seealso: [](ch:eps), `EPSSetStoppingTest()`, `EPSStoppingThreshold()`, `EPSConvergedReason`, `EPSGetConvergedReason()`
    285: @*/
    286: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
    287: {
    288:   PetscFunctionBegin;
    289:   *reason = EPS_CONVERGED_ITERATING;
    290:   if (nconv >= nev) {
    291:     PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
    292:     *reason = EPS_CONVERGED_TOL;
    293:   } else if (its >= max_it) {
    294:     *reason = EPS_DIVERGED_ITS;
    295:     PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
    296:   }
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: /*@C
    301:    EPSStoppingThreshold - Routine to determine whether the outer eigenvalue solver
    302:    iteration must be stopped, according to some threshold for the computed values.

    304:    Collective

    306:    Input Parameters:
    307: +  eps    - the linear eigensolver context
    308: .  its    - current number of iterations
    309: .  max_it - maximum number of iterations
    310: .  nconv  - number of currently converged eigenpairs
    311: .  nev    - number of requested eigenpairs
    312: -  ctx    - context containing additional data (`EPSStoppingCtx`)

    314:    Output Parameter:
    315: .  reason - result of the stopping test

    317:    Notes:
    318:    `EPSStoppingThreshold()` will stop when one of the computed eigenvalues is not
    319:    above/below the threshold given at `EPSSetThreshold()`. If a number of wanted
    320:    eigenvalues has been specified via `EPSSetDimensions()` then it is also taken into
    321:    account, and the solver will stop when one of the two conditions (threshold or
    322:    number of converged values) is met.

    324:    Use `EPSSetStoppingTest()` to provide your own test instead of using this one.

    326:    Level: advanced

    328: .seealso: [](ch:eps), `EPSSetStoppingTest()`, `EPSStoppingBasic()`, `EPSSetThreshold()`, `EPSSetDimensions()`, `EPSConvergedReason`, `EPSGetConvergedReason()`
    329: @*/
    330: PetscErrorCode EPSStoppingThreshold(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
    331: {
    332:   PetscReal thres,firstev,lastev,firstnc,errest;
    333:   PetscBool magnit,rel,isfilter;
    334:   PetscInt  napprox;
    335:   EPSWhich  which;

    337:   PetscFunctionBegin;
    338:   *reason = EPS_CONVERGED_ITERATING;
    339:   firstev = ((EPSStoppingCtx)ctx)->firstev;
    340:   lastev  = ((EPSStoppingCtx)ctx)->lastev;
    341:   firstnc = ((EPSStoppingCtx)ctx)->firstnc;
    342:   errest  = ((EPSStoppingCtx)ctx)->errest;
    343:   thres   = ((EPSStoppingCtx)ctx)->thres;
    344:   rel     = ((EPSStoppingCtx)ctx)->threlative;
    345:   napprox = ((EPSStoppingCtx)ctx)->napprox;
    346:   which   = ((EPSStoppingCtx)ctx)->which;

    348:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
    349:   if (isfilter) {
    350:     /*
    351:       ((EPSStoppingCtx)ctx)->threlative   should be PETSC_FALSE
    352:       ((EPSStoppingCtx)ctx)->which        should be EPS_LARGEST_MAGNITUDE
    353:     */
    354:     if (nconv && (nconv==napprox || firstnc+errest<thres)) {
    355:       if (its==((EPSStoppingCtx)ctx)->its+1) {
    356:         if (nconv<napprox) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g (plus error estimate %g) is below the threshold %g\n",(double)firstnc,(double)errest,(double)thres));
    357:         else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: all available eigenvalue approximations have converged\n"));
    358:         *reason = EPS_CONVERGED_TOL;
    359:       } else ((EPSStoppingCtx)ctx)->its = its;  /* wait until next iteration */
    360:     } else if (nev && nconv >= nev) {
    361:       PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
    362:       *reason = EPS_CONVERGED_TOL;
    363:     } else if (its >= max_it) {
    364:       *reason = EPS_DIVERGED_ITS;
    365:       PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
    366:     }
    367:   } else {
    368:     magnit  = (which==EPS_SMALLEST_MAGNITUDE || which==EPS_LARGEST_MAGNITUDE || which==EPS_TARGET_MAGNITUDE)? PETSC_TRUE: PETSC_FALSE;
    369:     if (nconv && magnit && which==EPS_TARGET_MAGNITUDE && ((rel && ((thres>1.0 && lastev>thres*firstev) || (thres<1.0 && lastev<thres*firstev))) || (!rel && lastev>thres))) {
    370:       if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
    371:       else if (thres>1.0) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is above the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
    372:       else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
    373:       *reason = EPS_CONVERGED_TOL;
    374:     } else if (nconv && magnit && ((which==EPS_LARGEST_MAGNITUDE && ((rel && lastev<thres*firstev) || (!rel && lastev<thres))) || (which==EPS_SMALLEST_MAGNITUDE && lastev>thres))) {
    375:       if (which==EPS_SMALLEST_MAGNITUDE) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
    376:       else if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is below the threshold %g\n",(double)lastev,(double)thres));
    377:       else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
    378:       *reason = EPS_CONVERGED_TOL;
    379:     } else if (nconv && !magnit && ((which==EPS_LARGEST_REAL && lastev<thres) || (which==EPS_SMALLEST_REAL && lastev>thres))) {
    380:       if (which==EPS_LARGEST_REAL) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is below the threshold %g\n",(double)lastev,(double)thres));
    381:       else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is above the threshold %g\n",(double)lastev,(double)thres));
    382:       *reason = EPS_CONVERGED_TOL;
    383:     } else if (nev && nconv >= nev) {
    384:       PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
    385:       *reason = EPS_CONVERGED_TOL;
    386:     } else if (its >= max_it) {
    387:       *reason = EPS_DIVERGED_ITS;
    388:       PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
    389:     }
    390:   }
    391:   PetscFunctionReturn(PETSC_SUCCESS);
    392: }

    394: /*
    395:   EPSComputeRitzVector - Computes the current Ritz vector.

    397:   Simple case (complex scalars or real scalars with Zi=NULL):
    398:     x = V*Zr  (V is a basis of nv vectors, Zr has length nv)

    400:   Split case:
    401:     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
    402: */
    403: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
    404: {
    405:   PetscInt       l,k;
    406:   PetscReal      norm;
    407: #if !defined(PETSC_USE_COMPLEX)
    408:   Vec            z;
    409: #endif

    411:   PetscFunctionBegin;
    412:   /* compute eigenvector */
    413:   PetscCall(BVGetActiveColumns(V,&l,&k));
    414:   PetscCall(BVSetActiveColumns(V,0,k));
    415:   PetscCall(BVMultVec(V,1.0,0.0,x,Zr));

    417:   /* purify eigenvector if necessary */
    418:   if (eps->purify) {
    419:     PetscCall(STApply(eps->st,x,y));
    420:     if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
    421:     else PetscCall(VecNorm(y,NORM_2,&norm));
    422:     PetscCall(VecScale(y,1.0/norm));
    423:     PetscCall(VecCopy(y,x));
    424:   }
    425:   /* fix eigenvector if balancing is used */
    426:   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
    427: #if !defined(PETSC_USE_COMPLEX)
    428:   /* compute imaginary part of eigenvector */
    429:   if (Zi) {
    430:     PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
    431:     if (eps->ispositive) {
    432:       PetscCall(BVCreateVec(V,&z));
    433:       PetscCall(STApply(eps->st,y,z));
    434:       PetscCall(VecNorm(z,NORM_2,&norm));
    435:       PetscCall(VecScale(z,1.0/norm));
    436:       PetscCall(VecCopy(z,y));
    437:       PetscCall(VecDestroy(&z));
    438:     }
    439:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
    440:   } else
    441: #endif
    442:     PetscCall(VecSet(y,0.0));

    444:   /* normalize eigenvectors (when using balancing) */
    445:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
    446: #if !defined(PETSC_USE_COMPLEX)
    447:     if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
    448:     else
    449: #endif
    450:     PetscCall(VecNormalize(x,NULL));
    451:   }
    452:   PetscCall(BVSetActiveColumns(V,l,k));
    453:   PetscFunctionReturn(PETSC_SUCCESS);
    454: }

    456: /*
    457:   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
    458:   diagonal matrix to be applied for balancing in non-Hermitian problems.
    459: */
    460: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
    461: {
    462:   Vec               z,p,r;
    463:   PetscInt          i,j;
    464:   PetscReal         norma;
    465:   PetscScalar       *pz,*pD;
    466:   const PetscScalar *pr,*pp;
    467:   PetscRandom       rand;

    469:   PetscFunctionBegin;
    470:   PetscCall(EPSSetWorkVecs(eps,3));
    471:   PetscCall(BVGetRandomContext(eps->V,&rand));
    472:   r = eps->work[0];
    473:   p = eps->work[1];
    474:   z = eps->work[2];
    475:   PetscCall(VecSet(eps->D,1.0));

    477:   for (j=0;j<eps->balance_its;j++) {

    479:     /* Build a random vector of +-1's */
    480:     PetscCall(VecSetRandom(z,rand));
    481:     PetscCall(VecGetArray(z,&pz));
    482:     for (i=0;i<eps->nloc;i++) {
    483:       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
    484:       else pz[i]=1.0;
    485:     }
    486:     PetscCall(VecRestoreArray(z,&pz));

    488:     /* Compute p=DA(D\z) */
    489:     PetscCall(VecPointwiseDivide(r,z,eps->D));
    490:     PetscCall(STApply(eps->st,r,p));
    491:     PetscCall(VecPointwiseMult(p,p,eps->D));
    492:     if (eps->balance == EPS_BALANCE_TWOSIDE) {
    493:       if (j==0) {
    494:         /* Estimate the matrix inf-norm */
    495:         PetscCall(VecAbs(p));
    496:         PetscCall(VecMax(p,NULL,&norma));
    497:       }
    498:       /* Compute r=D\(A'Dz) */
    499:       PetscCall(VecPointwiseMult(z,z,eps->D));
    500:       PetscCall(STApplyHermitianTranspose(eps->st,z,r));
    501:       PetscCall(VecPointwiseDivide(r,r,eps->D));
    502:     }

    504:     /* Adjust values of D */
    505:     PetscCall(VecGetArrayRead(r,&pr));
    506:     PetscCall(VecGetArrayRead(p,&pp));
    507:     PetscCall(VecGetArray(eps->D,&pD));
    508:     for (i=0;i<eps->nloc;i++) {
    509:       if (eps->balance == EPS_BALANCE_TWOSIDE) {
    510:         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
    511:           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
    512:       } else {
    513:         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
    514:       }
    515:     }
    516:     PetscCall(VecRestoreArrayRead(r,&pr));
    517:     PetscCall(VecRestoreArrayRead(p,&pp));
    518:     PetscCall(VecRestoreArray(eps->D,&pD));
    519:   }
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }