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