Actual source code: dvdcalcpairs.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: SLEPc eigensolver: "davidson"
13: Step: calculate the best eigenpairs in the subspace V
15: For that, performs these steps:
16: 1) Update W <- A * V
17: 2) Update H <- V' * W
18: 3) Obtain eigenpairs of H
19: 4) Select some eigenpairs
20: 5) Compute the Ritz pairs of the selected ones
21: */
23: #include "davidson.h"
24: #include <slepcblaslapack.h>
26: static PetscErrorCode dvd_calcpairs_qz_start(dvdDashboard *d)
27: {
28: BVSetActiveColumns(d->eps->V,0,0);
29: if (d->W) BVSetActiveColumns(d->W,0,0);
30: BVSetActiveColumns(d->AX,0,0);
31: if (d->BX) BVSetActiveColumns(d->BX,0,0);
32: return 0;
33: }
35: static PetscErrorCode dvd_calcpairs_qz_d(dvdDashboard *d)
36: {
37: BVDestroy(&d->W);
38: BVDestroy(&d->AX);
39: BVDestroy(&d->BX);
40: BVDestroy(&d->auxBV);
41: MatDestroy(&d->H);
42: MatDestroy(&d->G);
43: MatDestroy(&d->auxM);
44: SlepcVecPoolDestroy(&d->auxV);
45: PetscFree(d->nBds);
46: return 0;
47: }
49: /* in complex, d->size_H real auxiliary values are needed */
50: static PetscErrorCode dvd_calcpairs_projeig_solve(dvdDashboard *d)
51: {
52: Vec v;
53: Mat A,B,H0,G0;
54: PetscScalar *pA;
55: const PetscScalar *pv;
56: PetscInt i,lV,kV,n,ld;
58: BVGetActiveColumns(d->eps->V,&lV,&kV);
59: n = kV-lV;
60: DSSetDimensions(d->eps->ds,n,0,0);
61: DSGetMat(d->eps->ds,DS_MAT_A,&A);
62: MatDenseGetSubMatrix(d->H,lV,lV+n,lV,lV+n,&H0);
63: MatCopy(H0,A,SAME_NONZERO_PATTERN);
64: MatDenseRestoreSubMatrix(d->H,&H0);
65: DSRestoreMat(d->eps->ds,DS_MAT_A,&A);
66: if (d->G) {
67: DSGetMat(d->eps->ds,DS_MAT_B,&B);
68: MatDenseGetSubMatrix(d->G,lV,lV+n,lV,lV+n,&G0);
69: MatCopy(G0,B,SAME_NONZERO_PATTERN);
70: MatDenseRestoreSubMatrix(d->G,&G0);
71: DSRestoreMat(d->eps->ds,DS_MAT_B,&B);
72: }
73: /* Set the signature on projected matrix B */
74: if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
75: DSGetLeadingDimension(d->eps->ds,&ld);
76: DSGetArray(d->eps->ds,DS_MAT_B,&pA);
77: PetscArrayzero(pA,n*ld);
78: VecCreateSeq(PETSC_COMM_SELF,kV,&v);
79: BVGetSignature(d->eps->V,v);
80: VecGetArrayRead(v,&pv);
81: for (i=0;i<n;i++) {
82: pA[i+ld*i] = d->nBds[i] = PetscRealPart(pv[lV+i]);
83: }
84: VecRestoreArrayRead(v,&pv);
85: VecDestroy(&v);
86: DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);
87: }
88: DSSetState(d->eps->ds,DS_STATE_RAW);
89: DSSolve(d->eps->ds,d->eigr,d->eigi);
90: return 0;
91: }
93: /*
94: A(lA:kA-1,lA:kA-1) <- Z(l:k-1)'*A(l:k-1,l:k-1)*Q(l,k-1), where k=l+kA-lA
95: */
96: static PetscErrorCode EPSXDUpdateProj(Mat Q,Mat Z,PetscInt l,Mat A,PetscInt lA,PetscInt kA,Mat aux)
97: {
98: PetscScalar one=1.0,zero=0.0;
99: PetscInt i,j,dA_=kA-lA,m0,n0,ldA_,ldQ_,ldZ_,nQ_;
100: PetscBLASInt dA,nQ,ldA,ldQ,ldZ;
101: PetscScalar *pA,*pW;
102: const PetscScalar *pQ,*pZ;
103: PetscBool symm=PETSC_FALSE,set,flg;
105: MatGetSize(A,&m0,&n0);
106: MatDenseGetLDA(A,&ldA_);
107: PetscAssert(m0==n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"A should be square");
108: PetscAssert(lA>=0 && lA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial row, column in A");
109: PetscAssert(kA>=0 && kA>=lA && kA<=m0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid final row, column in A");
110: MatIsHermitianKnown(A,&set,&flg);
111: symm = set? flg: PETSC_FALSE;
112: MatGetSize(Q,&m0,&n0); nQ_=m0;
113: MatDenseGetLDA(Q,&ldQ_);
114: PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Q");
115: MatGetSize(Z,&m0,&n0);
116: MatDenseGetLDA(Z,&ldZ_);
117: PetscAssert(l>=0 && l<=n0 && l+dA_<=n0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid initial column in Z");
118: MatGetSize(aux,&m0,&n0);
119: PetscAssert(m0*n0>=nQ_*dA_,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"aux should be larger");
120: PetscBLASIntCast(dA_,&dA);
121: PetscBLASIntCast(nQ_,&nQ);
122: PetscBLASIntCast(ldA_,&ldA);
123: PetscBLASIntCast(ldQ_,&ldQ);
124: PetscBLASIntCast(ldZ_,&ldZ);
125: MatDenseGetArray(A,&pA);
126: MatDenseGetArrayRead(Q,&pQ);
127: if (Q!=Z) MatDenseGetArrayRead(Z,&pZ);
128: else pZ = pQ;
129: MatDenseGetArrayWrite(aux,&pW);
130: /* W = A*Q */
131: if (symm) {
132: /* symmetrize before multiplying */
133: for (i=lA+1;i<lA+nQ;i++) {
134: for (j=lA;j<i;j++) pA[i+j*ldA] = PetscConj(pA[j+i*ldA]);
135: }
136: }
137: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nQ,&dA,&nQ,&one,&pA[ldA*lA+lA],&ldA,&pQ[ldQ*l+l],&ldQ,&zero,pW,&nQ));
138: /* A = Q'*W */
139: PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&dA,&dA,&nQ,&one,&pZ[ldZ*l+l],&ldZ,pW,&nQ,&zero,&pA[ldA*lA+lA],&ldA));
140: MatDenseRestoreArray(A,&pA);
141: MatDenseRestoreArrayRead(Q,&pQ);
142: if (Q!=Z) MatDenseRestoreArrayRead(Z,&pZ);
143: MatDenseRestoreArrayWrite(aux,&pW);
144: return 0;
145: }
147: static PetscErrorCode dvd_calcpairs_updateproj(dvdDashboard *d)
148: {
149: Mat Q,Z;
150: PetscInt lV,kV;
151: PetscBool symm;
153: DSGetMat(d->eps->ds,DS_MAT_Q,&Q);
154: if (d->W) DSGetMat(d->eps->ds,DS_MAT_Z,&Z);
155: else Z = Q;
156: BVGetActiveColumns(d->eps->V,&lV,&kV);
157: EPSXDUpdateProj(Q,Z,0,d->H,lV,lV+d->V_tra_e,d->auxM);
158: if (d->G) EPSXDUpdateProj(Q,Z,0,d->G,lV,lV+d->V_tra_e,d->auxM);
159: DSRestoreMat(d->eps->ds,DS_MAT_Q,&Q);
160: if (d->W) DSRestoreMat(d->eps->ds,DS_MAT_Z,&Z);
162: PetscObjectTypeCompareAny((PetscObject)d->eps->ds,&symm,DSHEP,DSGHIEP,DSGHEP,"");
163: if (d->V_tra_s==0 || symm) return 0;
164: /* Compute upper part of H (and G): H(0:l-1,l:k-1) <- W(0:l-1)' * AV(l:k-1), where
165: k=l+d->V_tra_s */
166: BVSetActiveColumns(d->W?d->W:d->eps->V,0,lV);
167: BVSetActiveColumns(d->AX,lV,lV+d->V_tra_s);
168: BVDot(d->AX,d->W?d->W:d->eps->V,d->H);
169: if (d->G) {
170: BVSetActiveColumns(d->BX?d->BX:d->eps->V,lV,lV+d->V_tra_s);
171: BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);
172: }
173: PetscObjectTypeCompare((PetscObject)d->eps->ds,DSGHEP,&symm);
174: if (!symm) {
175: BVSetActiveColumns(d->W?d->W:d->eps->V,lV,lV+d->V_tra_s);
176: BVSetActiveColumns(d->AX,0,lV);
177: BVDot(d->AX,d->W?d->W:d->eps->V,d->H);
178: if (d->G) {
179: BVSetActiveColumns(d->BX?d->BX:d->eps->V,0,lV);
180: BVDot(d->BX?d->BX:d->eps->V,d->W?d->W:d->eps->V,d->G);
181: }
182: }
183: BVSetActiveColumns(d->eps->V,lV,kV);
184: BVSetActiveColumns(d->AX,lV,kV);
185: if (d->BX) BVSetActiveColumns(d->BX,lV,kV);
186: if (d->W) BVSetActiveColumns(d->W,lV,kV);
187: if (d->W) dvd_harm_updateproj(d);
188: return 0;
189: }
191: /*
192: BV <- BV*MT
193: */
194: static inline PetscErrorCode dvd_calcpairs_updateBV0_gen(dvdDashboard *d,BV bv,DSMatType mat)
195: {
196: PetscInt l,k,n;
197: Mat M,M0,auxM,auxM0;
199: BVGetActiveColumns(d->eps->V,&l,&k);
200: DSGetDimensions(d->eps->ds,&n,NULL,NULL,NULL);
201: PetscAssert(k-l==n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
202: DSGetMat(d->eps->ds,mat,&M);
203: MatDenseGetSubMatrix(M,0,n,0,d->V_tra_e,&M0);
204: MatCreateSeqDense(PETSC_COMM_SELF,k,k,NULL,&auxM);
205: MatDenseGetSubMatrix(auxM,l,l+n,l,l+d->V_tra_e,&auxM0);
206: MatCopy(M0,auxM0,SAME_NONZERO_PATTERN);
207: MatDenseRestoreSubMatrix(auxM,&auxM0);
208: MatDenseRestoreSubMatrix(M,&M0);
209: DSRestoreMat(d->eps->ds,mat,&M);
210: BVMultInPlace(bv,auxM,l,l+d->V_tra_e);
211: MatDestroy(&auxM);
212: return 0;
213: }
215: static PetscErrorCode dvd_calcpairs_proj(dvdDashboard *d)
216: {
217: PetscInt i,l,k;
218: Vec v1,v2;
219: PetscScalar *pv;
221: BVGetActiveColumns(d->eps->V,&l,&k);
222: /* Update AV, BV, W and the projected matrices */
223: /* 1. S <- S*MT */
224: if (d->V_tra_s != d->V_tra_e || d->V_tra_e > 0) {
225: dvd_calcpairs_updateBV0_gen(d,d->eps->V,DS_MAT_Q);
226: if (d->W) dvd_calcpairs_updateBV0_gen(d,d->W,DS_MAT_Z);
227: dvd_calcpairs_updateBV0_gen(d,d->AX,DS_MAT_Q);
228: if (d->BX) dvd_calcpairs_updateBV0_gen(d,d->BX,DS_MAT_Q);
229: dvd_calcpairs_updateproj(d);
230: /* Update signature */
231: if (d->nBds) {
232: VecCreateSeq(PETSC_COMM_SELF,l+d->V_tra_e,&v1);
233: BVSetActiveColumns(d->eps->V,0,l+d->V_tra_e);
234: BVGetSignature(d->eps->V,v1);
235: VecGetArray(v1,&pv);
236: for (i=0;i<d->V_tra_e;i++) pv[l+i] = d->nBds[i];
237: VecRestoreArray(v1,&pv);
238: BVSetSignature(d->eps->V,v1);
239: BVSetActiveColumns(d->eps->V,l,k);
240: VecDestroy(&v1);
241: }
242: k = l+d->V_tra_e;
243: l+= d->V_tra_s;
244: } else {
245: /* 2. V <- orth(V, V_new) */
246: dvd_orthV(d->eps->V,l+d->V_new_s,l+d->V_new_e);
247: /* 3. AV <- [AV A * V(V_new_s:V_new_e-1)] */
248: /* Check consistency */
249: PetscAssert(k-l==d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
250: for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
251: BVGetColumn(d->eps->V,i,&v1);
252: BVGetColumn(d->AX,i,&v2);
253: MatMult(d->A,v1,v2);
254: BVRestoreColumn(d->eps->V,i,&v1);
255: BVRestoreColumn(d->AX,i,&v2);
256: }
257: /* 4. BV <- [BV B * V(V_new_s:V_new_e-1)] */
258: if (d->BX) {
259: /* Check consistency */
260: PetscAssert(k-l==d->V_new_s,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Consistency broken");
261: for (i=l+d->V_new_s;i<l+d->V_new_e;i++) {
262: BVGetColumn(d->eps->V,i,&v1);
263: BVGetColumn(d->BX,i,&v2);
264: MatMult(d->B,v1,v2);
265: BVRestoreColumn(d->eps->V,i,&v1);
266: BVRestoreColumn(d->BX,i,&v2);
267: }
268: }
269: /* 5. W <- [W f(AV,BV)] */
270: if (d->W) {
271: d->calcpairs_W(d);
272: dvd_orthV(d->W,l+d->V_new_s,l+d->V_new_e);
273: }
274: /* 6. H <- W' * AX; G <- W' * BX */
275: BVSetActiveColumns(d->eps->V,l+d->V_new_s,l+d->V_new_e);
276: BVSetActiveColumns(d->AX,l+d->V_new_s,l+d->V_new_e);
277: if (d->BX) BVSetActiveColumns(d->BX,l+d->V_new_s,l+d->V_new_e);
278: if (d->W) BVSetActiveColumns(d->W,l+d->V_new_s,l+d->V_new_e);
279: BVMatProject(d->AX,NULL,d->W?d->W:d->eps->V,d->H);
280: if (d->G) BVMatProject(d->BX?d->BX:d->eps->V,NULL,d->W?d->W:d->eps->V,d->G);
281: BVSetActiveColumns(d->eps->V,l,k);
282: BVSetActiveColumns(d->AX,l,k);
283: if (d->BX) BVSetActiveColumns(d->BX,l,k);
284: if (d->W) BVSetActiveColumns(d->W,l,k);
286: /* Perform the transformation on the projected problem */
287: if (d->W) d->calcpairs_proj_trans(d);
288: k = l+d->V_new_e;
289: }
290: BVSetActiveColumns(d->eps->V,l,k);
291: BVSetActiveColumns(d->AX,l,k);
292: if (d->BX) BVSetActiveColumns(d->BX,l,k);
293: if (d->W) BVSetActiveColumns(d->W,l,k);
295: /* Solve the projected problem */
296: dvd_calcpairs_projeig_solve(d);
298: d->V_tra_s = d->V_tra_e = 0;
299: d->V_new_s = d->V_new_e;
300: return 0;
301: }
303: static PetscErrorCode dvd_calcpairs_apply_arbitrary(dvdDashboard *d,PetscInt r_s,PetscInt r_e,PetscScalar *rr,PetscScalar *ri)
304: {
305: PetscInt i,k,ld;
306: PetscScalar *pX;
307: Vec *X,xr,xi;
308: #if defined(PETSC_USE_COMPLEX)
309: PetscInt N=1;
310: #else
311: PetscInt N=2,j;
312: #endif
314: /* Quick exit without neither arbitrary selection nor harmonic extraction */
315: if (!d->eps->arbitrary && !d->calcpairs_eig_backtrans) return 0;
317: /* Quick exit without arbitrary selection, but with harmonic extraction */
318: if (d->calcpairs_eig_backtrans) {
319: for (i=r_s; i<r_e; i++) d->calcpairs_eig_backtrans(d,d->eigr[i],d->eigi[i],&rr[i-r_s],&ri[i-r_s]);
320: }
321: if (!d->eps->arbitrary) return 0;
323: SlepcVecPoolGetVecs(d->auxV,N,&X);
324: DSGetLeadingDimension(d->eps->ds,&ld);
325: for (i=r_s;i<r_e;i++) {
326: k = i;
327: DSVectors(d->eps->ds,DS_MAT_X,&k,NULL);
328: DSGetArray(d->eps->ds,DS_MAT_X,&pX);
329: dvd_improvex_compute_X(d,i,k+1,X,pX,ld);
330: DSRestoreArray(d->eps->ds,DS_MAT_X,&pX);
331: #if !defined(PETSC_USE_COMPLEX)
332: if (d->nX[i] != 1.0) {
333: for (j=i;j<k+1;j++) VecScale(X[j-i],1.0/d->nX[i]);
334: }
335: xr = X[0];
336: xi = X[1];
337: if (i == k) VecSet(xi,0.0);
338: #else
339: xr = X[0];
340: xi = NULL;
341: if (d->nX[i] != 1.0) VecScale(xr,1.0/d->nX[i]);
342: #endif
343: (d->eps->arbitrary)(rr[i-r_s],ri[i-r_s],xr,xi,&rr[i-r_s],&ri[i-r_s],d->eps->arbitraryctx);
344: #if !defined(PETSC_USE_COMPLEX)
345: if (i != k) {
346: rr[i+1-r_s] = rr[i-r_s];
347: ri[i+1-r_s] = ri[i-r_s];
348: i++;
349: }
350: #endif
351: }
352: SlepcVecPoolRestoreVecs(d->auxV,N,&X);
353: return 0;
354: }
356: static PetscErrorCode dvd_calcpairs_selectPairs(dvdDashboard *d,PetscInt n)
357: {
358: PetscInt k,lV,kV,nV;
359: PetscScalar *rr,*ri;
361: BVGetActiveColumns(d->eps->V,&lV,&kV);
362: nV = kV - lV;
363: n = PetscMin(n,nV);
364: if (n <= 0) return 0;
365: /* Put the best n pairs at the beginning. Useful for restarting */
366: if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
367: PetscMalloc1(nV,&rr);
368: PetscMalloc1(nV,&ri);
369: dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);
370: } else {
371: rr = d->eigr;
372: ri = d->eigi;
373: }
374: k = n;
375: DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);
376: /* Put the best pair at the beginning. Useful to check its residual */
377: #if !defined(PETSC_USE_COMPLEX)
378: if (n != 1 && (n != 2 || d->eigi[0] == 0.0))
379: #else
380: if (n != 1)
381: #endif
382: {
383: dvd_calcpairs_apply_arbitrary(d,0,nV,rr,ri);
384: k = 1;
385: DSSort(d->eps->ds,d->eigr,d->eigi,rr,ri,&k);
386: }
387: DSSynchronize(d->eps->ds,d->eigr,d->eigi);
389: if (d->calcpairs_eigs_trans) d->calcpairs_eigs_trans(d);
390: if (d->eps->arbitrary || d->calcpairs_eig_backtrans) {
391: PetscFree(rr);
392: PetscFree(ri);
393: }
394: return 0;
395: }
397: static PetscErrorCode EPSXDComputeDSConv(dvdDashboard *d)
398: {
399: PetscInt i,ld;
400: Vec v;
401: Mat A,B,H0,G0;
402: PetscScalar *pA;
403: const PetscScalar *pv;
404: PetscBool symm;
406: BVSetActiveColumns(d->eps->V,0,d->eps->nconv);
407: PetscObjectTypeCompare((PetscObject)d->eps->ds,DSHEP,&symm);
408: if (symm) return 0;
409: DSSetDimensions(d->eps->ds,d->eps->nconv,0,0);
410: DSGetMat(d->eps->ds,DS_MAT_A,&A);
411: MatDenseGetSubMatrix(d->H,0,d->eps->nconv,0,d->eps->nconv,&H0);
412: MatCopy(H0,A,SAME_NONZERO_PATTERN);
413: MatDenseRestoreSubMatrix(d->H,&H0);
414: DSRestoreMat(d->eps->ds,DS_MAT_A,&A);
415: if (d->G) {
416: DSGetMat(d->eps->ds,DS_MAT_B,&B);
417: MatDenseGetSubMatrix(d->G,0,d->eps->nconv,0,d->eps->nconv,&G0);
418: MatCopy(G0,B,SAME_NONZERO_PATTERN);
419: MatDenseRestoreSubMatrix(d->G,&G0);
420: DSRestoreMat(d->eps->ds,DS_MAT_B,&B);
421: }
422: /* Set the signature on projected matrix B */
423: if (DVD_IS(d->sEP,DVD_EP_INDEFINITE)) {
424: DSGetLeadingDimension(d->eps->ds,&ld);
425: DSGetArray(d->eps->ds,DS_MAT_B,&pA);
426: PetscArrayzero(pA,d->eps->nconv*ld);
427: VecCreateSeq(PETSC_COMM_SELF,d->eps->nconv,&v);
428: BVGetSignature(d->eps->V,v);
429: VecGetArrayRead(v,&pv);
430: for (i=0;i<d->eps->nconv;i++) pA[i+ld*i] = pv[i];
431: VecRestoreArrayRead(v,&pv);
432: VecDestroy(&v);
433: DSRestoreArray(d->eps->ds,DS_MAT_B,&pA);
434: }
435: DSSetState(d->eps->ds,DS_STATE_RAW);
436: DSSolve(d->eps->ds,d->eps->eigr,d->eps->eigi);
437: DSSynchronize(d->eps->ds,d->eps->eigr,d->eps->eigi);
438: if (d->W) {
439: for (i=0;i<d->eps->nconv;i++) d->calcpairs_eig_backtrans(d,d->eps->eigr[i],d->eps->eigi[i],&d->eps->eigr[i],&d->eps->eigi[i]);
440: }
441: return 0;
442: }
444: /*
445: Compute the residual vectors R(i) <- (AV - BV*eigr(i))*pX(i), and also
446: the norm associated to the Schur pair, where i = r_s..r_e
447: */
448: static PetscErrorCode dvd_calcpairs_res_0(dvdDashboard *d,PetscInt r_s,PetscInt r_e)
449: {
450: PetscInt i,ldpX;
451: PetscScalar *pX;
452: BV BX = d->BX?d->BX:d->eps->V;
453: Vec *R;
455: DSGetLeadingDimension(d->eps->ds,&ldpX);
456: DSGetArray(d->eps->ds,DS_MAT_Q,&pX);
457: /* nX(i) <- ||X(i)|| */
458: dvd_improvex_compute_X(d,r_s,r_e,NULL,pX,ldpX);
459: SlepcVecPoolGetVecs(d->auxV,r_e-r_s,&R);
460: for (i=r_s;i<r_e;i++) {
461: /* R(i-r_s) <- AV*pX(i) */
462: BVMultVec(d->AX,1.0,0.0,R[i-r_s],&pX[ldpX*i]);
463: /* R(i-r_s) <- R(i-r_s) - eigr(i)*BX*pX(i) */
464: BVMultVec(BX,-d->eigr[i],1.0,R[i-r_s],&pX[ldpX*i]);
465: }
466: DSRestoreArray(d->eps->ds,DS_MAT_Q,&pX);
467: d->calcpairs_proj_res(d,r_s,r_e,R);
468: SlepcVecPoolRestoreVecs(d->auxV,r_e-r_s,&R);
469: return 0;
470: }
472: static PetscErrorCode dvd_calcpairs_proj_res(dvdDashboard *d,PetscInt r_s,PetscInt r_e,Vec *R)
473: {
474: PetscInt i,l,k;
475: PetscBool lindep=PETSC_FALSE;
476: BV cX;
478: if (d->W) cX = d->W; /* If left subspace exists, R <- orth(cY, R), nR[i] <- ||R[i]|| */
479: else if (!(DVD_IS(d->sEP, DVD_EP_STD) && DVD_IS(d->sEP, DVD_EP_HERMITIAN))) cX = d->eps->V; /* If not HEP, R <- orth(cX, R), nR[i] <- ||R[i]|| */
480: else cX = NULL; /* Otherwise, nR[i] <- ||R[i]|| */
482: if (cX) {
483: BVGetActiveColumns(cX,&l,&k);
484: BVSetActiveColumns(cX,0,l);
485: for (i=0;i<r_e-r_s;i++) BVOrthogonalizeVec(cX,R[i],NULL,&d->nR[r_s+i],&lindep);
486: BVSetActiveColumns(cX,l,k);
487: if (lindep || (PetscAbs(d->nR[r_s+i]) < PETSC_MACHINE_EPSILON)) PetscInfo(d->eps,"The computed eigenvector residual %" PetscInt_FMT " is too low, %g!\n",r_s+i,(double)(d->nR[r_s+i]));
488: } else {
489: for (i=0;i<r_e-r_s;i++) VecNormBegin(R[i],NORM_2,&d->nR[r_s+i]);
490: for (i=0;i<r_e-r_s;i++) VecNormEnd(R[i],NORM_2,&d->nR[r_s+i]);
491: }
492: return 0;
493: }
495: PetscErrorCode dvd_calcpairs_qz(dvdDashboard *d,dvdBlackboard *b,PetscBool borth,PetscBool harm)
496: {
497: PetscBool std_probl,her_probl,ind_probl;
498: DSType dstype;
499: Vec v1;
501: std_probl = DVD_IS(d->sEP,DVD_EP_STD)? PETSC_TRUE: PETSC_FALSE;
502: her_probl = DVD_IS(d->sEP,DVD_EP_HERMITIAN)? PETSC_TRUE: PETSC_FALSE;
503: ind_probl = DVD_IS(d->sEP,DVD_EP_INDEFINITE)? PETSC_TRUE: PETSC_FALSE;
505: /* Setting configuration constrains */
506: b->max_size_proj = PetscMax(b->max_size_proj,b->max_size_V);
507: d->W_shift = d->B? PETSC_TRUE: PETSC_FALSE;
509: /* Setup the step */
510: if (b->state >= DVD_STATE_CONF) {
511: d->max_size_P = b->max_size_P;
512: d->max_size_proj = b->max_size_proj;
513: /* Create a DS if the method works with Schur decompositions */
514: d->calcPairs = dvd_calcpairs_proj;
515: d->calcpairs_residual = dvd_calcpairs_res_0;
516: d->calcpairs_proj_res = dvd_calcpairs_proj_res;
517: d->calcpairs_selectPairs = dvd_calcpairs_selectPairs;
518: /* Create and configure a DS for solving the projected problems */
519: if (d->W) dstype = DSGNHEP; /* If we use harmonics */
520: else {
521: if (ind_probl) dstype = DSGHIEP;
522: else if (std_probl) dstype = her_probl? DSHEP : DSNHEP;
523: else dstype = her_probl? DSGHEP : DSGNHEP;
524: }
525: DSSetType(d->eps->ds,dstype);
526: DSAllocate(d->eps->ds,d->eps->ncv);
527: /* Create various vector basis */
528: if (harm) {
529: BVDuplicateResize(d->eps->V,d->eps->ncv,&d->W);
530: BVSetMatrix(d->W,NULL,PETSC_FALSE);
531: } else d->W = NULL;
532: BVDuplicateResize(d->eps->V,d->eps->ncv,&d->AX);
533: BVSetMatrix(d->AX,NULL,PETSC_FALSE);
534: BVDuplicateResize(d->eps->V,d->eps->ncv,&d->auxBV);
535: BVSetMatrix(d->auxBV,NULL,PETSC_FALSE);
536: if (d->B) {
537: BVDuplicateResize(d->eps->V,d->eps->ncv,&d->BX);
538: BVSetMatrix(d->BX,NULL,PETSC_FALSE);
539: } else d->BX = NULL;
540: MatCreateVecsEmpty(d->A,&v1,NULL);
541: SlepcVecPoolCreate(v1,0,&d->auxV);
542: VecDestroy(&v1);
543: /* Create projected problem matrices */
544: MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->H);
545: if (!std_probl) MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->G);
546: else d->G = NULL;
547: if (her_probl) {
548: MatSetOption(d->H,MAT_HERMITIAN,PETSC_TRUE);
549: if (d->G) MatSetOption(d->G,MAT_HERMITIAN,PETSC_TRUE);
550: }
552: if (ind_probl) PetscMalloc1(d->eps->ncv,&d->nBds);
553: else d->nBds = NULL;
554: MatCreateSeqDense(PETSC_COMM_SELF,d->eps->ncv,d->eps->ncv,NULL,&d->auxM);
556: EPSDavidsonFLAdd(&d->startList,dvd_calcpairs_qz_start);
557: EPSDavidsonFLAdd(&d->endList,EPSXDComputeDSConv);
558: EPSDavidsonFLAdd(&d->destroyList,dvd_calcpairs_qz_d);
559: }
560: return 0;
561: }