Actual source code: pepkrylov.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 all Krylov-type PEP solvers
 12: */

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

 18: PetscErrorCode PEPExtractVectors_TOAR(PEP pep)
 19: {
 21:   PetscInt       i,j,nq,deg=pep->nmat-1,lds,idxcpy=0,ldds,k,ld;
 22:   PetscScalar    *X,*er,*ei,*SS,*vals,*ivals,sone=1.0,szero=0.0,*yi,*yr,*tr,*ti,alpha,t,*S,*pS0;
 23:   PetscBLASInt   k_,nq_,lds_,one=1,ldds_;
 24:   PetscBool      flg;
 25:   PetscReal      norm,max,factor=1.0;
 26:   Vec            xr,xi,w[4];
 27:   PEP_TOAR       *ctx = (PEP_TOAR*)pep->data;
 28:   Mat            S0,MS;

 31:   BVTensorGetFactors(ctx->V,NULL,&MS);
 32:   MatDenseGetArray(MS,&S);
 33:   BVGetSizes(pep->V,NULL,NULL,&ld);
 34:   BVGetActiveColumns(pep->V,NULL,&nq);
 35:   k = pep->nconv;
 36:   if (k==0) return(0);
 37:   lds = deg*ld;
 38:   DSGetLeadingDimension(pep->ds,&ldds);
 39:   PetscCalloc5(k,&er,k,&ei,nq*k,&SS,pep->nmat,&vals,pep->nmat,&ivals);
 40:   STGetTransform(pep->st,&flg);
 41:   if (flg) factor = pep->sfactor;
 42:   for (i=0;i<k;i++) {
 43:     er[i] = factor*pep->eigr[i];
 44:     ei[i] = factor*pep->eigi[i];
 45:   }
 46:   STBackTransform(pep->st,k,er,ei);

 48:   DSVectors(pep->ds,DS_MAT_X,NULL,NULL);
 49:   DSGetArray(pep->ds,DS_MAT_X,&X);

 51:   PetscBLASIntCast(k,&k_);
 52:   PetscBLASIntCast(nq,&nq_);
 53:   PetscBLASIntCast(lds,&lds_);
 54:   PetscBLASIntCast(ldds,&ldds_);

 56:   if (pep->extract==PEP_EXTRACT_NONE || pep->refine==PEP_REFINE_MULTIPLE) {
 57:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nq_,&k_,&k_,&sone,S,&lds_,X,&ldds_,&szero,SS,&nq_));
 58:   } else {
 59:     switch (pep->extract) {
 60:     case PEP_EXTRACT_NONE:
 61:       break;
 62:     case PEP_EXTRACT_NORM:
 63:       for (i=0;i<k;i++) {
 64:         PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
 65:         max = 1.0;
 66:         for (j=1;j<deg;j++) {
 67:           norm = SlepcAbsEigenvalue(vals[j],ivals[j]);
 68:           if (max < norm) { max = norm; idxcpy = j; }
 69:         }
 70:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 71: #if !defined(PETSC_USE_COMPLEX)
 72:         if (PetscRealPart(ei[i])!=0.0) {
 73:           i++;
 74:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 75:         }
 76: #endif
 77:       }
 78:       break;
 79:     case PEP_EXTRACT_RESIDUAL:
 80:       VecDuplicate(pep->work[0],&xr);
 81:       VecDuplicate(pep->work[0],&w[0]);
 82:       VecDuplicate(pep->work[0],&w[1]);
 83: #if !defined(PETSC_USE_COMPLEX)
 84:       VecDuplicate(pep->work[0],&w[2]);
 85:       VecDuplicate(pep->work[0],&w[3]);
 86:       VecDuplicate(pep->work[0],&xi);
 87: #else
 88:       xi = NULL;
 89: #endif
 90:       for (i=0;i<k;i++) {
 91:         max = 0.0;
 92:         for (j=0;j<deg;j++) {
 93:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 94:           BVMultVec(pep->V,1.0,0.0,xr,SS+i*nq);
 95: #if !defined(PETSC_USE_COMPLEX)
 96:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+(i+1)*ldds,&one,&szero,SS+i*nq,&one));
 97:           BVMultVec(pep->V,1.0,0.0,xi,SS+i*nq);
 98: #endif
 99:           PEPComputeResidualNorm_Private(pep,er[i],ei[i],xr,xi,w,&norm);
100:           if (norm>max) { max = norm; idxcpy=j; }
101:         }
102:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
103: #if !defined(PETSC_USE_COMPLEX)
104:         if (PetscRealPart(ei[i])!=0.0) {
105:           i++;
106:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
107:         }
108: #endif
109:       }
110:       VecDestroy(&xr);
111:       VecDestroy(&w[0]);
112:       VecDestroy(&w[1]);
113: #if !defined(PETSC_USE_COMPLEX)
114:       VecDestroy(&w[2]);
115:       VecDestroy(&w[3]);
116:       VecDestroy(&xi);
117: #endif
118:       break;
119:     case PEP_EXTRACT_STRUCTURED:
120:       PetscMalloc2(k,&tr,k,&ti);
121:       for (i=0;i<k;i++) {
122:         t = 0.0;
123:         PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
124:         yr = X+i*ldds; yi = NULL;
125:         for (j=0;j<deg;j++) {
126:           alpha = PetscConj(vals[j]);
127: #if !defined(PETSC_USE_COMPLEX)
128:           if (ei[i]!=0.0) {
129:             PetscMemzero(tr,k*sizeof(PetscScalar));
130:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
131:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
132:             yr = tr;
133:             PetscMemzero(ti,k*sizeof(PetscScalar));
134:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
135:             alpha = -ivals[j];
136:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
137:             yi = ti;
138:             alpha = 1.0;
139:           } else { yr = X+i*ldds; yi = NULL; }
140: #endif
141:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yr,&one,&sone,SS+i*nq,&one));
142:           t += SlepcAbsEigenvalue(vals[j],ivals[j])*SlepcAbsEigenvalue(vals[j],ivals[j]);
143:           if (yi) {
144:             PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
145:           }
146:         }
147:         t = 1.0/t;
148:         PetscStackCallBLAS("BLASscal",BLASscal_(&nq_,&t,SS+i*nq,&one));
149:         if (yi) {
150:           PetscStackCallBLAS("BLASscal",BLASscal_(&nq_,&t,SS+(i+1)*nq,&one));
151:           i++;
152:         }
153:       }
154:       PetscFree2(tr,ti);
155:       break;
156:     }
157:   }

