Actual source code: pepkrylov.c
slepc-3.21.1 2024-04-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: PetscFunctionBegin;
31: PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
32: PetscCall(MatDenseGetArrayRead(MS,&S));
33: PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
34: PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
35: k = pep->nconv;
36: if (k==0) PetscFunctionReturn(PETSC_SUCCESS);
37: lds = deg*ld;
38: PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
39: PetscCall(PetscCalloc5(k,&er,k,&ei,nq*k,&SS,pep->nmat,&vals,pep->nmat,&ivals));
40: PetscCall(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: PetscCall(STBackTransform(pep->st,k,er,ei));
48: PetscCall(DSVectors(pep->ds,DS_MAT_X,NULL,NULL));
49: PetscCall(DSGetArray(pep->ds,DS_MAT_X,&X));
51: PetscCall(PetscBLASIntCast(k,&k_));
52: PetscCall(PetscBLASIntCast(nq,&nq_));
53: PetscCall(PetscBLASIntCast(lds,&lds_));
54: PetscCall(PetscBLASIntCast(ldds,&ldds_));
56: if (pep->extract==PEP_EXTRACT_NONE || pep->refine==PEP_REFINE_MULTIPLE) {
57: PetscCallBLAS("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: PetscCall(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: PetscCallBLAS("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: PetscCallBLAS("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: PetscCall(VecDuplicate(pep->work[0],&xr));
81: PetscCall(VecDuplicate(pep->work[0],&w[0]));
82: PetscCall(VecDuplicate(pep->work[0],&w[1]));
83: #if !defined(PETSC_USE_COMPLEX)
84: PetscCall(VecDuplicate(pep->work[0],&w[2]));
85: PetscCall(VecDuplicate(pep->work[0],&w[3]));
86: PetscCall(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: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
94: PetscCall(BVMultVec(pep->V,1.0,0.0,xr,SS+i*nq));
95: #if !defined(PETSC_USE_COMPLEX)
96: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+j*ld,&lds_,X+(i+1)*ldds,&one,&szero,SS+i*nq,&one));
97: PetscCall(BVMultVec(pep->V,1.0,0.0,xi,SS+i*nq));
98: #endif
99: PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],xr,xi,w,&norm));
100: if (norm>max) { max = norm; idxcpy=j; }
101: }
102: PetscCallBLAS("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: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&sone,S+idxcpy*ld,&lds_,X+i*ldds,&one,&szero,SS+i*nq,&one));
107: }
108: #endif
109: }
110: PetscCall(VecDestroy(&xr));
111: PetscCall(VecDestroy(&w[0]));
112: PetscCall(VecDestroy(&w[1]));
113: #if !defined(PETSC_USE_COMPLEX)
114: PetscCall(VecDestroy(&w[2]));
115: PetscCall(VecDestroy(&w[3]));
116: PetscCall(VecDestroy(&xi));
117: #endif
118: break;
119: case PEP_EXTRACT_STRUCTURED:
120: PetscCall(PetscMalloc2(k,&tr,k,&ti));
121: for (i=0;i<k;i++) {
122: t = 0.0;
123: PetscCall(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: PetscCall(PetscArrayzero(tr,k));
133: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+i*ldds,&one,tr,&one));
134: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&ivals[j],X+(i+1)*ldds,&one,tr,&one));
135: PetscCall(PetscArrayzero(ti,k));
136: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&vals[j],X+(i+1)*ldds,&one,ti,&one));
137: alpha = -ivals[j];
138: PetscCallBLAS("BLASaxpy",BLASaxpy_(&k_,&alpha,X+i*ldds,&one,ti,&one));
139: alpha = 1.0;
140: }
141: #endif
142: PetscCallBLAS("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: PetscCallBLAS("BLASgemv",BLASgemv_("N",&nq_,&k_,&alpha,S+j*ld,&lds_,yi,&one,&sone,SS+(i+1)*nq,&one));
146: }
147: }
148: cols = yi? 2: 1;
149: PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&t,&done,&nq_,&cols,SS+i*nq,&nq_,&info));
150: SlepcCheckLapackInfo("lascl",info);
151: if (yi) i++;
152: }
153: PetscCall(PetscFree2(tr,ti));
154: break;
155: }
156: }
158: /* update vectors V = V*S */
159: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,nq,k,NULL,&S0));
160: PetscCall(MatDenseGetArrayWrite(S0,&pS0));
161: for (i=0;i<k;i++) PetscCall(PetscArraycpy(pS0+i*nq,SS+i*nq,nq));
162: PetscCall(MatDenseRestoreArrayWrite(S0,&pS0));
163: PetscCall(BVMultInPlace(pep->V,S0,0,k));
164: PetscCall(MatDestroy(&S0));
165: PetscCall(PetscFree5(er,ei,SS,vals,ivals));
166: PetscCall(MatDenseRestoreArrayRead(MS,&S));
167: PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*
172: PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
173: for polynomial Krylov methods.
175: Differences:
176: - Always non-symmetric
177: - Does not check for STSHIFT
178: - No correction factor
179: - No support for true residual
180: */
181: PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
182: {
183: PetscInt k,newk,marker,inside;
184: PetscScalar re,im;
185: PetscReal resnorm;
186: PetscBool istrivial;
188: PetscFunctionBegin;
189: PetscCall(RGIsTrivial(pep->rg,&istrivial));
190: marker = -1;
191: if (pep->trackall) getall = PETSC_TRUE;
192: for (k=kini;k<kini+nits;k++) {
193: /* eigenvalue */
194: re = pep->eigr[k];
195: im = pep->eigi[k];
196: if (PetscUnlikely(!istrivial)) {
197: PetscCall(STBackTransform(pep->st,1,&re,&im));
198: PetscCall(RGCheckInside(pep->rg,1,&re,&im,&inside));
199: if (marker==-1 && inside<0) marker = k;
200: re = pep->eigr[k];
201: im = pep->eigi[k];
202: }
203: newk = k;
204: PetscCall(DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm));
205: resnorm *= beta;
206: /* error estimate */
207: PetscCall((*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx));
208: if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
209: if (PetscUnlikely(newk==k+1)) {
210: pep->errest[k+1] = pep->errest[k];
211: k++;
212: }
213: if (marker!=-1 && !getall) break;
214: }
215: if (marker!=-1) k = marker;
216: *kout = k;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }