Actual source code: pepkrylov.c
slepc-3.13.2 2020-05-12
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: {
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: #if !defined(PETSC_USE_COMPLEX)
126: if (ei[i]!=0.0) { yr = tr; yi = ti; }
127: #endif
128: for (j=0;j<deg;j++) {
129: alpha = PetscConj(vals[j]);
130: #if !defined(PETSC_USE_COMPLEX)
131: if (ei[i]!=0.0) {
132: PetscArrayzero(tr,k);
133: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
134: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
135: PetscArrayzero(ti,k);
136: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
137: alpha = -ivals[j];
138: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
139: alpha = 1.0;
140: }
141: #endif
142: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yr,&one,&sone,SS+i*nq,&one));
143: t += SlepcAbsEigenvalue(vals[j],ivals[j])*SlepcAbsEigenvalue(vals[j],ivals[j]);
144: if (yi) {
145: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
146: }
147: }
148: t = 1.0/t;
149: PetscStackCallBLAS("BLASscal",BLASscal_(&nq_,&t,SS+i*nq,&one));
150: if (yi) {
151: PetscStackCallBLAS("BLASscal",BLASscal_(&nq_,&t,SS+(i+1)*nq,&one));
152: i++;
153: }
154: }
155: PetscFree2(tr,ti);
156: break;
157: }
158: }
160: /* update vectors V = V*S */
161: MatCreateSeqDense(PETSC_COMM_SELF,nq,k,NULL,&S0);
162: MatDenseGetArray(S0,&pS0);
163: for (i=0;i<k;i++) {
164: PetscArraycpy(pS0+i*nq,SS+i*nq,nq);
165: }
166: MatDenseRestoreArray(S0,&pS0);
167: BVMultInPlace(pep->V,S0,0,k);
168: MatDestroy(&S0);
169: PetscFree5(er,ei,SS,vals,ivals);
170: if (ctx->V) {
171: MatDenseRestoreArray(MS,&S);
172: BVTensorRestoreFactors(ctx->V,NULL,&MS);
173: }
174: return(0);
175: }
177: /*
178: PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
179: for polynomial Krylov methods.
181: Differences:
182: - Always non-symmetric
183: - Does not check for STSHIFT
184: - No correction factor
185: - No support for true residual
186: */
187: PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
188: {
190: PetscInt k,newk,marker,inside;
191: PetscScalar re,im;
192: PetscReal resnorm;
193: PetscBool istrivial;
196: RGIsTrivial(pep->rg,&istrivial);
197: marker = -1;
198: if (pep->trackall) getall = PETSC_TRUE;
199: for (k=kini;k<kini+nits;k++) {
200: /* eigenvalue */
201: re = pep->eigr[k];
202: im = pep->eigi[k];
203: if (!istrivial) {
204: STBackTransform(pep->st,1,&re,&im);
205: RGCheckInside(pep->rg,1,&re,&im,&inside);
206: if (marker==-1 && inside<0) marker = k;
207: re = pep->eigr[k];
208: im = pep->eigi[k];
209: }
210: newk = k;
211: DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm);
212: resnorm *= beta;
213: /* error estimate */
214: (*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx);
215: if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
216: if (newk==k+1) {
217: pep->errest[k+1] = pep->errest[k];
218: k++;
219: }
220: if (marker!=-1 && !getall) break;
221: }
222: if (marker!=-1) k = marker;
223: *kout = k;
224: return(0);
225: }