Actual source code: pepkrylov.c

slepc-3.16.0 2021-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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>
 15: #include <slepcblaslapack.h>
 16: #include "pepkrylov.h"

 18: PetscErrorCode PEPExtractVectors_TOAR(PEP pep)
 19: {
 20:   PetscErrorCode    ierr;
 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,*pS0;
 23:   const PetscScalar *S;
 24:   PetscBLASInt      k_,nq_,lds_,one=1,ldds_,cols,info,zero=0;
 25:   PetscBool         flg;
 26:   PetscReal         norm,max,t,factor=1.0,done=1.0;
 27:   Vec               xr,xi,w[4];
 28:   PEP_TOAR          *ctx = (PEP_TOAR*)pep->data;
 29:   Mat               S0,MS;

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

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

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

 57:   if (pep->extract==PEP_EXTRACT_NONE || pep->refine==PEP_REFINE_MULTIPLE) {
 58:     PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nq_,&k_,&k_,&sone,S,&lds_,X,&ldds_,&szero,SS,&nq_));
 59:   } else {
 60:     switch (pep->extract) {
 61:     case PEP_EXTRACT_NONE:
 62:       break;
 63:     case PEP_EXTRACT_NORM:
 64:       for (i=0;i<k;i++) {
 65:         PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
 66:         max = 1.0;
 67:         for (j=1;j<deg;j++) {
 68:           norm = SlepcAbsEigenvalue(vals[j],ivals[j]);
 69:           if (max < norm) { max = norm; idxcpy = j; }
 70:         }
 71:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 72: #if !defined(PETSC_USE_COMPLEX)
 73:         if (PetscRealPart(ei[i])!=0.0) {
 74:           i++;
 75:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 76:         }
 77: #endif
 78:       }
 79:       break;
 80:     case PEP_EXTRACT_RESIDUAL:
 81:       VecDuplicate(pep->work[0],&xr);
 82:       VecDuplicate(pep->work[0],&w[0]);
 83:       VecDuplicate(pep->work[0],&w[1]);
 84: #if !defined(PETSC_USE_COMPLEX)
 85:       VecDuplicate(pep->work[0],&w[2]);
 86:       VecDuplicate(pep->work[0],&w[3]);
 87:       VecDuplicate(pep->work[0],&xi);
 88: #else
 89:       xi = NULL;
 90: #endif
 91:       for (i=0;i<k;i++) {
 92:         max = 0.0;
 93:         for (j=0;j<deg;j++) {
 94:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
 95:           BVMultVec(pep->V,1.0,0.0,xr,SS+i*nq);
 96: #if !defined(PETSC_USE_COMPLEX)
 97:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+(i+1)*ldds,&one,&szero,SS+i*nq,&one));
 98:           BVMultVec(pep->V,1.0,0.0,xi,SS+i*nq);
 99: #endif
100:           PEPComputeResidualNorm_Private(pep,er[i],ei[i],xr,xi,w,&norm);
101:           if (norm>max) { max = norm; idxcpy=j; }
102:         }
103:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
104: #if !defined(PETSC_USE_COMPLEX)
105:         if (PetscRealPart(ei[i])!=0.0) {
106:           i++;
107:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
108:         }
109: #endif
110:       }
111:       VecDestroy(&xr);
112:       VecDestroy(&w[0]);
113:       VecDestroy(&w[1]);
114: #if !defined(PETSC_USE_COMPLEX)
115:       VecDestroy(&w[2]);
116:       VecDestroy(&w[3]);
117:       VecDestroy(&xi);
118: #endif
119:       break;
120:     case PEP_EXTRACT_STRUCTURED:
121:       PetscMalloc2(k,&tr,k,&ti);
122:       for (i=0;i<k;i++) {
123:         t = 0.0;
124:         PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
125:         yr = X+i*ldds; yi = NULL;
126: #if !defined(PETSC_USE_COMPLEX)
127:         if (ei[i]!=0.0) { yr = tr; yi = ti; }
128: #endif
129:         for (j=0;j<deg;j++) {
130:           alpha = PetscConj(vals[j]);
131: #if !defined(PETSC_USE_COMPLEX)
132:           if (ei[i]!=0.0) {
133:             PetscArrayzero(tr,k);
134:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
135:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
136:             PetscArrayzero(ti,k);
137:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
138:             alpha = -ivals[j];
139:             PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
140:             alpha = 1.0;
141:           }
142: #endif
143:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yr,&one,&sone,SS+i*nq,&one));
144:           t += SlepcAbsEigenvalue(vals[j],ivals[j])*SlepcAbsEigenvalue(vals[j],ivals[j]);
145:           if (yi) {
146:             PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
147:           }
148:         }
149:         cols = yi? 2: 1;
150:         PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&t,&done,&nq_,&cols,SS+i*nq,&nq_,&info));
151:         SlepcCheckLapackInfo("lascl",info);
152:         if (yi) i++;
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:   MatDenseGetArrayWrite(S0,&pS0);
162:   for (i=0;i<k;i++) {
163:     PetscArraycpy(pS0+i*nq,SS+i*nq,nq);
164:   }
165:   MatDenseRestoreArrayWrite(S0,&pS0);
166:   BVMultInPlace(pep->V,S0,0,k);
167:   MatDestroy(&S0);
168:   PetscFree5(er,ei,SS,vals,ivals);
169:   MatDenseRestoreArrayRead(MS,&S);
170:   BVTensorRestoreFactors(ctx->V,NULL,&MS);
171:   return(0);
172: }

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

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

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