Actual source code: pepkrylov.c
slepc-3.18.2 2023-01-26
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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: PetscInt i,j,nq,deg=pep->nmat-1,lds,idxcpy=0,ldds,k,ld;
21: PetscScalar *X,*er,*ei,*SS,*vals,*ivals,sone=1.0,szero=0.0,*yi,*yr,*tr,*ti,alpha,*pS0;
22: const PetscScalar *S;
23: PetscBLASInt k_,nq_,lds_,one=1,ldds_,cols,info,zero=0;
24: PetscBool flg;
25: PetscReal norm,max,t,factor=1.0,done=1.0;
26: Vec xr,xi,w[4];
27: PEP_TOAR *ctx = (PEP_TOAR*)pep->data;
28: Mat S0,MS;
30: BVTensorGetFactors(ctx->V,NULL,&MS);
31: MatDenseGetArrayRead(MS,&S);
32: BVGetSizes(pep->V,NULL,NULL,&ld);
33: BVGetActiveColumns(pep->V,NULL,&nq);
34: k = pep->nconv;
35: if (k==0) return 0;
36: lds = deg*ld;
37: DSGetLeadingDimension(pep->ds,&ldds);
38: PetscCalloc5(k,&er,k,&ei,nq*k,&SS,pep->nmat,&vals,pep->nmat,&ivals);
39: STGetTransform(pep->st,&flg);
40: if (flg) factor = pep->sfactor;
41: for (i=0;i<k;i++) {
42: er[i] = factor*pep->eigr[i];
43: ei[i] = factor*pep->eigi[i];
44: }
45: STBackTransform(pep->st,k,er,ei);
47: DSVectors(pep->ds,DS_MAT_X,NULL,NULL);
48: DSGetArray(pep->ds,DS_MAT_X,&X);
50: PetscBLASIntCast(k,&k_);
51: PetscBLASIntCast(nq,&nq_);
52: PetscBLASIntCast(lds,&lds_);
53: PetscBLASIntCast(ldds,&ldds_);
55: if (pep->extract==PEP_EXTRACT_NONE || pep->refine==PEP_REFINE_MULTIPLE) {
56: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nq_,&k_,&k_,&sone,S,&lds_,X,&ldds_,&szero,SS,&nq_));
57: } else {
58: switch (pep->extract) {
59: case PEP_EXTRACT_NONE:
60: break;
61: case PEP_EXTRACT_NORM:
62: for (i=0;i<k;i++) {
63: PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
64: max = 1.0;
65: for (j=1;j<deg;j++) {
66: norm = SlepcAbsEigenvalue(vals[j],ivals[j]);
67: if (max < norm) { max = norm; idxcpy = j; }
68: }
69: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
70: #if !defined(PETSC_USE_COMPLEX)
71: if (PetscRealPart(ei[i])!=0.0) {
72: i++;
73: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
74: }
75: #endif
76: }
77: break;
78: case PEP_EXTRACT_RESIDUAL:
79: VecDuplicate(pep->work[0],&xr);
80: VecDuplicate(pep->work[0],&w[0]);
81: VecDuplicate(pep->work[0],&w[1]);
82: #if !defined(PETSC_USE_COMPLEX)
83: VecDuplicate(pep->work[0],&w[2]);
84: VecDuplicate(pep->work[0],&w[3]);
85: VecDuplicate(pep->work[0],&xi);
86: #else
87: xi = NULL;
88: #endif
89: for (i=0;i<k;i++) {
90: max = 0.0;
91: for (j=0;j<deg;j++) {
92: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
93: BVMultVec(pep->V,1.0,0.0,xr,SS+i*nq);
94: #if !defined(PETSC_USE_COMPLEX)
95: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+(i+1)*ldds,&one,&szero,SS+i*nq,&one));
96: BVMultVec(pep->V,1.0,0.0,xi,SS+i*nq);
97: #endif
98: PEPComputeResidualNorm_Private(pep,er[i],ei[i],xr,xi,w,&norm);
99: if (norm>max) { max = norm; idxcpy=j; }
100: }
101: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
102: #if !defined(PETSC_USE_COMPLEX)
103: if (PetscRealPart(ei[i])!=0.0) {
104: i++;
105: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
106: }
107: #endif
108: }
109: VecDestroy(&xr);
110: VecDestroy(&w[0]);
111: VecDestroy(&w[1]);
112: #if !defined(PETSC_USE_COMPLEX)
113: VecDestroy(&w[2]);
114: VecDestroy(&w[3]);
115: VecDestroy(&xi);
116: #endif
117: break;
118: case PEP_EXTRACT_STRUCTURED:
119: PetscMalloc2(k,&tr,k,&ti);
120: for (i=0;i<k;i++) {
121: t = 0.0;
122: PEPEvaluateBasis(pep,er[i],ei[i],vals,ivals);
123: yr = X+i*ldds; yi = NULL;
124: #if !defined(PETSC_USE_COMPLEX)
125: if (ei[i]!=0.0) { yr = tr; yi = ti; }
126: #endif
127: for (j=0;j<deg;j++) {
128: alpha = PetscConj(vals[j]);
129: #if !defined(PETSC_USE_COMPLEX)
130: if (ei[i]!=0.0) {
131: PetscArrayzero(tr,k);
132: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
133: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
134: PetscArrayzero(ti,k);
135: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
136: alpha = -ivals[j];
137: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
138: alpha = 1.0;
139: }
140: #endif
141: PetscCallBLAS("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: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
145: }
146: }
147: cols = yi? 2: 1;
148: PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&t,&done,&nq_,&cols,SS+i*nq,&nq_,&info));
149: SlepcCheckLapackInfo("lascl",info);
150: if (yi) i++;
151: }
152: PetscFree2(tr,ti);
153: break;
154: }
155: }
157: /* update vectors V = V*S */
158: MatCreateSeqDense(PETSC_COMM_SELF,nq,k,NULL,&S0);
159: MatDenseGetArrayWrite(S0,&pS0);
160: for (i=0;i<k;i++) PetscArraycpy(pS0+i*nq,SS+i*nq,nq);
161: MatDenseRestoreArrayWrite(S0,&pS0);
162: BVMultInPlace(pep->V,S0,0,k);
163: MatDestroy(&S0);
164: PetscFree5(er,ei,SS,vals,ivals);
165: MatDenseRestoreArrayRead(MS,&S);
166: BVTensorRestoreFactors(ctx->V,NULL,&MS);
167: return 0;
168: }
170: /*
171: PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
172: for polynomial Krylov methods.
174: Differences:
175: - Always non-symmetric
176: - Does not check for STSHIFT
177: - No correction factor
178: - No support for true residual
179: */
180: PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
181: {
182: PetscInt k,newk,marker,inside;
183: PetscScalar re,im;
184: PetscReal resnorm;
185: PetscBool istrivial;
187: RGIsTrivial(pep->rg,&istrivial);
188: marker = -1;
189: if (pep->trackall) getall = PETSC_TRUE;
190: for (k=kini;k<kini+nits;k++) {
191: /* eigenvalue */
192: re = pep->eigr[k];
193: im = pep->eigi[k];
194: if (PetscUnlikely(!istrivial)) {
195: STBackTransform(pep->st,1,&re,&im);
196: RGCheckInside(pep->rg,1,&re,&im,&inside);
197: if (marker==-1 && inside<0) marker = k;
198: re = pep->eigr[k];
199: im = pep->eigi[k];
200: }
201: newk = k;
202: DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm);
203: resnorm *= beta;
204: /* error estimate */
205: (*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx);
206: if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
207: if (PetscUnlikely(newk==k+1)) {
208: pep->errest[k+1] = pep->errest[k];
209: k++;
210: }
211: if (marker!=-1 && !getall) break;
212: }
213: if (marker!=-1) k = marker;
214: *kout = k;
215: return 0;
216: }