Actual source code: dsghiep.c
slepc-main 2025-03-25
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: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: static PetscErrorCode DSAllocate_GHIEP(DS ds,PetscInt ld)
15: {
16: PetscFunctionBegin;
17: PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
18: PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
19: PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
20: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
21: PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
22: PetscCall(PetscFree(ds->perm));
23: PetscCall(PetscMalloc1(ld,&ds->perm));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: PetscErrorCode DSSwitchFormat_GHIEP(DS ds,PetscBool tocompact)
28: {
29: PetscReal *T,*S;
30: PetscScalar *A,*B;
31: PetscInt i,n,ld;
33: PetscFunctionBegin;
34: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
35: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
36: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
37: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
38: n = ds->n;
39: ld = ds->ld;
40: if (tocompact) { /* switch from dense (arrow) to compact storage */
41: PetscCall(PetscArrayzero(T,n));
42: PetscCall(PetscArrayzero(T+ld,n));
43: PetscCall(PetscArrayzero(T+2*ld,n));
44: PetscCall(PetscArrayzero(S,n));
45: for (i=0;i<n-1;i++) {
46: T[i] = PetscRealPart(A[i+i*ld]);
47: T[ld+i] = PetscRealPart(A[i+1+i*ld]);
48: S[i] = PetscRealPart(B[i+i*ld]);
49: }
50: T[n-1] = PetscRealPart(A[n-1+(n-1)*ld]);
51: S[n-1] = PetscRealPart(B[n-1+(n-1)*ld]);
52: for (i=ds->l;i<ds->k;i++) T[2*ld+i] = PetscRealPart(A[ds->k+i*ld]);
53: } else { /* switch from compact (arrow) to dense storage */
54: for (i=0;i<n;i++) {
55: PetscCall(PetscArrayzero(A+i*ld,n));
56: PetscCall(PetscArrayzero(B+i*ld,n));
57: }
58: for (i=0;i<n-1;i++) {
59: A[i+i*ld] = T[i];
60: A[i+1+i*ld] = T[ld+i];
61: A[i+(i+1)*ld] = T[ld+i];
62: B[i+i*ld] = S[i];
63: }
64: A[n-1+(n-1)*ld] = T[n-1];
65: B[n-1+(n-1)*ld] = S[n-1];
66: for (i=ds->l;i<ds->k;i++) {
67: A[ds->k+i*ld] = T[2*ld+i];
68: A[i+ds->k*ld] = T[2*ld+i];
69: }
70: }
71: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
72: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
73: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
74: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode DSView_GHIEP(DS ds,PetscViewer viewer)
79: {
80: PetscViewerFormat format;
81: PetscInt i,j;
82: PetscReal *T,*S,value;
83: const char *methodname[] = {
84: "QR + Inverse Iteration",
85: "HZ method",
86: "QR"
87: };
88: const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
90: PetscFunctionBegin;
91: PetscCall(PetscViewerGetFormat(viewer,&format));
92: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
93: if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
96: if (ds->compact) {
97: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
98: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
99: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
100: if (format == PETSC_VIEWER_ASCII_MATLAB) {
101: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
102: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
103: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
104: for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
105: for (i=0;i<ds->n-1;i++) {
106: if (T[i+ds->ld] !=0 && i!=ds->k-1) {
107: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+2,i+1,(double)T[i+ds->ld]));
108: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+2,(double)T[i+ds->ld]));
109: }
110: }
111: for (i = ds->l;i<ds->k;i++) {
112: if (T[i+2*ds->ld]) {
113: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",ds->k+1,i+1,(double)T[i+2*ds->ld]));
114: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,ds->k+1,(double)T[i+2*ds->ld]));
115: }
116: }
117: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_A]));
119: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
120: PetscCall(PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
121: PetscCall(PetscViewerASCIIPrintf(viewer,"omega = [\n"));
122: for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
123: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_B]));
125: } else {
126: PetscCall(PetscViewerASCIIPrintf(viewer,"T\n"));
127: for (i=0;i<ds->n;i++) {
128: for (j=0;j<ds->n;j++) {
129: if (i==j) value = T[i];
130: else if (i==j+1 || j==i+1) value = T[PetscMin(i,j)+ds->ld];
131: else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = T[PetscMin(i,j)+2*ds->ld];
132: else value = 0.0;
133: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
134: }
135: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
136: }
137: PetscCall(PetscViewerASCIIPrintf(viewer,"omega\n"));
138: for (i=0;i<ds->n;i++) {
139: for (j=0;j<ds->n;j++) {
140: if (i==j) value = S[i];
141: else value = 0.0;
142: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
143: }
144: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
145: }
146: }
147: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
148: PetscCall(PetscViewerFlush(viewer));
149: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
150: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
151: } else {
152: PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
153: PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
154: }
155: if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode DSVectors_GHIEP_Eigen_Some(DS ds,PetscInt *idx,PetscReal *rnorm)
160: {
161: PetscReal b[4],M[4],*T,*S,d1,d2,s1,s2,e,scal1,scal2,wr1,wr2,wi,ep,norm;
162: PetscScalar *X,Y[4],alpha,szero=0.0;
163: const PetscScalar *A,*B,*Q;
164: PetscInt k;
165: PetscBLASInt two=2,n_,ld,one=1;
166: #if !defined(PETSC_USE_COMPLEX)
167: PetscBLASInt four=4;
168: #endif
170: PetscFunctionBegin;
171: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
172: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
173: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
174: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
175: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
176: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
177: k = *idx;
178: PetscCall(PetscBLASIntCast(ds->n,&n_));
179: PetscCall(PetscBLASIntCast(ds->ld,&ld));
180: if (k < ds->n-1) e = (ds->compact)?T[k+ld]:PetscRealPart(A[(k+1)+ld*k]);
181: else e = 0.0;
182: if (e == 0.0) { /* Real */
183: if (ds->state>=DS_STATE_CONDENSED) PetscCall(PetscArraycpy(X+k*ld,Q+k*ld,ld));
184: else {
185: PetscCall(PetscArrayzero(X+k*ds->ld,ds->ld));
186: X[k+k*ds->ld] = 1.0;
187: }
188: if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
189: } else { /* 2x2 block */
190: if (ds->compact) {
191: s1 = S[k];
192: d1 = T[k];
193: s2 = S[k+1];
194: d2 = T[k+1];
195: } else {
196: s1 = PetscRealPart(B[k*ld+k]);
197: d1 = PetscRealPart(A[k+k*ld]);
198: s2 = PetscRealPart(B[(k+1)*ld+k+1]);
199: d2 = PetscRealPart(A[k+1+(k+1)*ld]);
200: }
201: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
202: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
203: ep = LAPACKlamch_("S");
204: /* Compute eigenvalues of the block */
205: PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
206: PetscCheck(wi!=0.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Real block in DSVectors_GHIEP");
207: /* Complex eigenvalues */
208: PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
209: wr1 /= scal1;
210: wi /= scal1;
211: #if !defined(PETSC_USE_COMPLEX)
212: if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
213: Y[0] = wr1-s2*d2; Y[1] = s2*e; Y[2] = wi; Y[3] = 0.0;
214: } else {
215: Y[0] = s1*e; Y[1] = wr1-s1*d1; Y[2] = 0.0; Y[3] = wi;
216: }
217: norm = BLASnrm2_(&four,Y,&one);
218: norm = 1.0/norm;
219: if (ds->state >= DS_STATE_CONDENSED) {
220: alpha = norm;
221: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&two,&two,&alpha,Q+k*ld,&ld,Y,&two,&szero,X+k*ld,&ld));
222: if (rnorm) *rnorm = SlepcAbsEigenvalue(X[ds->n-1+k*ld],X[ds->n-1+(k+1)*ld]);
223: } else {
224: PetscCall(PetscArrayzero(X+k*ld,2*ld));
225: X[k*ld+k] = Y[0]*norm;
226: X[k*ld+k+1] = Y[1]*norm;
227: X[(k+1)*ld+k] = Y[2]*norm;
228: X[(k+1)*ld+k+1] = Y[3]*norm;
229: }
230: #else
231: if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
232: Y[0] = PetscCMPLX(wr1-s2*d2,wi);
233: Y[1] = s2*e;
234: } else {
235: Y[0] = s1*e;
236: Y[1] = PetscCMPLX(wr1-s1*d1,wi);
237: }
238: norm = BLASnrm2_(&two,Y,&one);
239: norm = 1.0/norm;
240: if (ds->state >= DS_STATE_CONDENSED) {
241: alpha = norm;
242: PetscCallBLAS("BLASgemv",BLASgemv_("N",&n_,&two,&alpha,Q+k*ld,&ld,Y,&one,&szero,X+k*ld,&one));
243: if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
244: } else {
245: PetscCall(PetscArrayzero(X+k*ld,2*ld));
246: X[k*ld+k] = Y[0]*norm;
247: X[k*ld+k+1] = Y[1]*norm;
248: }
249: X[(k+1)*ld+k] = PetscConj(X[k*ld+k]);
250: X[(k+1)*ld+k+1] = PetscConj(X[k*ld+k+1]);
251: #endif
252: (*idx)++;
253: }
254: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
255: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
256: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
257: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
258: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
259: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
264: {
265: PetscScalar *Z;
266: const PetscScalar *A,*Q;
267: PetscInt i;
268: PetscReal e,*T,*S;
270: PetscFunctionBegin;
271: switch (mat) {
272: case DS_MAT_X:
273: case DS_MAT_Y:
274: if (k) PetscCall(DSVectors_GHIEP_Eigen_Some(ds,k,rnorm));
275: else {
276: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
277: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
278: PetscCall(MatDenseGetArray(ds->omat[mat],&Z));
279: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
280: for (i=0; i<ds->n; i++) {
281: e = (ds->compact)?T[i+ds->ld]:PetscRealPart(A[(i+1)+ds->ld*i]);
282: if (e == 0.0) { /* real */
283: if (ds->state >= DS_STATE_CONDENSED) PetscCall(PetscArraycpy(Z+i*ds->ld,Q+i*ds->ld,ds->ld));
284: else {
285: PetscCall(PetscArrayzero(Z+i*ds->ld,ds->ld));
286: Z[i+i*ds->ld] = 1.0;
287: }
288: } else {
289: PetscCall(DSVectors_GHIEP_Eigen_Some(ds,&i,rnorm));
290: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
291: S[i]=S[i-1]=1.0;
292: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
293: }
294: }
295: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
296: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
297: PetscCall(MatDenseRestoreArray(ds->omat[mat],&Z));
298: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
299: }
300: break;
301: case DS_MAT_U:
302: case DS_MAT_V:
303: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
304: default:
305: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
306: }
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*
311: Extract the eigenvalues contained in the block-diagonal of the indefinite problem.
312: Only the index range n0..n1 is processed.
313: */
314: PetscErrorCode DSGHIEPComplexEigs(DS ds,PetscInt n0,PetscInt n1,PetscScalar *wr,PetscScalar *wi)
315: {
316: PetscInt k,ld;
317: PetscBLASInt two=2;
318: const PetscScalar *A,*B;
319: PetscReal *D,*T,b[4],M[4],d1,d2,s1,s2,e,scal1,scal2,ep,wr1,wr2,wi1;
321: PetscFunctionBegin;
322: ld = ds->ld;
323: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
324: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
325: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
326: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
327: for (k=n0;k<n1;k++) {
328: if (k < n1-1) e = (ds->compact)?T[ld+k]:PetscRealPart(A[(k+1)+ld*k]);
329: else e = 0.0;
330: if (e==0.0) { /* real eigenvalue */
331: wr[k] = (ds->compact)?T[k]/D[k]:A[k+k*ld]/B[k+k*ld];
332: #if !defined(PETSC_USE_COMPLEX)
333: wi[k] = 0.0 ;
334: #endif
335: } else { /* diagonal block */
336: if (ds->compact) {
337: s1 = D[k];
338: d1 = T[k];
339: s2 = D[k+1];
340: d2 = T[k+1];
341: } else {
342: s1 = PetscRealPart(B[k*ld+k]);
343: d1 = PetscRealPart(A[k+k*ld]);
344: s2 = PetscRealPart(B[(k+1)*ld+k+1]);
345: d2 = PetscRealPart(A[k+1+(k+1)*ld]);
346: }
347: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
348: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
349: ep = LAPACKlamch_("S");
350: /* Compute eigenvalues of the block */
351: PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi1));
352: PetscCheck(scal1>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
353: if (wi1==0.0) { /* Real eigenvalues */
354: PetscCheck(scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
355: wr[k] = wr1/scal1; wr[k+1] = wr2/scal2;
356: #if !defined(PETSC_USE_COMPLEX)
357: wi[k] = wi[k+1] = 0.0;
358: #endif
359: } else { /* Complex eigenvalues */
360: #if !defined(PETSC_USE_COMPLEX)
361: wr[k] = wr1/scal1;
362: wr[k+1] = wr[k];
363: wi[k] = wi1/scal1;
364: wi[k+1] = -wi[k];
365: #else
366: wr[k] = PetscCMPLX(wr1,wi1)/scal1;
367: wr[k+1] = PetscConj(wr[k]);
368: #endif
369: }
370: k++;
371: }
372: }
373: #if defined(PETSC_USE_COMPLEX)
374: if (wi) {
375: for (k=n0;k<n1;k++) wi[k] = 0.0;
376: }
377: #endif
378: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
379: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
380: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
381: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: static PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
386: {
387: PetscInt n,i,*perm;
388: PetscReal *d,*e,*s;
390: PetscFunctionBegin;
391: #if !defined(PETSC_USE_COMPLEX)
392: PetscAssertPointer(wi,3);
393: #endif
394: n = ds->n;
395: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
396: e = d + ds->ld;
397: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
398: PetscCall(DSAllocateWork_Private(ds,ds->ld,ds->ld,0));
399: perm = ds->perm;
400: if (!rr) {
401: rr = wr;
402: ri = wi;
403: }
404: PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_TRUE));
405: if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_TRUE));
406: PetscCall(PetscArraycpy(ds->work,wr,n));
407: for (i=ds->l;i<n;i++) wr[i] = *(ds->work+perm[i]);
408: #if !defined(PETSC_USE_COMPLEX)
409: PetscCall(PetscArraycpy(ds->work,wi,n));
410: for (i=ds->l;i<n;i++) wi[i] = *(ds->work+perm[i]);
411: #endif
412: PetscCall(PetscArraycpy(ds->rwork,s,n));
413: for (i=ds->l;i<n;i++) s[i] = *(ds->rwork+perm[i]);
414: PetscCall(PetscArraycpy(ds->rwork,d,n));
415: for (i=ds->l;i<n;i++) d[i] = *(ds->rwork+perm[i]);
416: PetscCall(PetscArraycpy(ds->rwork,e,n-1));
417: PetscCall(PetscArrayzero(e+ds->l,n-1-ds->l));
418: for (i=ds->l;i<n-1;i++) {
419: if (perm[i]<n-1) e[i] = *(ds->rwork+perm[i]);
420: }
421: if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
422: PetscCall(DSPermuteColumns_Private(ds,ds->l,n,n,DS_MAT_Q,perm));
423: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
424: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode DSUpdateExtraRow_GHIEP(DS ds)
429: {
430: PetscInt i;
431: PetscBLASInt n,ld,incx=1;
432: PetscScalar *A,*x,*y,one=1.0,zero=0.0;
433: const PetscScalar *Q;
434: PetscReal *T,*b,*r,beta;
436: PetscFunctionBegin;
437: PetscCall(PetscBLASIntCast(ds->n,&n));
438: PetscCall(PetscBLASIntCast(ds->ld,&ld));
439: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
440: if (ds->compact) {
441: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
442: b = T+ld;
443: r = T+2*ld;
444: beta = b[n-1]; /* in compact, we assume all entries are zero except the last one */
445: for (i=0;i<n;i++) r[i] = PetscRealPart(beta*Q[n-1+i*ld]);
446: ds->k = n;
447: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
448: } else {
449: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
450: PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
451: x = ds->work;
452: y = ds->work+ld;
453: for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
454: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
455: for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
456: ds->k = n;
457: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
458: }
459: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*
464: Get eigenvectors with inverse iteration.
465: The system matrix is in Hessenberg form.
466: */
467: PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
468: {
469: PetscInt i,off;
470: PetscBLASInt *select,*infoC,ld,n1,mout,info;
471: const PetscScalar *A,*B;
472: PetscScalar *H,*X;
473: PetscReal *s,*d,*e;
474: #if defined(PETSC_USE_COMPLEX)
475: PetscInt j;
476: #endif
478: PetscFunctionBegin;
479: PetscCall(PetscBLASIntCast(ds->ld,&ld));
480: PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
481: PetscCall(DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld));
482: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
483: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
484: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
485: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
486: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
487: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
488: e = d + ld;
489: select = ds->iwork;
490: infoC = ds->iwork + ld;
491: off = ds->l+ds->l*ld;
492: if (ds->compact) {
493: H[off] = d[ds->l]*s[ds->l];
494: H[off+ld] = e[ds->l]*s[ds->l];
495: for (i=ds->l+1;i<ds->n-1;i++) {
496: H[i+(i-1)*ld] = e[i-1]*s[i];
497: H[i+i*ld] = d[i]*s[i];
498: H[i+(i+1)*ld] = e[i]*s[i];
499: }
500: H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
501: H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
502: } else {
503: s[ds->l] = PetscRealPart(B[off]);
504: H[off] = A[off]*s[ds->l];
505: H[off+ld] = A[off+ld]*s[ds->l];
506: for (i=ds->l+1;i<ds->n-1;i++) {
507: s[i] = PetscRealPart(B[i+i*ld]);
508: H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
509: H[i+i*ld] = A[i+i*ld]*s[i];
510: H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
511: }
512: s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
513: H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
514: H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
515: }
516: PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
517: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
518: for (i=0;i<n1;i++) select[i] = 1;
519: #if !defined(PETSC_USE_COMPLEX)
520: PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,NULL,infoC,&info));
521: #else
522: PetscCallBLAS("LAPACKhsein",LAPACKhsein_("R","N","N",select,&n1,H+off,&ld,wr+ds->l,NULL,&ld,X+off,&ld,&n1,&mout,ds->work,ds->rwork,NULL,infoC,&info));
524: /* Separate real and imaginary part of complex eigenvectors */
525: for (j=ds->l;j<ds->n;j++) {
526: if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
527: for (i=ds->l;i<ds->n;i++) {
528: X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
529: X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
530: }
531: j++;
532: }
533: }
534: #endif
535: SlepcCheckLapackInfo("hsein",info);
536: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
537: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
538: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
539: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
540: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
541: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
542: PetscCall(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE));
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*
547: Undo 2x2 blocks that have real eigenvalues.
548: */
549: PetscErrorCode DSGHIEPRealBlocks(DS ds)
550: {
551: PetscInt i;
552: PetscReal e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
553: PetscReal maxy,ep,scal1,scal2,snorm;
554: PetscReal *T,*D,b[4],M[4],wr1,wr2,wi;
555: PetscScalar *A,*B,*Q,Y[4],sone=1.0,szero=0.0;
556: PetscBLASInt m,two=2,ld;
557: PetscBool isreal;
559: PetscFunctionBegin;
560: PetscCall(PetscBLASIntCast(ds->ld,&ld));
561: PetscCall(PetscBLASIntCast(ds->n-ds->l,&m));
562: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
563: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
564: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
565: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
566: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
567: PetscCall(DSAllocateWork_Private(ds,2*m,0,0));
568: for (i=ds->l;i<ds->n-1;i++) {
569: e = (ds->compact)?T[ld+i]:PetscRealPart(A[(i+1)+ld*i]);
570: if (e != 0.0) { /* 2x2 block */
571: if (ds->compact) {
572: s1 = D[i];
573: d1 = T[i];
574: s2 = D[i+1];
575: d2 = T[i+1];
576: } else {
577: s1 = PetscRealPart(B[i*ld+i]);
578: d1 = PetscRealPart(A[i*ld+i]);
579: s2 = PetscRealPart(B[(i+1)*ld+i+1]);
580: d2 = PetscRealPart(A[(i+1)*ld+i+1]);
581: }
582: isreal = PETSC_FALSE;
583: if (s1==s2) { /* apply a Jacobi rotation to compute the eigendecomposition */
584: dd = d1-d2;
585: if (2*PetscAbsReal(e) <= dd) {
586: t = 2*e/dd;
587: t = t/(1 + PetscSqrtReal(1+t*t));
588: } else {
589: t = dd/(2*e);
590: ss = (t>=0)?1.0:-1.0;
591: t = ss/(PetscAbsReal(t)+PetscSqrtReal(1+t*t));
592: }
593: Y[0] = 1/PetscSqrtReal(1 + t*t); Y[3] = Y[0]; /* c */
594: Y[1] = Y[0]*t; Y[2] = -Y[1]; /* s */
595: wr1 = d1+t*e; wr2 = d2-t*e;
596: ss1 = s1; ss2 = s2;
597: isreal = PETSC_TRUE;
598: } else {
599: ss1 = 1.0; ss2 = 1.0,
600: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
601: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
602: ep = LAPACKlamch_("S");
604: /* Compute eigenvalues of the block */
605: PetscCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
606: if (wi==0.0) { /* Real eigenvalues */
607: isreal = PETSC_TRUE;
608: PetscCheck(scal1>=ep && scal2>=ep,PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
609: wr1 /= scal1;
610: wr2 /= scal2;
611: if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
612: Y[0] = wr1-s2*d2;
613: Y[1] = s2*e;
614: } else {
615: Y[0] = s1*e;
616: Y[1] = wr1-s1*d1;
617: }
618: /* normalize with a signature*/
619: maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
620: scal1 = PetscRealPart(Y[0])/maxy;
621: scal2 = PetscRealPart(Y[1])/maxy;
622: snorm = scal1*scal1*s1 + scal2*scal2*s2;
623: if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
624: snorm = maxy*PetscSqrtReal(snorm);
625: Y[0] = Y[0]/snorm;
626: Y[1] = Y[1]/snorm;
627: if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
628: Y[2] = wr2-s2*d2;
629: Y[3] = s2*e;
630: } else {
631: Y[2] = s1*e;
632: Y[3] = wr2-s1*d1;
633: }
634: maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
635: scal1 = PetscRealPart(Y[2])/maxy;
636: scal2 = PetscRealPart(Y[3])/maxy;
637: snorm = scal1*scal1*s1 + scal2*scal2*s2;
638: if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
639: snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
640: }
641: wr1 *= ss1; wr2 *= ss2;
642: }
643: if (isreal) {
644: if (ds->compact) {
645: D[i] = ss1;
646: T[i] = wr1;
647: D[i+1] = ss2;
648: T[i+1] = wr2;
649: T[ld+i] = 0.0;
650: } else {
651: B[i*ld+i] = ss1;
652: A[i*ld+i] = wr1;
653: B[(i+1)*ld+i+1] = ss2;
654: A[(i+1)*ld+i+1] = wr2;
655: A[(i+1)+ld*i] = 0.0;
656: A[i+ld*(i+1)] = 0.0;
657: }
658: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&two,&two,&sone,Q+ds->l+i*ld,&ld,Y,&two,&szero,ds->work,&m));
659: PetscCall(PetscArraycpy(Q+ds->l+i*ld,ds->work,m));
660: PetscCall(PetscArraycpy(Q+ds->l+(i+1)*ld,ds->work+m,m));
661: }
662: i++;
663: }
664: }
665: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
666: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
667: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
668: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
669: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: static PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
674: {
675: PetscInt i,off;
676: PetscBLASInt n1,ld,one=1,info,lwork;
677: const PetscScalar *A,*B;
678: PetscScalar *H,*Q;
679: PetscReal *d,*e,*s;
680: #if defined(PETSC_USE_COMPLEX)
681: PetscInt j;
682: #endif
684: PetscFunctionBegin;
685: #if !defined(PETSC_USE_COMPLEX)
686: PetscAssertPointer(wi,3);
687: #endif
688: PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
689: PetscCall(PetscBLASIntCast(ds->ld,&ld));
690: off = ds->l + ds->l*ld;
691: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
692: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
693: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
694: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
695: e = d + ld;
696: #if defined(PETSC_USE_DEBUG)
697: /* Check signature */
698: for (i=0;i<ds->n;i++) {
699: PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
700: PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
701: }
702: #endif
704: /* Quick return if possible */
705: if (n1 == 1) {
706: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
707: for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
708: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
709: PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
710: if (!ds->compact) {
711: d[ds->l] = PetscRealPart(A[off]);
712: s[ds->l] = PetscRealPart(B[off]);
713: }
714: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
715: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
716: wr[ds->l] = d[ds->l]/s[ds->l];
717: if (wi) wi[ds->l] = 0.0;
718: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
719: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: PetscCall(DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2));
724: lwork = ld*ld;
726: /* Reduce to pseudotriadiagonal form */
727: PetscCall(DSIntermediate_GHIEP(ds));
729: /* Compute Eigenvalues (QR) */
730: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
731: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&H));
732: if (ds->compact) {
733: H[off] = d[ds->l]*s[ds->l];
734: H[off+ld] = e[ds->l]*s[ds->l];
735: for (i=ds->l+1;i<ds->n-1;i++) {
736: H[i+(i-1)*ld] = e[i-1]*s[i];
737: H[i+i*ld] = d[i]*s[i];
738: H[i+(i+1)*ld] = e[i]*s[i];
739: }
740: H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
741: H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
742: } else {
743: s[ds->l] = PetscRealPart(B[off]);
744: H[off] = A[off]*s[ds->l];
745: H[off+ld] = A[off+ld]*s[ds->l];
746: for (i=ds->l+1;i<ds->n-1;i++) {
747: s[i] = PetscRealPart(B[i+i*ld]);
748: H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
749: H[i+i*ld] = A[i+i*ld]*s[i];
750: H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
751: }
752: s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
753: H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
754: H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
755: }
757: #if !defined(PETSC_USE_COMPLEX)
758: PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
759: #else
760: PetscCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
761: for (i=ds->l;i<ds->n;i++) if (PetscAbsReal(PetscImaginaryPart(wr[i]))<10*PETSC_MACHINE_EPSILON) wr[i] = PetscRealPart(wr[i]);
762: /* Sort to have consecutive conjugate pairs */
763: for (i=ds->l;i<ds->n;i++) {
764: j=i+1;
765: while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
766: if (j==ds->n) {
767: PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
768: wr[i] = PetscRealPart(wr[i]);
769: } else { /* complex eigenvalue */
770: wr[j] = wr[i+1];
771: if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
772: wr[i+1] = PetscConj(wr[i]);
773: i++;
774: }
775: }
776: #endif
777: SlepcCheckLapackInfo("hseqr",info);
778: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&H));
779: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
780: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
781: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
782: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
784: /* Compute Eigenvectors with Inverse Iteration */
785: PetscCall(DSGHIEPInverseIteration(ds,wr,wi));
787: /* Recover eigenvalues from diagonal */
788: PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
789: #if defined(PETSC_USE_COMPLEX)
790: if (wi) {
791: for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
792: }
793: #endif
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: static PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
798: {
799: PetscInt i,j,off,nwu=0,n,lw,lwr,nwru=0;
800: PetscBLASInt n_,ld,info,lwork,ilo,ihi;
801: const PetscScalar *A,*B;
802: PetscScalar *H,*Q,*X;
803: PetscReal *d,*s,*scale,nrm,*rcde,*rcdv;
804: #if defined(PETSC_USE_COMPLEX)
805: PetscInt k;
806: #endif
808: PetscFunctionBegin;
809: #if !defined(PETSC_USE_COMPLEX)
810: PetscAssertPointer(wi,3);
811: #endif
812: n = ds->n-ds->l;
813: PetscCall(PetscBLASIntCast(n,&n_));
814: PetscCall(PetscBLASIntCast(ds->ld,&ld));
815: off = ds->l + ds->l*ld;
816: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
817: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
818: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
819: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
820: #if defined(PETSC_USE_DEBUG)
821: /* Check signature */
822: for (i=0;i<ds->n;i++) {
823: PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
824: PetscCheck(de==1.0 || de==-1.0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diagonal elements of the signature matrix must be 1 or -1");
825: }
826: #endif
828: /* Quick return if possible */
829: if (n_ == 1) {
830: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
831: for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
832: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
833: PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
834: if (!ds->compact) {
835: d[ds->l] = PetscRealPart(A[off]);
836: s[ds->l] = PetscRealPart(B[off]);
837: }
838: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
839: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
840: wr[ds->l] = d[ds->l]/s[ds->l];
841: if (wi) wi[ds->l] = 0.0;
842: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
843: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: lw = 14*ld+ld*ld;
848: lwr = 7*ld;
849: PetscCall(DSAllocateWork_Private(ds,lw,lwr,0));
850: scale = ds->rwork+nwru;
851: nwru += ld;
852: rcde = ds->rwork+nwru;
853: nwru += ld;
854: rcdv = ds->rwork+nwru;
856: /* Form pseudo-symmetric matrix */
857: H = ds->work+nwu;
858: nwu += n*n;
859: PetscCall(PetscArrayzero(H,n*n));
860: if (ds->compact) {
861: for (i=0;i<n-1;i++) {
862: H[i+i*n] = s[ds->l+i]*d[ds->l+i];
863: H[i+1+i*n] = s[ds->l+i+1]*d[ld+ds->l+i];
864: H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
865: }
866: H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
867: for (i=0;i<ds->k-ds->l;i++) {
868: H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
869: H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
870: }
871: } else {
872: for (j=0;j<n;j++) {
873: for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
874: }
875: }
877: /* Compute eigenpairs */
878: PetscCall(PetscBLASIntCast(lw-nwu,&lwork));
879: PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
880: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_X],&X));
881: #if !defined(PETSC_USE_COMPLEX)
882: PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,wi+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,NULL,&info));
883: #else
884: PetscCallBLAS("LAPACKgeevx",LAPACKgeevx_("B","N","V","N",&n_,H,&n_,wr+ds->l,NULL,&ld,X+off,&ld,&ilo,&ihi,scale,&nrm,rcde,rcdv,ds->work+nwu,&lwork,ds->rwork+nwru,&info));
886: /* Sort to have consecutive conjugate pairs
887: Separate real and imaginary part of complex eigenvectors*/
888: for (i=ds->l;i<ds->n;i++) {
889: j=i+1;
890: while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
891: if (j==ds->n) {
892: PetscCheck(PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON,PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
893: wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
894: for (k=ds->l;k<ds->n;k++) {
895: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
896: }
897: } else { /* complex eigenvalue */
898: if (j!=i+1) {
899: wr[j] = wr[i+1];
900: PetscCall(PetscArraycpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld));
901: }
902: if (PetscImaginaryPart(wr[i])<0) {
903: wr[i] = PetscConj(wr[i]);
904: for (k=ds->l;k<ds->n;k++) {
905: X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
906: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
907: }
908: } else {
909: for (k=ds->l;k<ds->n;k++) {
910: X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
911: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
912: }
913: }
914: wr[i+1] = PetscConj(wr[i]);
915: i++;
916: }
917: }
918: #endif
919: SlepcCheckLapackInfo("geevx",info);
920: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_X],&X));
921: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
922: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
923: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
924: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
926: /* Compute real s-orthonormal basis */
927: PetscCall(DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_FALSE));
929: /* Recover eigenvalues from diagonal */
930: PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
931: #if defined(PETSC_USE_COMPLEX)
932: if (wi) {
933: for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
934: }
935: #endif
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: static PetscErrorCode DSGetTruncateSize_GHIEP(DS ds,PetscInt l,PetscInt n,PetscInt *k)
940: {
941: PetscReal *T;
943: PetscFunctionBegin;
944: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
945: if (T[l+(*k)-1+ds->ld] !=0.0) {
946: if (l+(*k)<n-1) (*k)++;
947: else (*k)--;
948: }
949: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: static PetscErrorCode DSTruncate_GHIEP(DS ds,PetscInt n,PetscBool trim)
954: {
955: PetscInt i,ld=ds->ld,l=ds->l;
956: PetscScalar *A;
957: PetscReal *T,*b,*r,*omega;
959: PetscFunctionBegin;
960: if (ds->compact) {
961: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
962: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&omega));
963: #if defined(PETSC_USE_DEBUG)
964: /* make sure diagonal 2x2 block is not broken */
965: PetscCheck(ds->state<DS_STATE_CONDENSED || n==0 || n==ds->n || T[n-1+ld]==0.0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"The given size would break a 2x2 block, call DSGetTruncateSize() first");
966: #endif
967: }
968: if (trim) {
969: if (!ds->compact && ds->extrarow) { /* clean extra row */
970: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
971: for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
972: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
973: }
974: ds->l = 0;
975: ds->k = 0;
976: ds->n = n;
977: ds->t = ds->n; /* truncated length equal to the new dimension */
978: } else {
979: if (!ds->compact && ds->extrarow && ds->k==ds->n) {
980: /* copy entries of extra row to the new position, then clean last row */
981: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
982: for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
983: for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
984: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
985: }
986: if (ds->compact) {
987: b = T+ld;
988: r = T+2*ld;
989: b[n-1] = r[n-1];
990: b[n] = b[ds->n];
991: omega[n] = omega[ds->n];
992: }
993: ds->k = (ds->extrarow)? n: 0;
994: ds->t = ds->n; /* truncated length equal to previous dimension */
995: ds->n = n;
996: }
997: if (ds->compact) {
998: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
999: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&omega));
1000: }
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: #if !defined(PETSC_HAVE_MPIUNI)
1005: static PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
1006: {
1007: PetscScalar *A,*B,*Q;
1008: PetscReal *T,*D;
1009: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
1010: PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
1012: PetscFunctionBegin;
1013: if (ds->compact) kr = 4*ld;
1014: else k = 2*(ds->n-l)*ld;
1015: if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
1016: if (eigr) k += (ds->n-l);
1017: if (eigi) k += (ds->n-l);
1018: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
1019: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
1020: PetscCall(PetscMPIIntCast(ds->n-l,&n));
1021: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
1022: PetscCall(PetscMPIIntCast(ld*3,&ld3));
1023: PetscCall(PetscMPIIntCast(ld,&ld_));
1024: if (ds->compact) {
1025: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
1026: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
1027: } else {
1028: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
1029: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
1030: }
1031: if (ds->state>DS_STATE_RAW) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
1032: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
1033: if (!rank) {
1034: if (ds->compact) {
1035: PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1036: PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1037: } else {
1038: PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1039: PetscCallMPI(MPI_Pack(B+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1040: }
1041: if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1042: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1043: #if !defined(PETSC_USE_COMPLEX)
1044: if (eigi) PetscCallMPI(MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
1045: #endif
1046: }
1047: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
1048: if (rank) {
1049: if (ds->compact) {
1050: PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1051: PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
1052: } else {
1053: PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1054: PetscCallMPI(MPI_Unpack(ds->work,size,&off,B+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1055: }
1056: if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1057: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1058: #if !defined(PETSC_USE_COMPLEX)
1059: if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
1060: #endif
1061: }
1062: if (ds->compact) {
1063: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
1064: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
1065: } else {
1066: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
1067: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
1068: }
1069: if (ds->state>DS_STATE_RAW) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1072: #endif
1074: static PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
1075: {
1076: PetscFunctionBegin;
1077: if ((m==DS_MAT_A && !ds->extrarow) || m==DS_MAT_B) *flg = PETSC_TRUE;
1078: else *flg = PETSC_FALSE;
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }
1082: /*MC
1083: DSGHIEP - Dense Generalized Hermitian Indefinite Eigenvalue Problem.
1085: Level: beginner
1087: Notes:
1088: The problem is expressed as A*X = B*X*Lambda, where both A and B are
1089: real symmetric (or complex Hermitian) and possibly indefinite. Lambda
1090: is a diagonal matrix whose diagonal elements are the arguments of DSSolve().
1091: After solve, A is overwritten with Lambda. Note that in the case of real
1092: scalars, A is overwritten with a real representation of Lambda, i.e.,
1093: complex conjugate eigenvalue pairs are stored as a 2x2 block in the
1094: quasi-diagonal matrix.
1096: In the intermediate state A is reduced to tridiagonal form and B is
1097: transformed into a signature matrix. In compact storage format, these
1098: matrices are stored in T and D, respectively.
1100: Used DS matrices:
1101: + DS_MAT_A - first problem matrix
1102: . DS_MAT_B - second problem matrix
1103: . DS_MAT_T - symmetric tridiagonal matrix of the reduced pencil
1104: . DS_MAT_D - diagonal matrix (signature) of the reduced pencil
1105: - DS_MAT_Q - pseudo-orthogonal transformation that reduces (A,B) to
1106: tridiagonal-diagonal form (intermediate step) or a real basis of eigenvectors
1108: Implemented methods:
1109: + 0 - QR iteration plus inverse iteration for the eigenvectors
1110: . 1 - HZ iteration
1111: - 2 - QR iteration plus pseudo-orthogonalization for the eigenvectors
1113: References:
1114: . 1. - C. Campos and J. E. Roman, "Restarted Q-Arnoldi-type methods exploiting
1115: symmetry in quadratic eigenvalue problems", BIT Numer. Math. 56(4):1213-1236, 2016.
1117: .seealso: DSCreate(), DSSetType(), DSType
1118: M*/
1119: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
1120: {
1121: PetscFunctionBegin;
1122: ds->ops->allocate = DSAllocate_GHIEP;
1123: ds->ops->view = DSView_GHIEP;
1124: ds->ops->vectors = DSVectors_GHIEP;
1125: ds->ops->solve[0] = DSSolve_GHIEP_QR_II;
1126: ds->ops->solve[1] = DSSolve_GHIEP_HZ;
1127: ds->ops->solve[2] = DSSolve_GHIEP_QR;
1128: ds->ops->sort = DSSort_GHIEP;
1129: #if !defined(PETSC_HAVE_MPIUNI)
1130: ds->ops->synchronize = DSSynchronize_GHIEP;
1131: #endif
1132: ds->ops->gettruncatesize = DSGetTruncateSize_GHIEP;
1133: ds->ops->truncate = DSTruncate_GHIEP;
1134: ds->ops->update = DSUpdateExtraRow_GHIEP;
1135: ds->ops->hermitian = DSHermitian_GHIEP;
1136: PetscFunctionReturn(PETSC_SUCCESS);
1137: }