Actual source code: pepkrylov.c
slepc-3.11.2 2019-07-30
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: }