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