Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    Simple default routines for common PEP operations
      12             : */
      13             : 
      14             : #include <slepc/private/pepimpl.h>     /*I "slepcpep.h" I*/
      15             : 
      16             : /*@
      17             :    PEPSetWorkVecs - Sets a number of work vectors into a PEP object.
      18             : 
      19             :    Collective
      20             : 
      21             :    Input Parameters:
      22             : +  pep - polynomial eigensolver context
      23             : -  nw  - number of work vectors to allocate
      24             : 
      25             :    Developer Notes:
      26             :    This is SLEPC_EXTERN because it may be required by user plugin PEP
      27             :    implementations.
      28             : 
      29             :    Level: developer
      30             : 
      31             : .seealso: PEPSetUp()
      32             : @*/
      33        1217 : PetscErrorCode PEPSetWorkVecs(PEP pep,PetscInt nw)
      34             : {
      35        1217 :   Vec            t;
      36             : 
      37        1217 :   PetscFunctionBegin;
      38        1217 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
      39        4868 :   PetscValidLogicalCollectiveInt(pep,nw,2);
      40        1217 :   PetscCheck(nw>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
      41        1217 :   if (pep->nwork < nw) {
      42         206 :     PetscCall(VecDestroyVecs(pep->nwork,&pep->work));
      43         206 :     pep->nwork = nw;
      44         206 :     PetscCall(BVGetColumn(pep->V,0,&t));
      45         206 :     PetscCall(VecDuplicateVecs(t,nw,&pep->work));
      46         206 :     PetscCall(BVRestoreColumn(pep->V,0,&t));
      47             :   }
      48        1217 :   PetscFunctionReturn(PETSC_SUCCESS);
      49             : }
      50             : 
      51             : /*
      52             :   PEPConvergedRelative - Checks convergence relative to the eigenvalue.
      53             : */
      54        2820 : PetscErrorCode PEPConvergedRelative(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
      55             : {
      56        2820 :   PetscReal w;
      57             : 
      58        2820 :   PetscFunctionBegin;
      59        2820 :   w = SlepcAbsEigenvalue(eigr,eigi);
      60        2820 :   *errest = res/w;
      61        2820 :   PetscFunctionReturn(PETSC_SUCCESS);
      62             : }
      63             : 
      64             : /*
      65             :   PEPConvergedNorm - Checks convergence relative to the matrix norms.
      66             : */
      67          11 : PetscErrorCode PEPConvergedNorm(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
      68             : {
      69          11 :   PetscReal      w=0.0,t;
      70          11 :   PetscInt       j;
      71          11 :   PetscBool      flg;
      72             : 
      73          11 :   PetscFunctionBegin;
      74             :   /* initialization of matrix norms */
      75          11 :   if (!pep->nrma[pep->nmat-1]) {
      76           4 :     for (j=0;j<pep->nmat;j++) {
      77           3 :       PetscCall(MatHasOperation(pep->A[j],MATOP_NORM,&flg));
      78           3 :       PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
      79           3 :       PetscCall(MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]));
      80             :     }
      81             :   }
      82          11 :   t = SlepcAbsEigenvalue(eigr,eigi);
      83          44 :   for (j=pep->nmat-1;j>=0;j--) {
      84          33 :     w = w*t+pep->nrma[j];
      85             :   }
      86          11 :   *errest = res/w;
      87          11 :   PetscFunctionReturn(PETSC_SUCCESS);
      88             : }
      89             : 
      90             : /*
      91             :   PEPSetWhichEigenpairs_Default - Sets the default value for which,
      92             :   depending on the ST.
      93             :  */
      94          67 : PetscErrorCode PEPSetWhichEigenpairs_Default(PEP pep)
      95             : {
      96          67 :   PetscBool      target;
      97             : 
      98          67 :   PetscFunctionBegin;
      99          67 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&target));
     100          67 :   if (target) pep->which = PEP_TARGET_MAGNITUDE;
     101          49 :   else pep->which = PEP_LARGEST_MAGNITUDE;
     102          67 :   PetscFunctionReturn(PETSC_SUCCESS);
     103             : }
     104             : 
     105             : /*
     106             :   PEPConvergedAbsolute - Checks convergence absolutely.
     107             : */
     108          25 : PetscErrorCode PEPConvergedAbsolute(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     109             : {
     110          25 :   PetscFunctionBegin;
     111          25 :   *errest = res;
     112          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     113             : }
     114             : 
     115             : /*@C
     116             :    PEPStoppingBasic - Default routine to determine whether the outer eigensolver
     117             :    iteration must be stopped.
     118             : 
     119             :    Collective
     120             : 
     121             :    Input Parameters:
     122             : +  pep    - eigensolver context obtained from PEPCreate()
     123             : .  its    - current number of iterations
     124             : .  max_it - maximum number of iterations
     125             : .  nconv  - number of currently converged eigenpairs
     126             : .  nev    - number of requested eigenpairs
     127             : -  ctx    - context (not used here)
     128             : 
     129             :    Output Parameter:
     130             : .  reason - result of the stopping test
     131             : 
     132             :    Notes:
     133             :    A positive value of reason indicates that the iteration has finished successfully
     134             :    (converged), and a negative value indicates an error condition (diverged). If
     135             :    the iteration needs to be continued, reason must be set to PEP_CONVERGED_ITERATING
     136             :    (zero).
     137             : 
     138             :    PEPStoppingBasic() will stop if all requested eigenvalues are converged, or if
     139             :    the maximum number of iterations has been reached.
     140             : 
     141             :    Use PEPSetStoppingTest() to provide your own test instead of using this one.
     142             : 
     143             :    Level: advanced
     144             : 
     145             : .seealso: PEPSetStoppingTest(), PEPConvergedReason, PEPGetConvergedReason()
     146             : @*/
     147        1701 : PetscErrorCode PEPStoppingBasic(PEP pep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,PEPConvergedReason *reason,void *ctx)
     148             : {
     149        1701 :   PetscFunctionBegin;
     150        1701 :   *reason = PEP_CONVERGED_ITERATING;
     151        1701 :   if (nconv >= nev) {
     152         160 :     PetscCall(PetscInfo(pep,"Polynomial eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
     153         160 :     *reason = PEP_CONVERGED_TOL;
     154        1541 :   } else if (its >= max_it) {
     155           1 :     *reason = PEP_DIVERGED_ITS;
     156           1 :     PetscCall(PetscInfo(pep,"Polynomial eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
     157             :   }
     158        1701 :   PetscFunctionReturn(PETSC_SUCCESS);
     159             : }
     160             : 
     161         115 : PetscErrorCode PEPBackTransform_Default(PEP pep)
     162             : {
     163         115 :   PetscFunctionBegin;
     164         115 :   PetscCall(STBackTransform(pep->st,pep->nconv,pep->eigr,pep->eigi));
     165         115 :   PetscFunctionReturn(PETSC_SUCCESS);
     166             : }
     167             : 
     168         166 : PetscErrorCode PEPComputeVectors_Default(PEP pep)
     169             : {
     170         166 :   PetscInt       i;
     171         166 :   Vec            v;
     172             : 
     173         166 :   PetscFunctionBegin;
     174         166 :   PetscCall(PEPExtractVectors(pep));
     175             : 
     176             :   /* Fix eigenvectors if balancing was used */
     177         166 :   if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr && (pep->refine!=PEP_REFINE_MULTIPLE)) {
     178          89 :     for (i=0;i<pep->nconv;i++) {
     179          75 :       PetscCall(BVGetColumn(pep->V,i,&v));
     180          75 :       PetscCall(VecPointwiseMult(v,v,pep->Dr));
     181          75 :       PetscCall(BVRestoreColumn(pep->V,i,&v));
     182             :     }
     183             :   }
     184             : 
     185             :   /* normalization */
     186         166 :   PetscCall(BVNormalize(pep->V,pep->eigi));
     187         166 :   PetscFunctionReturn(PETSC_SUCCESS);
     188             : }
     189             : 
     190             : /*
     191             :   PEPBuildDiagonalScaling - compute two diagonal matrices to be applied for balancing
     192             :   in polynomial eigenproblems.
     193             : */
     194          18 : PetscErrorCode PEPBuildDiagonalScaling(PEP pep)
     195             : {
     196          18 :   PetscInt       it,i,j,k,nmat,nr,e,nz,lst,lend,nc=0,*cols,emax,emin,emaxl,eminl;
     197          18 :   const PetscInt *cidx,*ridx;
     198          18 :   Mat            M,*T,A;
     199          18 :   PetscMPIInt    n;
     200          18 :   PetscBool      cont=PETSC_TRUE,flg=PETSC_FALSE;
     201          18 :   PetscScalar    *array,*Dr,*Dl,t;
     202          18 :   PetscReal      l2,d,*rsum,*aux,*csum,w=1.0;
     203          18 :   MatStructure   str;
     204          18 :   MatInfo        info;
     205             : 
     206          18 :   PetscFunctionBegin;
     207          18 :   l2 = 2*PetscLogReal(2.0);
     208          18 :   nmat = pep->nmat;
     209          18 :   PetscCall(PetscMPIIntCast(pep->n,&n));
     210          18 :   PetscCall(STGetMatStructure(pep->st,&str));
     211          18 :   PetscCall(PetscMalloc1(nmat,&T));
     212          72 :   for (k=0;k<nmat;k++) PetscCall(STGetMatrixTransformed(pep->st,k,&T[k]));
     213             :   /* Form local auxiliary matrix M */
     214          18 :   PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)T[0],&cont,MATMPIAIJ,MATSEQAIJ,""));
     215          18 :   PetscCheck(cont,PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"Only for MPIAIJ or SEQAIJ matrix types");
     216          18 :   PetscCall(PetscObjectBaseTypeCompare((PetscObject)T[0],MATMPIAIJ,&cont));
     217          18 :   if (cont) {
     218          12 :     PetscCall(MatMPIAIJGetLocalMat(T[0],MAT_INITIAL_MATRIX,&M));
     219             :     flg = PETSC_TRUE;
     220           6 :   } else PetscCall(MatDuplicate(T[0],MAT_COPY_VALUES,&M));
     221          18 :   PetscCall(MatGetInfo(M,MAT_LOCAL,&info));
     222          18 :   nz = (PetscInt)info.nz_used;
     223          18 :   PetscCall(MatSeqAIJGetArray(M,&array));
     224         800 :   for (i=0;i<nz;i++) {
     225         782 :     t = PetscAbsScalar(array[i]);
     226         782 :     array[i] = t*t;
     227             :   }
     228          18 :   PetscCall(MatSeqAIJRestoreArray(M,&array));
     229          54 :   for (k=1;k<nmat;k++) {
     230          36 :     if (flg) PetscCall(MatMPIAIJGetLocalMat(T[k],MAT_INITIAL_MATRIX,&A));
     231             :     else {
     232          12 :       if (str==SAME_NONZERO_PATTERN) PetscCall(MatCopy(T[k],A,SAME_NONZERO_PATTERN));
     233          12 :       else PetscCall(MatDuplicate(T[k],MAT_COPY_VALUES,&A));
     234             :     }
     235          36 :     PetscCall(MatGetInfo(A,MAT_LOCAL,&info));
     236          36 :     nz = (PetscInt)info.nz_used;
     237          36 :     PetscCall(MatSeqAIJGetArray(A,&array));
     238        1206 :     for (i=0;i<nz;i++) {
     239        1170 :       t = PetscAbsScalar(array[i]);
     240        1170 :       array[i] = t*t;
     241             :     }
     242          36 :     PetscCall(MatSeqAIJRestoreArray(A,&array));
     243          36 :     w *= pep->slambda*pep->slambda*pep->sfactor;
     244          36 :     PetscCall(MatAXPY(M,w,A,str));
     245          36 :     if (flg || str!=SAME_NONZERO_PATTERN || k==nmat-2) PetscCall(MatDestroy(&A));
     246             :   }
     247          18 :   PetscCall(MatGetRowIJ(M,0,PETSC_FALSE,PETSC_FALSE,&nr,&ridx,&cidx,&cont));
     248          18 :   PetscCheck(cont,PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"It is not possible to compute scaling diagonals for these PEP matrices");
     249          18 :   PetscCall(MatGetInfo(M,MAT_LOCAL,&info));
     250          18 :   nz = (PetscInt)info.nz_used;
     251          18 :   PetscCall(VecGetOwnershipRange(pep->Dl,&lst,&lend));
     252          18 :   PetscCall(PetscMalloc4(nr,&rsum,pep->n,&csum,pep->n,&aux,PetscMin(pep->n-lend+lst,nz),&cols));
     253          18 :   PetscCall(VecSet(pep->Dr,1.0));
     254          18 :   PetscCall(VecSet(pep->Dl,1.0));
     255          18 :   PetscCall(VecGetArray(pep->Dl,&Dl));
     256          18 :   PetscCall(VecGetArray(pep->Dr,&Dr));
     257          18 :   PetscCall(MatSeqAIJGetArray(M,&array));
     258          18 :   PetscCall(PetscArrayzero(aux,pep->n));
     259         846 :   for (j=0;j<nz;j++) {
     260             :     /* Search non-zero columns outsize lst-lend */
     261         828 :     if (aux[cidx[j]]==0 && (cidx[j]<lst || lend<=cidx[j])) cols[nc++] = cidx[j];
     262             :     /* Local column sums */
     263         828 :     aux[cidx[j]] += PetscAbsScalar(array[j]);
     264             :   }
     265          54 :   for (it=0;it<pep->sits && cont;it++) {
     266          36 :     emaxl = 0; eminl = 0;
     267             :     /* Column sum  */
     268          36 :     if (it>0) { /* it=0 has been already done*/
     269          18 :       PetscCall(MatSeqAIJGetArray(M,&array));
     270          18 :       PetscCall(PetscArrayzero(aux,pep->n));
     271         846 :       for (j=0;j<nz;j++) aux[cidx[j]] += PetscAbsScalar(array[j]);
     272          18 :       PetscCall(MatSeqAIJRestoreArray(M,&array));
     273             :     }
     274         144 :     PetscCall(MPIU_Allreduce(aux,csum,n,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)pep->Dr)));
     275             :     /* Update Dr */
     276         628 :     for (j=lst;j<lend;j++) {
     277         592 :       d = PetscLogReal(csum[j])/l2;
     278         592 :       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
     279         592 :       d = PetscPowReal(2.0,e);
     280         592 :       Dr[j-lst] *= d;
     281         592 :       aux[j] = d*d;
     282         592 :       emaxl = PetscMax(emaxl,e);
     283         592 :       eminl = PetscMin(eminl,e);
     284             :     }
     285          64 :     for (j=0;j<nc;j++) {
     286          28 :       d = PetscLogReal(csum[cols[j]])/l2;
     287          28 :       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
     288          28 :       d = PetscPowReal(2.0,e);
     289          28 :       aux[cols[j]] = d*d;
     290          28 :       emaxl = PetscMax(emaxl,e);
     291          28 :       eminl = PetscMin(eminl,e);
     292             :     }
     293             :     /* Scale M */
     294          36 :     PetscCall(MatSeqAIJGetArray(M,&array));
     295        1692 :     for (j=0;j<nz;j++) {
     296        1656 :       array[j] *= aux[cidx[j]];
     297             :     }
     298          36 :     PetscCall(MatSeqAIJRestoreArray(M,&array));
     299             :     /* Row sum */
     300          36 :     PetscCall(PetscArrayzero(rsum,nr));
     301          36 :     PetscCall(MatSeqAIJGetArray(M,&array));
     302         628 :     for (i=0;i<nr;i++) {
     303        2248 :       for (j=ridx[i];j<ridx[i+1];j++) rsum[i] += PetscAbsScalar(array[j]);
     304             :       /* Update Dl */
     305         592 :       d = PetscLogReal(rsum[i])/l2;
     306         592 :       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
     307         592 :       d = PetscPowReal(2.0,e);
     308         592 :       Dl[i] *= d;
     309             :       /* Scale M */
     310        2248 :       for (j=ridx[i];j<ridx[i+1];j++) array[j] *= d*d;
     311         592 :       emaxl = PetscMax(emaxl,e);
     312         592 :       eminl = PetscMin(eminl,e);
     313             :     }
     314          36 :     PetscCall(MatSeqAIJRestoreArray(M,&array));
     315             :     /* Compute global max and min */
     316         144 :     PetscCall(MPIU_Allreduce(&emaxl,&emax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pep->Dl)));
     317         144 :     PetscCall(MPIU_Allreduce(&eminl,&emin,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)pep->Dl)));
     318          36 :     if (emax<=emin+2) cont = PETSC_FALSE;
     319             :   }
     320          18 :   PetscCall(VecRestoreArray(pep->Dr,&Dr));
     321          18 :   PetscCall(VecRestoreArray(pep->Dl,&Dl));
     322             :   /* Free memory*/
     323          18 :   PetscCall(MatDestroy(&M));
     324          18 :   PetscCall(PetscFree4(rsum,csum,aux,cols));
     325          18 :   PetscCall(PetscFree(T));
     326          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     327             : }
     328             : 
     329             : /*
     330             :    PEPComputeScaleFactor - compute sfactor as described in [Betcke 2008].
     331             : */
     332         228 : PetscErrorCode PEPComputeScaleFactor(PEP pep)
     333             : {
     334         228 :   PetscBool      has0,has1,flg;
     335         228 :   PetscReal      norm0,norm1;
     336         228 :   Mat            T[2];
     337         228 :   PEPBasis       basis;
     338         228 :   PetscInt       i;
     339             : 
     340         228 :   PetscFunctionBegin;
     341         228 :   if (pep->scale==PEP_SCALE_NONE || pep->scale==PEP_SCALE_DIAGONAL) {  /* no scalar scaling */
     342         219 :     pep->sfactor = 1.0;
     343         219 :     pep->dsfactor = 1.0;
     344         219 :     PetscFunctionReturn(PETSC_SUCCESS);
     345             :   }
     346           9 :   if (pep->sfactor_set) PetscFunctionReturn(PETSC_SUCCESS);  /* user provided value */
     347           8 :   pep->sfactor = 1.0;
     348           8 :   pep->dsfactor = 1.0;
     349           8 :   PetscCall(PEPGetBasis(pep,&basis));
     350           8 :   if (basis==PEP_BASIS_MONOMIAL) {
     351           8 :     PetscCall(STGetTransform(pep->st,&flg));
     352           8 :     if (flg) {
     353           1 :       PetscCall(STGetMatrixTransformed(pep->st,0,&T[0]));
     354           1 :       PetscCall(STGetMatrixTransformed(pep->st,pep->nmat-1,&T[1]));
     355             :     } else {
     356           7 :       T[0] = pep->A[0];
     357           7 :       T[1] = pep->A[pep->nmat-1];
     358             :     }
     359           8 :     if (pep->nmat>2) {
     360           8 :       PetscCall(MatHasOperation(T[0],MATOP_NORM,&has0));
     361           8 :       PetscCall(MatHasOperation(T[1],MATOP_NORM,&has1));
     362           8 :       if (has0 && has1) {
     363           8 :         PetscCall(MatNorm(T[0],NORM_INFINITY,&norm0));
     364           8 :         PetscCall(MatNorm(T[1],NORM_INFINITY,&norm1));
     365           8 :         pep->sfactor = PetscPowReal(norm0/norm1,1.0/(pep->nmat-1));
     366           8 :         pep->dsfactor = norm1;
     367          16 :         for (i=pep->nmat-2;i>0;i--) {
     368           8 :           PetscCall(STGetMatrixTransformed(pep->st,i,&T[1]));
     369           8 :           PetscCall(MatHasOperation(T[1],MATOP_NORM,&has1));
     370           8 :           if (has1) {
     371           8 :             PetscCall(MatNorm(T[1],NORM_INFINITY,&norm1));
     372           8 :             pep->dsfactor = pep->dsfactor*pep->sfactor+norm1;
     373             :           } else break;
     374             :         }
     375           8 :         if (has1) {
     376           8 :           pep->dsfactor = pep->dsfactor*pep->sfactor+norm0;
     377           8 :           pep->dsfactor = pep->nmat/pep->dsfactor;
     378           0 :         } else pep->dsfactor = 1.0;
     379             :       }
     380             :     }
     381             :   }
     382           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     383             : }
     384             : 
     385             : /*
     386             :    PEPBasisCoefficients - compute polynomial basis coefficients
     387             : */
     388         144 : PetscErrorCode PEPBasisCoefficients(PEP pep,PetscReal *pbc)
     389             : {
     390         144 :   PetscReal *ca,*cb,*cg;
     391         144 :   PetscInt  k,nmat=pep->nmat;
     392             : 
     393         144 :   PetscFunctionBegin;
     394         144 :   ca = pbc;
     395         144 :   cb = pbc+nmat;
     396         144 :   cg = pbc+2*nmat;
     397         144 :   switch (pep->basis) {
     398             :   case PEP_BASIS_MONOMIAL:
     399         486 :     for (k=0;k<nmat;k++) {
     400         368 :       ca[k] = 1.0; cb[k] = 0.0; cg[k] = 0.0;
     401             :     }
     402             :     break;
     403          26 :   case PEP_BASIS_CHEBYSHEV1:
     404          26 :     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
     405         204 :     for (k=1;k<nmat;k++) {
     406         178 :       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
     407             :     }
     408             :     break;
     409           0 :   case PEP_BASIS_CHEBYSHEV2:
     410           0 :     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
     411           0 :     for (k=1;k<nmat;k++) {
     412           0 :       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
     413             :     }
     414             :     break;
     415           0 :   case PEP_BASIS_LEGENDRE:
     416           0 :     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
     417           0 :     for (k=1;k<nmat;k++) {
     418           0 :       ca[k] = k+1; cb[k] = -2*k; cg[k] = k;
     419             :     }
     420             :     break;
     421           0 :   case PEP_BASIS_LAGUERRE:
     422           0 :     ca[0] = -1.0; cb[0] = 0.0; cg[0] = 0.0;
     423           0 :     for (k=1;k<nmat;k++) {
     424           0 :       ca[k] = -(k+1); cb[k] = 2*k+1; cg[k] = -k;
     425             :     }
     426             :     break;
     427           0 :   case PEP_BASIS_HERMITE:
     428           0 :     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
     429           0 :     for (k=1;k<nmat;k++) {
     430           0 :       ca[k] = .5; cb[k] = 0.0; cg[k] = -k;
     431             :     }
     432             :     break;
     433             :   }
     434         144 :   PetscFunctionReturn(PETSC_SUCCESS);
     435             : }