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 : HZ iteration for generalized symmetric-indefinite eigenproblem.
12 : Based on Matlab code from David Watkins.
13 :
14 : References:
15 :
16 : [1] D.S. Watkins, The Matrix Eigenvalue Problem: GR and Krylov Subspace
17 : Methods, SIAM, 2007.
18 :
19 : [2] M.A. Brebner, J. Grad, "Eigenvalues of Ax = lambda Bx for real
20 : symmetric matrices A and B computed by reduction to pseudosymmetric
21 : form and the HR process", Linear Alg. Appl. 43:99-118, 1982.
22 : */
23 :
24 : #include <slepc/private/dsimpl.h>
25 : #include <slepcblaslapack.h>
26 :
27 : /*
28 : Sets up a 2-by-2 matrix to eliminate y in the vector [x y]'.
29 : Transformation is rotator if sygn = 1 and hyperbolic if sygn = -1.
30 : */
31 355 : static PetscErrorCode UnifiedRotation(PetscReal x,PetscReal y,PetscReal sygn,PetscReal *rot,PetscReal *rcond,PetscBool *swap)
32 : {
33 355 : PetscReal nrm,c,s;
34 :
35 355 : PetscFunctionBegin;
36 355 : *swap = PETSC_FALSE;
37 355 : if (y == 0) {
38 0 : rot[0] = 1.0; rot[1] = 0.0; rot[2] = 0.0; rot[3] = 1.0;
39 0 : *rcond = 1.0;
40 : } else {
41 355 : nrm = PetscMax(PetscAbs(x),PetscAbs(y));
42 355 : c = x/nrm; s = y/nrm;
43 355 : PetscCheck(sygn==1.0 || sygn==-1,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Value of sygn sent to transetup must be 1 or -1");
44 355 : if (sygn == 1.0) { /* set up a rotator */
45 239 : nrm = SlepcAbs(c,s);
46 239 : c = c/nrm; s = s/nrm;
47 : /* rot = [c s; -s c]; */
48 239 : rot[0] = c; rot[1] = -s; rot[2] = s; rot[3] = c;
49 239 : *rcond = 1.0;
50 : } else { /* sygn == -1, set up a hyperbolic transformation */
51 116 : nrm = c*c-s*s;
52 116 : if (nrm > 0) nrm = PetscSqrtReal(nrm);
53 : else {
54 59 : PetscCheck(nrm<0,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Breakdown in construction of hyperbolic transformation");
55 59 : nrm = PetscSqrtReal(-nrm);
56 59 : *swap = PETSC_TRUE;
57 : }
58 116 : c = c/nrm; s = s/nrm;
59 : /* rot = [c -s; -s c]; */
60 116 : rot[0] = c; rot[1] = -s; rot[2] = -s; rot[3] = c;
61 121 : *rcond = PetscAbs(PetscAbs(s)-PetscAbs(c))/(PetscAbs(s)+PetscAbs(c));
62 : }
63 : }
64 355 : PetscFunctionReturn(PETSC_SUCCESS);
65 : }
66 :
67 37 : static PetscErrorCode HZStep(PetscBLASInt ntop,PetscBLASInt nn,PetscReal tr,PetscReal dt,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscInt n,PetscInt ld,PetscBool *flag)
68 : {
69 37 : PetscBLASInt one=1;
70 37 : PetscInt k,jj,ii;
71 37 : PetscBLASInt n_;
72 37 : PetscReal bulge10,bulge20,bulge30,bulge31,bulge41,bulge42;
73 37 : PetscReal sygn,rcond=1.0,worstcond,rot[4],buf[2],t;
74 37 : PetscScalar rtmp;
75 37 : PetscBool swap;
76 :
77 37 : PetscFunctionBegin;
78 37 : *flag = PETSC_FALSE;
79 37 : worstcond = 1.0;
80 37 : PetscCall(PetscBLASIntCast(n,&n_));
81 :
82 : /* Build initial bulge that sets step in motion */
83 37 : bulge10 = dd[ntop+1]*(aa[ntop]*(aa[ntop] - dd[ntop]*tr) + dt*dd[ntop]*dd[ntop]) + dd[ntop]*bb[ntop]*bb[ntop];
84 37 : bulge20 = bb[ntop]*(dd[ntop+1]*aa[ntop] + dd[ntop]*aa[ntop+1] - dd[ntop]*dd[ntop+1]*tr);
85 37 : bulge30 = bb[ntop]*bb[ntop+1]*dd[ntop];
86 37 : bulge31 = 0.0;
87 37 : bulge41 = 0.0;
88 37 : bulge42 = 0.0;
89 :
90 : /* Chase the bulge */
91 237 : for (jj=ntop;jj<nn-1;jj++) {
92 :
93 : /* Check for trivial bulge */
94 487 : if (jj>ntop && PetscMax(PetscMax(PetscAbs(bulge10),PetscAbs(bulge20)),PetscAbs(bulge30))<PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj]) + PetscAbs(aa[jj+1]))) {
95 8 : bb[jj-1] = 0.0; /* deflate and move on */
96 :
97 : } else { /* carry out the step */
98 :
99 : /* Annihilate tip entry bulge30 */
100 192 : if (bulge30 != 0.0) {
101 :
102 : /* Make an interchange if necessary to ensure that the
103 : first transformation is othogonal, not hyperbolic. */
104 163 : if (dd[jj+1] != dd[jj+2]) { /* make an interchange */
105 45 : if (dd[jj] != dd[jj+1]) { /* interchange 1st and 2nd */
106 20 : buf[0] = bulge20; bulge20 = bulge10; bulge10 = buf[0];
107 20 : buf[0] = aa[jj]; aa[jj] = aa[jj+1]; aa[jj+1] = buf[0];
108 20 : buf[0] = bb[jj+1]; bb[jj+1] = bulge31; bulge31 = buf[0];
109 20 : buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
110 220 : for (k=0;k<n;k++) {
111 200 : rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+1)*ld]; uu[k+(jj+1)*ld] = rtmp;
112 : }
113 : } else { /* interchange 1st and 3rd */
114 25 : buf[0] = bulge30; bulge30 = bulge10; bulge10 = buf[0];
115 25 : buf[0] = aa[jj]; aa[jj] = aa[jj+2]; aa[jj+2] = buf[0];
116 25 : buf[0] = bb[jj]; bb[jj] = bb[jj+1]; bb[jj+1] = buf[0];
117 25 : buf[0] = dd[jj]; dd[jj] = dd[jj+2]; dd[jj+2] = buf[0];
118 25 : if (jj + 2 < nn-1) {
119 21 : bulge41 = bb[jj+2];
120 21 : bb[jj+2] = 0;
121 : }
122 275 : for (k=0;k<n;k++) {
123 250 : rtmp = uu[k+jj*ld]; uu[k+jj*ld] = uu[k+(jj+2)*ld]; uu[k+(jj+2)*ld] = rtmp;
124 : }
125 : }
126 : }
127 :
128 : /* Set up transforming matrix rot. */
129 163 : PetscCall(UnifiedRotation(bulge20,bulge30,1,rot,&rcond,&swap));
130 :
131 : /* Apply transforming matrix rot to T. */
132 163 : bulge20 = rot[0]*bulge20 + rot[2]*bulge30;
133 163 : buf[0] = rot[0]*bb[jj] + rot[2]*bulge31;
134 163 : buf[1] = rot[1]*bb[jj] + rot[3]*bulge31;
135 163 : bb[jj] = buf[0];
136 163 : bulge31 = buf[1];
137 163 : buf[0] = rot[0]*rot[0]*aa[jj+1] + 2.0*rot[0]*rot[2]*bb[jj+1] + rot[2]*rot[2]*aa[jj+2];
138 163 : buf[1] = rot[1]*rot[1]*aa[jj+1] + 2.0*rot[3]*rot[1]*bb[jj+1] + rot[3]*rot[3]*aa[jj+2];
139 163 : bb[jj+1] = rot[1]*rot[0]*aa[jj+1] + rot[3]*rot[2]*aa[jj+2] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj+1];
140 163 : aa[jj+1] = buf[0];
141 163 : aa[jj+2] = buf[1];
142 163 : if (jj + 2 < nn-1) {
143 126 : bulge42 = bb[jj+2]*rot[2];
144 126 : bb[jj+2] = bb[jj+2]*rot[3];
145 : }
146 :
147 : /* Accumulate transforming matrix */
148 163 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&n_,uu+(jj+1)*ld,&one,uu+(jj+2)*ld,&one,&rot[0],&rot[2]));
149 : }
150 :
151 : /* Annihilate inner entry bulge20 */
152 192 : if (bulge20 != 0.0) {
153 :
154 : /* Begin by setting up transforming matrix rot */
155 192 : sygn = dd[jj]*dd[jj+1];
156 192 : PetscCall(UnifiedRotation(bulge10,bulge20,sygn,rot,&rcond,&swap));
157 192 : if (rcond<PETSC_MACHINE_EPSILON) {
158 0 : *flag = PETSC_TRUE;
159 0 : PetscFunctionReturn(PETSC_SUCCESS);
160 : }
161 192 : if (rcond < worstcond) worstcond = rcond;
162 :
163 : /* Apply transforming matrix rot to T */
164 192 : if (jj > ntop) bb[jj-1] = rot[0]*bulge10 + rot[2]*bulge20;
165 192 : buf[0] = rot[0]*rot[0]*aa[jj] + 2*rot[0]*rot[2]*bb[jj] + rot[2]*rot[2]*aa[jj+1];
166 192 : buf[1] = rot[1]*rot[1]*aa[jj] + 2*rot[3]*rot[1]*bb[jj] + rot[3]*rot[3]*aa[jj+1];
167 192 : bb[jj] = rot[1]*rot[0]*aa[jj] + rot[3]*rot[2]*aa[jj+1] + (rot[3]*rot[0] + rot[1]*rot[2])*bb[jj];
168 192 : aa[jj] = buf[0];
169 192 : aa[jj+1] = buf[1];
170 192 : if (jj + 1 < nn-1) {
171 : /* buf = [ bulge31 bb(jj+1) ] * rot' */
172 163 : buf[0] = rot[0]*bulge31 + rot[2]*bb[jj+1];
173 163 : buf[1] = rot[1]*bulge31 + rot[3]*bb[jj+1];
174 163 : bulge31 = buf[0];
175 163 : bb[jj+1] = buf[1];
176 : }
177 192 : if (jj + 2 < nn-1) {
178 : /* buf = [bulge41 bulge42] * rot' */
179 126 : buf[0] = rot[0]*bulge41 + rot[2]*bulge42;
180 126 : buf[1] = rot[1]*bulge41 + rot[3]*bulge42;
181 126 : bulge41 = buf[0];
182 126 : bulge42 = buf[1];
183 : }
184 :
185 : /* Apply transforming matrix rot to D */
186 192 : if (swap == 1) {
187 59 : buf[0] = dd[jj]; dd[jj] = dd[jj+1]; dd[jj+1] = buf[0];
188 : }
189 :
190 : /* Accumulate transforming matrix, uu(jj:jj+1,:) = rot*uu(jj:jj+1,:) */
191 192 : if (sygn==1) {
192 76 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&n_,uu+jj*ld,&one,uu+(jj+1)*ld,&one,&rot[0],&rot[2]));
193 : } else {
194 116 : if (PetscAbsReal(rot[0])>PetscAbsReal(rot[1])) { /* Type I */
195 57 : t = rot[1]/rot[0];
196 627 : for (ii=0;ii<n;ii++) {
197 570 : uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
198 570 : uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + uu[(jj+1)*ld+ii]/rot[0];
199 : }
200 : } else { /* Type II */
201 59 : t = rot[0]/rot[1];
202 649 : for (ii=0;ii<n;ii++) {
203 590 : rtmp = uu[jj*ld+ii];
204 590 : uu[jj*ld+ii] = rot[0]*uu[jj*ld+ii] + rot[1]*uu[(jj+1)*ld+ii];
205 590 : uu[(jj+1)*ld+ii] = t*uu[jj*ld+ii] + rtmp/rot[1];
206 : }
207 : }
208 : }
209 : }
210 : }
211 :
212 : /* Adjust bulge for next step */
213 200 : bulge10 = bb[jj];
214 200 : bulge20 = bulge31;
215 200 : bulge30 = bulge41;
216 200 : bulge31 = bulge42;
217 200 : bulge41 = 0.0;
218 200 : bulge42 = 0.0;
219 : }
220 37 : PetscFunctionReturn(PETSC_SUCCESS);
221 : }
222 :
223 3 : static PetscErrorCode HZIteration(PetscBLASInt nn,PetscBLASInt cgd,PetscReal *aa,PetscReal *bb,PetscReal *dd,PetscScalar *uu,PetscBLASInt ld)
224 : {
225 3 : PetscBLASInt j2,one=1;
226 3 : PetscInt its,nits,nstop,jj,ntop,nbot,ntry;
227 3 : PetscReal htr,det,dis,dif,tn,kt,c,s,tr,dt;
228 3 : PetscBool flag=PETSC_FALSE;
229 :
230 3 : PetscFunctionBegin;
231 3 : its = 0;
232 3 : nbot = nn-1;
233 3 : nits = 0;
234 3 : nstop = 40*(nn - cgd);
235 :
236 57 : while (nbot >= cgd && nits < nstop) {
237 :
238 : /* Check for zeros on the subdiagonal */
239 54 : jj = nbot - 1;
240 265 : while (jj>=cgd && PetscAbs(bb[jj])>PETSC_MACHINE_EPSILON*(PetscAbs(aa[jj])+PetscAbs(aa[jj+1]))) jj = jj-1;
241 54 : if (jj>=cgd) bb[jj]=0;
242 54 : ntop = jj + 1; /* starting point for step */
243 54 : if (ntop == nbot) { /* isolate single eigenvalue */
244 : nbot = ntop - 1;
245 : its = 0;
246 48 : } else if (ntop+1 == nbot) { /* isolate pair of eigenvalues */
247 11 : htr = 0.5*(aa[ntop]*dd[ntop] + aa[nbot]*dd[nbot]);
248 11 : det = dd[ntop]*dd[nbot]*(aa[ntop]*aa[nbot]-bb[ntop]*bb[ntop]);
249 11 : dis = htr*htr - det;
250 11 : if (dis > 0) { /* distinct real eigenvalues */
251 8 : if (dd[ntop] == dd[nbot]) { /* separate the eigenvalues by a Jacobi rotator */
252 6 : dif = aa[ntop]-aa[nbot];
253 6 : if (2.0*PetscAbs(bb[ntop])<=dif) {
254 3 : tn = 2*bb[ntop]/dif;
255 3 : tn = tn/(1.0 + PetscSqrtReal(1.0+tn*tn));
256 : } else {
257 3 : kt = dif/(2.0*bb[ntop]);
258 3 : tn = PetscSign(kt)/(PetscAbsReal(kt)+PetscSqrtReal(1.0+kt*kt));
259 : }
260 6 : c = 1.0/PetscSqrtReal(1.0 + tn*tn);
261 6 : s = c*tn;
262 6 : aa[ntop] = aa[ntop] + tn*bb[ntop];
263 6 : aa[nbot] = aa[nbot] - tn*bb[ntop];
264 6 : bb[ntop] = 0;
265 6 : j2 = nn-cgd;
266 6 : PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,uu+ntop*ld+cgd,&one,uu+nbot*ld+cgd,&one,&c,&s));
267 : }
268 : }
269 : nbot = ntop - 1;
270 : } else { /* Do an HZ iteration */
271 37 : its = its + 1;
272 37 : nits = nits + 1;
273 37 : tr = aa[nbot-1]*dd[nbot-1] + aa[nbot]*dd[nbot];
274 37 : dt = dd[nbot-1]*dd[nbot]*(aa[nbot-1]*aa[nbot]-bb[nbot-1]*bb[nbot-1]);
275 37 : for (ntry=1;ntry<=6;ntry++) {
276 37 : PetscCall(HZStep(ntop,nbot+1,tr,dt,aa,bb,dd,uu,nn,ld,&flag));
277 37 : if (!flag) break;
278 0 : PetscCheck(ntry<6,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Unable to complete hz step after six tries");
279 0 : tr = 0.9*tr; dt = 0.81*dt;
280 : }
281 : }
282 : }
283 3 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 4 : PetscErrorCode DSSolve_GHIEP_HZ(DS ds,PetscScalar *wr,PetscScalar *wi)
287 : {
288 4 : PetscInt i,off;
289 4 : PetscBLASInt n1,ld = 0;
290 4 : const PetscScalar *A,*B;
291 4 : PetscScalar *Q;
292 4 : PetscReal *d,*e,*s;
293 :
294 4 : PetscFunctionBegin;
295 : #if !defined(PETSC_USE_COMPLEX)
296 : PetscAssertPointer(wi,3);
297 : #endif
298 4 : PetscCall(PetscBLASIntCast(ds->ld,&ld));
299 4 : n1 = ds->n - ds->l;
300 4 : off = ds->l + ds->l*ld;
301 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
302 4 : PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
303 4 : e = d + ld;
304 : #if defined(PETSC_USE_DEBUG)
305 : /* Check signature */
306 4 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
307 39 : for (i=0;i<ds->n;i++) {
308 35 : PetscReal de = (ds->compact)?s[i]:PetscRealPart(B[i*ld+i]);
309 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");
310 : }
311 4 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
312 : #endif
313 : /* Quick return */
314 4 : if (n1 == 1) {
315 1 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
316 6 : for (i=0;i<=ds->l;i++) Q[i+i*ld] = 1.0;
317 1 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
318 1 : PetscCall(DSGHIEPComplexEigs(ds,0,ds->l,wr,wi));
319 1 : if (ds->compact) {
320 1 : wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
321 : } else {
322 0 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
323 0 : PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&B));
324 0 : d[ds->l] = PetscRealPart(A[off]);
325 0 : s[ds->l] = PetscRealPart(B[off]);
326 0 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
327 0 : PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&B));
328 0 : wr[ds->l] = d[ds->l]/s[ds->l]; wi[ds->l] = 0.0;
329 : }
330 1 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
331 1 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
332 1 : PetscFunctionReturn(PETSC_SUCCESS);
333 : }
334 : /* Reduce to pseudotriadiagonal form */
335 3 : PetscCall(DSIntermediate_GHIEP(ds));
336 3 : PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
337 3 : PetscCall(HZIteration(ds->n,ds->l,d,e,s,Q,ld));
338 3 : PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
339 3 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
340 3 : PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
341 3 : if (!ds->compact) PetscCall(DSSwitchFormat_GHIEP(ds,PETSC_FALSE));
342 : /* Undo from diagonal the blocks with real eigenvalues*/
343 3 : PetscCall(DSGHIEPRealBlocks(ds));
344 :
345 : /* Recover eigenvalues from diagonal */
346 3 : PetscCall(DSGHIEPComplexEigs(ds,0,ds->n,wr,wi));
347 : #if defined(PETSC_USE_COMPLEX)
348 3 : if (wi) {
349 31 : for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
350 : }
351 : #endif
352 3 : ds->t = ds->n;
353 3 : PetscFunctionReturn(PETSC_SUCCESS);
354 : }
|