159:   /* update vectors V = V*S */
160:   MatCreateSeqDense(PETSC_COMM_SELF,nq,k,NULL,&S0);
161:   MatDenseGetArray(S0,&pS0);
162:   for (i=0;i<k;i++) {
163:     PetscMemcpy(pS0+i*nq,SS+i*nq,nq*sizeof(PetscScalar));
164:   }
165:   MatDenseRestoreArray(S0,&pS0);
166:   BVMultInPlace(pep->V,S0,0,k);
167:   MatDestroy(&S0);
168:   PetscFree5(er,ei,SS,vals,ivals);
169:   if (ctx->V) {
170:     MatDenseRestoreArray(MS,&S);
171:     BVTensorRestoreFactors(ctx->V,NULL,&MS);
172:   }
173:   return(0);
174: }

176: /*
177:    PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
178:    for polynomial Krylov methods.

180:    Differences:
181:    - Always non-symmetric
182:    - Does not check for STSHIFT
183:    - No correction factor
184:    - No support for true residual
185: */
186: PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
187: {
189:   PetscInt       k,newk,marker,inside;
190:   PetscScalar    re,im;
191:   PetscReal      resnorm;
192:   PetscBool      istrivial;

195:   RGIsTrivial(pep->rg,&istrivial);
196:   marker = -1;
197:   if (pep->trackall) getall = PETSC_TRUE;
198:   for (k=kini;k<kini+nits;k++) {
199:     /* eigenvalue */
200:     re = pep->eigr[k];
201:     im = pep->eigi[k];
202:     if (!istrivial) {
203:       STBackTransform(pep->st,1,&re,&im);
204:       RGCheckInside(pep->rg,1,&re,&im,&inside);
205:       if (marker==-1 && inside<0) marker = k;
206:       re = pep->eigr[k];
207:       im = pep->eigi[k];
208:     }
209:     newk = k;
210:     DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm);
211:     resnorm *= beta;
212:     /* error estimate */
213:     (*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx);
214:     if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
215:     if (newk==k+1) {
216:       pep->errest[k+1] = pep->errest[k];
217:       k++;
218:     }
219:     if (marker!=-1 && !getall) break;
220:   }
221:   if (marker!=-1) k = marker;
222:   *kout = k;
223:   return(0);
224: }