Actual source code: peputils.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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:    Common subroutines for several PEP solvers
 12: */

 14: #include <slepc/private/pepimpl.h>    /*I "slepcpep.h" I*/
 15: #include <slepcblaslapack.h>

 17: /* 
 18:   Computes T_j=phy_idx(T) for a matrix T.
 19:     Tp (if j>0) and Tpp (if j>1) are the evaluation 
 20:     of phy_(j-1) and phy(j-2)respectively. 
 21: */
 22: PetscErrorCode PEPEvaluateBasisMat(PEP pep,PetscInt k,PetscScalar *T,PetscInt ldt,PetscInt idx,PetscScalar *Tpp,PetscInt ldtpp,PetscScalar *Tp,PetscInt ldtp,PetscScalar *Tj,PetscInt ldtj)
 23: {
 25:   PetscInt       i;
 26:   PetscReal      *ca,*cb,*cg;
 27:   PetscScalar    g,a;
 28:   PetscBLASInt   k_,ldt_,ldtj_,ldtp_;

 31:   if (idx==0) {
 32:     for (i=0;i<k;i++) {
 33:       PetscMemzero(Tj+i*ldtj,k*sizeof(PetscScalar));
 34:       Tj[i+i*ldtj] = 1.0;
 35:     }
 36:   } else {
 37:     PetscBLASIntCast(ldt,&ldt_);
 38:     PetscBLASIntCast(ldtj,&ldtj_);
 39:     PetscBLASIntCast(ldtp,&ldtp_);
 40:     PetscBLASIntCast(k,&k_);
 41:     ca = pep->pbc; cb = pep->pbc+pep->nmat; cg = pep->pbc+2*pep->nmat;
 42:     for (i=0;i<k;i++) T[i*ldt+i] -= cb[idx-1];
 43:     if (idx>1) { 
 44:       for (i=0;i<k;i++) {
 45:         PetscMemcpy(Tj+i*ldtj,Tpp+i*ldtpp,k*sizeof(PetscScalar));
 46:       }
 47:     }
 48:     a = 1/ca[idx-1];
 49:     g = (idx==1)?0.0:-cg[idx-1]/ca[idx-1];
 50:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&k_,&k_,&k_,&a,T,&ldt_,Tp,&ldtp_,&g,Tj,&ldtj_));
 51:     for (i=0;i<k;i++) T[i*ldt+i] += cb[idx-1];
 52:   }
 53:   return(0);
 54: }

 56: /*
 57:    PEPEvaluateBasis - evaluate the polynomial basis on a given parameter sigma
 58: */
 59: PetscErrorCode PEPEvaluateBasis(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
 60: {
 61:   PetscInt   nmat=pep->nmat,k;
 62:   PetscReal  *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;

 65:   if (ivals) for (k=0;k<nmat;k++) ivals[k] = 0.0;
 66:   vals[0] = 1.0;
 67:   vals[1] = (sigma-b[0])/a[0];
 68: #if !defined(PETSC_USE_COMPLEX)
 69:   if (ivals) ivals[1] = isigma/a[0];
 70: #endif
 71:   for (k=2;k<nmat;k++) {
 72:     vals[k] = ((sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2])/a[k-1];
 73:     if (ivals) vals[k] -= isigma*ivals[k-1]/a[k-1];
 74: #if !defined(PETSC_USE_COMPLEX)
 75:     if (ivals) ivals[k] = ((sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2])/a[k-1];
 76: #endif
 77:   }
 78:   return(0);
 79: }

 81: /*
 82:    PEPEvaluateBasisDerivative - evaluate the derivative of the polynomial basis on a given parameter sigma
 83: */
 84: PetscErrorCode PEPEvaluateBasisDerivative(PEP pep,PetscScalar sigma,PetscScalar isigma,PetscScalar *vals,PetscScalar *ivals)
 85: {
 87:   PetscInt       nmat=pep->nmat,k;
 88:   PetscReal      *a=pep->pbc,*b=pep->pbc+nmat,*g=pep->pbc+2*nmat;

 91:   PEPEvaluateBasis(pep,sigma,isigma,vals,ivals);
 92:   for (k=nmat-1;k>0;k--) {
 93:     vals[k] = vals[k-1];
 94:     if (ivals) ivals[k] = ivals[k-1];
 95:   }
 96:   vals[0] = 0.0;
 97:   vals[1] = vals[1]/a[0];
 98: #if !defined(PETSC_USE_COMPLEX)
 99:   if (ivals) ivals[1] = ivals[1]/a[0];
100: #endif
101:   for (k=2;k<nmat;k++) {
102:     vals[k] += (sigma-b[k-1])*vals[k-1]-g[k-1]*vals[k-2];
103:     if (ivals) vals[k] -= isigma*ivals[k-1];
104:     vals[k] /= a[k-1];
105: #if !defined(PETSC_USE_COMPLEX)
106:     if (ivals) {
107:       ivals[k] += (sigma-b[k-1])*ivals[k-1]+isigma*vals[k-1]-g[k-1]*ivals[k-2];
108:       ivals[k] /= a[k-1];
109:     }
110: #endif
111:   }
112:   return(0);
113: }