Actual source code: dsghiep.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: PetscErrorCode DSAllocate_GHIEP(DS ds,PetscInt ld)
15: {
19: DSAllocateMat_Private(ds,DS_MAT_A);
20: DSAllocateMat_Private(ds,DS_MAT_B);
21: DSAllocateMat_Private(ds,DS_MAT_Q);
22: DSAllocateMatReal_Private(ds,DS_MAT_T);
23: DSAllocateMatReal_Private(ds,DS_MAT_D);
24: PetscFree(ds->perm);
25: PetscMalloc1(ld,&ds->perm);
26: PetscLogObjectMemory((PetscObject)ds,ld*sizeof(PetscInt));
27: return(0);
28: }
30: PetscErrorCode DSSwitchFormat_GHIEP(DS ds,PetscBool tocompact)
31: {
33: PetscReal *T,*S;
34: PetscScalar *A,*B;
35: PetscInt i,n,ld;
38: A = ds->mat[DS_MAT_A];
39: B = ds->mat[DS_MAT_B];
40: T = ds->rmat[DS_MAT_T];
41: S = ds->rmat[DS_MAT_D];
42: n = ds->n;
43: ld = ds->ld;
44: if (tocompact) { /* switch from dense (arrow) to compact storage */
45: PetscMemzero(T,3*ld*sizeof(PetscReal));
46: PetscMemzero(S,ld*sizeof(PetscReal));
47: for (i=0;i<n-1;i++) {
48: T[i] = PetscRealPart(A[i+i*ld]);
49: T[ld+i] = PetscRealPart(A[i+1+i*ld]);
50: S[i] = PetscRealPart(B[i+i*ld]);
51: }
52: T[n-1] = PetscRealPart(A[n-1+(n-1)*ld]);
53: S[n-1] = PetscRealPart(B[n-1+(n-1)*ld]);
54: for (i=ds->l;i<ds->k;i++) T[2*ld+i] = PetscRealPart(A[ds->k+i*ld]);
55: } else { /* switch from compact (arrow) to dense storage */
56: PetscMemzero(A,ld*ld*sizeof(PetscScalar));
57: PetscMemzero(B,ld*ld*sizeof(PetscScalar));
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: return(0);
72: }
74: PetscErrorCode DSView_GHIEP(DS ds,PetscViewer viewer)
75: {
76: PetscErrorCode ierr;
77: PetscViewerFormat format;
78: PetscInt i,j;
79: PetscReal value;
80: const char *methodname[] = {
81: "QR + Inverse Iteration",
82: "HZ method",
83: "QR"
84: };
85: const int nmeth=sizeof(methodname)/sizeof(methodname[0]);
88: PetscViewerGetFormat(viewer,&format);
89: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
90: if (ds->method<nmeth) {
91: PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]);
92: }
93: return(0);
94: }
95: if (ds->compact) {
96: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
97: if (format == PETSC_VIEWER_ASCII_MATLAB) {
98: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",ds->n,ds->n);
99: PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",3*ds->n);
100: PetscViewerASCIIPrintf(viewer,"zzz = [\n");
101: for (i=0;i<ds->n;i++) {
102: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,i+1,(double)*(ds->rmat[DS_MAT_T]+i));
103: }
104: for (i=0;i<ds->n-1;i++) {
105: if (*(ds->rmat[DS_MAT_T]+ds->ld+i) !=0 && i!=ds->k-1) {
106: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+2,i+1,(double)*(ds->rmat[DS_MAT_T]+ds->ld+i));
107: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,i+2,(double)*(ds->rmat[DS_MAT_T]+ds->ld+i));
108: }
109: }
110: for (i = ds->l;i<ds->k;i++) {
111: if (*(ds->rmat[DS_MAT_T]+2*ds->ld+i)) {
112: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",ds->k+1,i+1,(double)*(ds->rmat[DS_MAT_T]+2*ds->ld+i));
113: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,ds->k+1,(double)*(ds->rmat[DS_MAT_T]+2*ds->ld+i));
114: }
115: }
116: PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_A]);
118: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",ds->n,ds->n);
119: PetscViewerASCIIPrintf(viewer,"omega = zeros(%D,3);\n",3*ds->n);
120: PetscViewerASCIIPrintf(viewer,"omega = [\n");
121: for (i=0;i<ds->n;i++) {
122: PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,i+1,(double)*(ds->rmat[DS_MAT_D]+i));
123: }
124: PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_B]);
126: } else {
127: PetscViewerASCIIPrintf(viewer,"T\n");
128: for (i=0;i<ds->n;i++) {
129: for (j=0;j<ds->n;j++) {
130: if (i==j) value = *(ds->rmat[DS_MAT_T]+i);
131: else if (i==j+1 || j==i+1) value = *(ds->rmat[DS_MAT_T]+ds->ld+PetscMin(i,j));
132: else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = *(ds->rmat[DS_MAT_T]+2*ds->ld+PetscMin(i,j));
133: else value = 0.0;
134: PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value);
135: }
136: PetscViewerASCIIPrintf(viewer,"\n");
137: }
138: PetscViewerASCIIPrintf(viewer,"omega\n");
139: for (i=0;i<ds->n;i++) {
140: for (j=0;j<ds->n;j++) {
141: if (i==j) value = *(ds->rmat[DS_MAT_D]+i);
142: else value = 0.0;
143: PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value);
144: }
145: PetscViewerASCIIPrintf(viewer,"\n");
146: }
147: }
148: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
149: PetscViewerFlush(viewer);
150: } else {
151: DSViewMat(ds,viewer,DS_MAT_A);
152: DSViewMat(ds,viewer,DS_MAT_B);
153: }
154: if (ds->state>DS_STATE_INTERMEDIATE) {
155: DSViewMat(ds,viewer,DS_MAT_Q);
156: }
157: return(0);
158: }
160: static PetscErrorCode DSVectors_GHIEP_Eigen_Some(DS ds,PetscInt *idx,PetscReal *rnorm)
161: {
162: #if defined(SLEPC_MISSING_LAPACK_LAG2)
164: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LAG2 - Lapack routine is unavailable");
165: #else
167: PetscReal b[4],M[4],d1,d2,s1,s2,e;
168: PetscReal scal1,scal2,wr1,wr2,wi,ep,norm;
169: PetscScalar *Q,*X,Y[4],alpha,zeroS = 0.0;
170: PetscInt k;
171: PetscBLASInt two = 2,n_,ld,one=1;
172: #if !defined(PETSC_USE_COMPLEX)
173: PetscBLASInt four=4;
174: #endif
177: X = ds->mat[DS_MAT_X];
178: Q = ds->mat[DS_MAT_Q];
179: k = *idx;
180: PetscBLASIntCast(ds->n,&n_);
181: PetscBLASIntCast(ds->ld,&ld);
182: if (k < ds->n-1) e = (ds->compact)?*(ds->rmat[DS_MAT_T]+ld+k):PetscRealPart(*(ds->mat[DS_MAT_A]+(k+1)+ld*k));
183: else e = 0.0;
184: if (e == 0.0) { /* Real */
185: if (ds->state>=DS_STATE_CONDENSED) {
186: PetscMemcpy(X+k*ld,Q+k*ld,ld*sizeof(PetscScalar));
187: } else {
188: PetscMemzero(X+k*ds->ld,ds->ld*sizeof(PetscScalar));
189: X[k+k*ds->ld] = 1.0;
190: }
191: if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
192: } else { /* 2x2 block */
193: if (ds->compact) {
194: s1 = *(ds->rmat[DS_MAT_D]+k);
195: d1 = *(ds->rmat[DS_MAT_T]+k);
196: s2 = *(ds->rmat[DS_MAT_D]+k+1);
197: d2 = *(ds->rmat[DS_MAT_T]+k+1);
198: } else {
199: s1 = PetscRealPart(*(ds->mat[DS_MAT_B]+k*ld+k));
200: d1 = PetscRealPart(*(ds->mat[DS_MAT_A]+k+k*ld));
201: s2 = PetscRealPart(*(ds->mat[DS_MAT_B]+(k+1)*ld+k+1));
202: d2 = PetscRealPart(*(ds->mat[DS_MAT_A]+k+1+(k+1)*ld));
203: }
204: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
205: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
206: ep = LAPACKlamch_("S");
207: /* Compute eigenvalues of the block */
208: PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
209: if (wi==0.0) SETERRQ(PETSC_COMM_SELF,1,"Real block in DSVectors_GHIEP");
210: else { /* Complex eigenvalues */
211: if (scal1<ep) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
212: wr1 /= scal1;
213: wi /= scal1;
214: #if !defined(PETSC_USE_COMPLEX)
215: if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
216: Y[0] = wr1-s2*d2; Y[1] = s2*e; Y[2] = wi; Y[3] = 0.0;
217: } else {
218: Y[0] = s1*e; Y[1] = wr1-s1*d1; Y[2] = 0.0; Y[3] = wi;
219: }
220: norm = BLASnrm2_(&four,Y,&one);
221: norm = 1.0/norm;
222: if (ds->state >= DS_STATE_CONDENSED) {
223: alpha = norm;
224: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n_,&two,&two,&alpha,ds->mat[DS_MAT_Q]+k*ld,&ld,Y,&two,&zeroS,X+k*ld,&ld));
225: if (rnorm) *rnorm = SlepcAbsEigenvalue(X[ds->n-1+k*ld],X[ds->n-1+(k+1)*ld]);
226: } else {
227: PetscMemzero(X+k*ld,2*ld*sizeof(PetscScalar));
228: X[k*ld+k] = Y[0]*norm;
229: X[k*ld+k+1] = Y[1]*norm;
230: X[(k+1)*ld+k] = Y[2]*norm;
231: X[(k+1)*ld+k+1] = Y[3]*norm;
232: }
233: #else
234: if (SlepcAbs(s1*d1-wr1,wi)<SlepcAbs(s2*d2-wr1,wi)) {
235: Y[0] = wr1-s2*d2+PETSC_i*wi;
236: Y[1] = s2*e;
237: } else {
238: Y[0] = s1*e;
239: Y[1] = wr1-s1*d1+PETSC_i*wi;
240: }
241: norm = BLASnrm2_(&two,Y,&one);
242: norm = 1.0/norm;
243: if (ds->state >= DS_STATE_CONDENSED) {
244: alpha = norm;
245: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n_,&two,&alpha,ds->mat[DS_MAT_Q]+k*ld,&ld,Y,&one,&zeroS,X+k*ld,&one));
246: if (rnorm) *rnorm = PetscAbsScalar(X[ds->n-1+k*ld]);
247: } else {
248: PetscMemzero(X+k*ld,2*ld*sizeof(PetscScalar));
249: X[k*ld+k] = Y[0]*norm;
250: X[k*ld+k+1] = Y[1]*norm;
251: }
252: X[(k+1)*ld+k] = PetscConj(X[k*ld+k]);
253: X[(k+1)*ld+k+1] = PetscConj(X[k*ld+k+1]);
254: #endif
255: (*idx)++;
256: }
257: }
258: return(0);
259: #endif
260: }
262: PetscErrorCode DSVectors_GHIEP(DS ds,DSMatType mat,PetscInt *k,PetscReal *rnorm)
263: {
264: PetscInt i;
265: PetscReal e;
269: switch (mat) {
270: case DS_MAT_X:
271: case DS_MAT_Y:
272: if (k) {
273: DSVectors_GHIEP_Eigen_Some(ds,k,rnorm);
274: } else {
275: for (i=0; i<ds->n; i++) {
276: e = (ds->compact)?*(ds->rmat[DS_MAT_T]+ds->ld+i):PetscRealPart(*(ds->mat[DS_MAT_A]+(i+1)+ds->ld*i));
277: if (e == 0.0) { /* real */
278: if (ds->state >= DS_STATE_CONDENSED) {
279: PetscMemcpy(ds->mat[mat]+i*ds->ld,ds->mat[DS_MAT_Q]+i*ds->ld,ds->ld*sizeof(PetscScalar));
280: } else {
281: PetscMemzero(ds->mat[mat]+i*ds->ld,ds->ld*sizeof(PetscScalar));
282: *(ds->mat[mat]+i+i*ds->ld) = 1.0;
283: }
284: } else {
285: DSVectors_GHIEP_Eigen_Some(ds,&i,rnorm);
286: }
287: }
288: }
289: break;
290: case DS_MAT_U:
291: case DS_MAT_VT:
292: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
293: break;
294: default:
295: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
296: }
297: return(0);
298: }
300: /*
301: Extract the eigenvalues contained in the block-diagonal of the indefinite problem.
302: Only the index range n0..n1 is processed.
303: */
304: PetscErrorCode DSGHIEPComplexEigs(DS ds,PetscInt n0,PetscInt n1,PetscScalar *wr,PetscScalar *wi)
305: {
306: #if defined(SLEPC_MISSING_LAPACK_LAG2)
308: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LAG2 - Lapack routine is unavailable");
309: #else
310: PetscInt k,ld;
311: PetscBLASInt two=2;
312: PetscScalar *A,*B;
313: PetscReal *D,*T;
314: PetscReal b[4],M[4],d1,d2,s1,s2,e;
315: PetscReal scal1,scal2,ep,wr1,wr2,wi1;
318: ld = ds->ld;
319: A = ds->mat[DS_MAT_A];
320: B = ds->mat[DS_MAT_B];
321: D = ds->rmat[DS_MAT_D];
322: T = ds->rmat[DS_MAT_T];
323: for (k=n0;k<n1;k++) {
324: if (k < n1-1) e = (ds->compact)?T[ld+k]:PetscRealPart(A[(k+1)+ld*k]);
325: else e = 0.0;
326: if (e==0.0) { /* real eigenvalue */
327: wr[k] = (ds->compact)?T[k]/D[k]:A[k+k*ld]/B[k+k*ld];
328: #if !defined(PETSC_USE_COMPLEX)
329: wi[k] = 0.0 ;
330: #endif
331: } else { /* diagonal block */
332: if (ds->compact) {
333: s1 = D[k];
334: d1 = T[k];
335: s2 = D[k+1];
336: d2 = T[k+1];
337: } else {
338: s1 = PetscRealPart(B[k*ld+k]);
339: d1 = PetscRealPart(A[k+k*ld]);
340: s2 = PetscRealPart(B[(k+1)*ld+k+1]);
341: d2 = PetscRealPart(A[k+1+(k+1)*ld]);
342: }
343: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
344: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
345: ep = LAPACKlamch_("S");
346: /* Compute eigenvalues of the block */
347: PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi1));
348: if (scal1<ep) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
349: wr[k] = wr1/scal1;
350: if (wi1==0.0) { /* Real eigenvalues */
351: if (scal2<ep) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
352: wr[k+1] = wr2/scal2;
353: #if !defined(PETSC_USE_COMPLEX)
354: wi[k] = 0.0;
355: wi[k+1] = 0.0;
356: #endif
357: } else { /* Complex eigenvalues */
358: #if !defined(PETSC_USE_COMPLEX)
359: wr[k+1] = wr[k];
360: wi[k] = wi1/scal1;
361: wi[k+1] = -wi[k];
362: #else
363: wr[k] += PETSC_i*wi1/scal1;
364: wr[k+1] = PetscConj(wr[k]);
365: #endif
366: }
367: k++;
368: }
369: }
370: #if defined(PETSC_USE_COMPLEX)
371: if (wi) {
372: for (k=n0;k<n1;k++) wi[k] = 0.0;
373: }
374: #endif
375: return(0);
376: #endif
377: }
379: PetscErrorCode DSSort_GHIEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
380: {
382: PetscInt n,i,*perm;
383: PetscReal *d,*e,*s;
386: #if !defined(PETSC_USE_COMPLEX)
388: #endif
389: n = ds->n;
390: d = ds->rmat[DS_MAT_T];
391: e = d + ds->ld;
392: s = ds->rmat[DS_MAT_D];
393: DSAllocateWork_Private(ds,ds->ld,ds->ld,0);
394: perm = ds->perm;
395: if (!rr) {
396: rr = wr;
397: ri = wi;
398: }
399: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_TRUE);
400: if (!ds->compact) { DSSwitchFormat_GHIEP(ds,PETSC_TRUE); }
401: PetscMemcpy(ds->work,wr,n*sizeof(PetscScalar));
402: for (i=ds->l;i<n;i++) wr[i] = *(ds->work+perm[i]);
403: #if !defined(PETSC_USE_COMPLEX)
404: PetscMemcpy(ds->work,wi,n*sizeof(PetscScalar));
405: for (i=ds->l;i<n;i++) wi[i] = *(ds->work+perm[i]);
406: #endif
407: PetscMemcpy(ds->rwork,s,n*sizeof(PetscReal));
408: for (i=ds->l;i<n;i++) s[i] = *(ds->rwork+perm[i]);
409: PetscMemcpy(ds->rwork,d,n*sizeof(PetscReal));
410: for (i=ds->l;i<n;i++) d[i] = *(ds->rwork+perm[i]);
411: PetscMemcpy(ds->rwork,e,(n-1)*sizeof(PetscReal));
412: PetscMemzero(e+ds->l,(n-1-ds->l)*sizeof(PetscScalar));
413: for (i=ds->l;i<n-1;i++) {
414: if (perm[i]<n-1) e[i] = *(ds->rwork+perm[i]);
415: }
416: if (!ds->compact) { DSSwitchFormat_GHIEP(ds,PETSC_FALSE); }
417: DSPermuteColumns_Private(ds,ds->l,n,DS_MAT_Q,perm);
418: return(0);
419: }
421: /*
422: Get eigenvectors with inverse iteration.
423: The system matrix is in Hessenberg form.
424: */
425: PetscErrorCode DSGHIEPInverseIteration(DS ds,PetscScalar *wr,PetscScalar *wi)
426: {
427: #if defined(SLEPC_MISSING_LAPACK_HSEIN)
429: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"HSEIN - Lapack routine is unavailable");
430: #else
432: PetscInt i,off;
433: PetscBLASInt *select,*infoC,ld,n1,mout,info;
434: PetscScalar *A,*B,*H,*X;
435: PetscReal *s,*d,*e;
436: #if defined(PETSC_USE_COMPLEX)
437: PetscInt j;
438: #endif
441: PetscBLASIntCast(ds->ld,&ld);
442: PetscBLASIntCast(ds->n-ds->l,&n1);
443: DSAllocateWork_Private(ds,ld*ld+2*ld,ld,2*ld);
444: DSAllocateMat_Private(ds,DS_MAT_W);
445: A = ds->mat[DS_MAT_A];
446: B = ds->mat[DS_MAT_B];
447: H = ds->mat[DS_MAT_W];
448: s = ds->rmat[DS_MAT_D];
449: d = ds->rmat[DS_MAT_T];
450: e = d + ld;
451: select = ds->iwork;
452: infoC = ds->iwork + ld;
453: off = ds->l+ds->l*ld;
454: if (ds->compact) {
455: H[off] = d[ds->l]*s[ds->l];
456: H[off+ld] = e[ds->l]*s[ds->l];
457: for (i=ds->l+1;i<ds->n-1;i++) {
458: H[i+(i-1)*ld] = e[i-1]*s[i];
459: H[i+i*ld] = d[i]*s[i];
460: H[i+(i+1)*ld] = e[i]*s[i];
461: }
462: H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
463: H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
464: } else {
465: s[ds->l] = PetscRealPart(B[off]);
466: H[off] = A[off]*s[ds->l];
467: H[off+ld] = A[off+ld]*s[ds->l];
468: for (i=ds->l+1;i<ds->n-1;i++) {
469: s[i] = PetscRealPart(B[i+i*ld]);
470: H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
471: H[i+i*ld] = A[i+i*ld]*s[i];
472: H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
473: }
474: s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
475: H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
476: H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
477: }
478: DSAllocateMat_Private(ds,DS_MAT_X);
479: X = ds->mat[DS_MAT_X];
480: for (i=0;i<n1;i++) select[i] = 1;
481: #if !defined(PETSC_USE_COMPLEX)
482: PetscStackCallBLAS("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));
483: #else
484: PetscStackCallBLAS("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));
486: /* Separate real and imaginary part of complex eigenvectors */
487: for (j=ds->l;j<ds->n;j++) {
488: if (PetscAbsReal(PetscImaginaryPart(wr[j])) > PetscAbsScalar(wr[j])*PETSC_SQRT_MACHINE_EPSILON) {
489: for (i=ds->l;i<ds->n;i++) {
490: X[i+(j+1)*ds->ld] = PetscImaginaryPart(X[i+j*ds->ld]);
491: X[i+j*ds->ld] = PetscRealPart(X[i+j*ds->ld]);
492: }
493: j++;
494: }
495: }
496: #endif
497: SlepcCheckLapackInfo("hsein",info);
498: DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_TRUE);
499: return(0);
500: #endif
501: }
503: /*
504: Undo 2x2 blocks that have real eigenvalues.
505: */
506: PetscErrorCode DSGHIEPRealBlocks(DS ds)
507: {
508: #if defined(SLEPC_MISSING_LAPACK_LAG2)
510: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"LAG2 - Lapack routine is unavailable");
511: #else
513: PetscInt i;
514: PetscReal e,d1,d2,s1,s2,ss1,ss2,t,dd,ss;
515: PetscReal maxy,ep,scal1,scal2,snorm;
516: PetscReal *T,*D,b[4],M[4],wr1,wr2,wi;
517: PetscScalar *A,*B,Y[4],oneS = 1.0,zeroS = 0.0;
518: PetscBLASInt m,two=2,ld;
519: PetscBool isreal;
522: PetscBLASIntCast(ds->ld,&ld);
523: PetscBLASIntCast(ds->n-ds->l,&m);
524: A = ds->mat[DS_MAT_A];
525: B = ds->mat[DS_MAT_B];
526: T = ds->rmat[DS_MAT_T];
527: D = ds->rmat[DS_MAT_D];
528: DSAllocateWork_Private(ds,2*m,0,0);
529: for (i=ds->l;i<ds->n-1;i++) {
530: e = (ds->compact)?T[ld+i]:PetscRealPart(A[(i+1)+ld*i]);
531: if (e != 0.0) { /* 2x2 block */
532: if (ds->compact) {
533: s1 = D[i];
534: d1 = T[i];
535: s2 = D[i+1];
536: d2 = T[i+1];
537: } else {
538: s1 = PetscRealPart(B[i*ld+i]);
539: d1 = PetscRealPart(A[i*ld+i]);
540: s2 = PetscRealPart(B[(i+1)*ld+i+1]);
541: d2 = PetscRealPart(A[(i+1)*ld+i+1]);
542: }
543: isreal = PETSC_FALSE;
544: if (s1==s2) { /* apply a Jacobi rotation to compute the eigendecomposition */
545: dd = d1-d2;
546: if (2*PetscAbsReal(e) <= dd) {
547: t = 2*e/dd;
548: t = t/(1 + PetscSqrtReal(1+t*t));
549: } else {
550: t = dd/(2*e);
551: ss = (t>=0)?1.0:-1.0;
552: t = ss/(PetscAbsReal(t)+PetscSqrtReal(1+t*t));
553: }
554: Y[0] = 1/PetscSqrtReal(1 + t*t); Y[3] = Y[0]; /* c */
555: Y[1] = Y[0]*t; Y[2] = -Y[1]; /* s */
556: wr1 = d1+t*e;
557: wr2 = d2-t*e;
558: ss1 = s1; ss2 = s2;
559: isreal = PETSC_TRUE;
560: } else {
561: ss1 = 1.0; ss2 = 1.0,
562: M[0] = d1; M[1] = e; M[2] = e; M[3]= d2;
563: b[0] = s1; b[1] = 0.0; b[2] = 0.0; b[3] = s2;
564: ep = LAPACKlamch_("S");
566: /* Compute eigenvalues of the block */
567: PetscStackCallBLAS("LAPACKlag2",LAPACKlag2_(M,&two,b,&two,&ep,&scal1,&scal2,&wr1,&wr2,&wi));
568: if (wi==0.0) { /* Real eigenvalues */
569: isreal = PETSC_TRUE;
570: if (scal1<ep||scal2<ep) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Nearly infinite eigenvalue");
571: wr1 /= scal1;
572: wr2 /= scal2;
573: if (PetscAbsReal(s1*d1-wr1)<PetscAbsReal(s2*d2-wr1)) {
574: Y[0] = wr1-s2*d2;
575: Y[1] = s2*e;
576: } else {
577: Y[0] = s1*e;
578: Y[1] = wr1-s1*d1;
579: }
580: /* normalize with a signature*/
581: maxy = PetscMax(PetscAbsScalar(Y[0]),PetscAbsScalar(Y[1]));
582: scal1 = PetscRealPart(Y[0])/maxy;
583: scal2 = PetscRealPart(Y[1])/maxy;
584: snorm = scal1*scal1*s1 + scal2*scal2*s2;
585: if (snorm<0) { ss1 = -1.0; snorm = -snorm; }
586: snorm = maxy*PetscSqrtReal(snorm);
587: Y[0] = Y[0]/snorm;
588: Y[1] = Y[1]/snorm;
589: if (PetscAbsReal(s1*d1-wr2)<PetscAbsReal(s2*d2-wr2)) {
590: Y[2] = wr2-s2*d2;
591: Y[3] = s2*e;
592: } else {
593: Y[2] = s1*e;
594: Y[3] = wr2-s1*d1;
595: }
596: maxy = PetscMax(PetscAbsScalar(Y[2]),PetscAbsScalar(Y[3]));
597: scal1 = PetscRealPart(Y[2])/maxy;
598: scal2 = PetscRealPart(Y[3])/maxy;
599: snorm = scal1*scal1*s1 + scal2*scal2*s2;
600: if (snorm<0) { ss2 = -1.0; snorm = -snorm; }
601: snorm = maxy*PetscSqrtReal(snorm); Y[2] = Y[2]/snorm; Y[3] = Y[3]/snorm;
602: }
603: wr1 *= ss1; wr2 *= ss2;
604: }
605: if (isreal) {
606: if (ds->compact) {
607: D[i] = ss1;
608: T[i] = wr1;
609: D[i+1] = ss2;
610: T[i+1] = wr2;
611: T[ld+i] = 0.0;
612: } else {
613: B[i*ld+i] = ss1;
614: A[i*ld+i] = wr1;
615: B[(i+1)*ld+i+1] = ss2;
616: A[(i+1)*ld+i+1] = wr2;
617: A[(i+1)+ld*i] = 0.0;
618: A[i+ld*(i+1)] = 0.0;
619: }
620: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&m,&two,&two,&oneS,ds->mat[DS_MAT_Q]+ds->l+i*ld,&ld,Y,&two,&zeroS,ds->work,&m));
621: PetscMemcpy(ds->mat[DS_MAT_Q]+ds->l+i*ld,ds->work,m*sizeof(PetscScalar));
622: PetscMemcpy(ds->mat[DS_MAT_Q]+ds->l+(i+1)*ld,ds->work+m,m*sizeof(PetscScalar));
623: }
624: i++;
625: }
626: }
627: return(0);
628: #endif
629: }
631: PetscErrorCode DSSolve_GHIEP_QR_II(DS ds,PetscScalar *wr,PetscScalar *wi)
632: {
633: #if defined(PETSC_MISSING_LAPACK_HSEQR)
635: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"HSEQR - Lapack routine is unavailable");
636: #else
638: PetscInt i,off;
639: PetscBLASInt n1,ld,one,info,lwork;
640: PetscScalar *H,*A,*B,*Q;
641: PetscReal *d,*e,*s;
642: #if defined(PETSC_USE_COMPLEX)
643: PetscInt j;
644: #endif
647: #if !defined(PETSC_USE_COMPLEX)
649: #endif
650: one = 1;
651: PetscBLASIntCast(ds->n-ds->l,&n1);
652: PetscBLASIntCast(ds->ld,&ld);
653: off = ds->l + ds->l*ld;
654: A = ds->mat[DS_MAT_A];
655: B = ds->mat[DS_MAT_B];
656: Q = ds->mat[DS_MAT_Q];
657: d = ds->rmat[DS_MAT_T];
658: e = ds->rmat[DS_MAT_T] + ld;
659: s = ds->rmat[DS_MAT_D];
660: DSAllocateWork_Private(ds,ld*ld,2*ld,ld*2);
661: lwork = ld*ld;
663: /* Quick return if possible */
664: if (n1 == 1) {
665: Q[off] = 1.0;
666: if (!ds->compact) {
667: d[ds->l] = PetscRealPart(A[off]);
668: s[ds->l] = PetscRealPart(B[off]);
669: }
670: wr[ds->l] = d[ds->l]/s[ds->l];
671: if (wi) wi[ds->l] = 0.0;
672: return(0);
673: }
674: /* Reduce to pseudotriadiagonal form */
675: DSIntermediate_GHIEP(ds);
677: /* Compute Eigenvalues (QR)*/
678: DSAllocateMat_Private(ds,DS_MAT_W);
679: H = ds->mat[DS_MAT_W];
680: if (ds->compact) {
681: H[off] = d[ds->l]*s[ds->l];
682: H[off+ld] = e[ds->l]*s[ds->l];
683: for (i=ds->l+1;i<ds->n-1;i++) {
684: H[i+(i-1)*ld] = e[i-1]*s[i];
685: H[i+i*ld] = d[i]*s[i];
686: H[i+(i+1)*ld] = e[i]*s[i];
687: }
688: H[ds->n-1+(ds->n-2)*ld] = e[ds->n-2]*s[ds->n-1];
689: H[ds->n-1+(ds->n-1)*ld] = d[ds->n-1]*s[ds->n-1];
690: } else {
691: s[ds->l] = PetscRealPart(B[off]);
692: H[off] = A[off]*s[ds->l];
693: H[off+ld] = A[off+ld]*s[ds->l];
694: for (i=ds->l+1;i<ds->n-1;i++) {
695: s[i] = PetscRealPart(B[i+i*ld]);
696: H[i+(i-1)*ld] = A[i+(i-1)*ld]*s[i];
697: H[i+i*ld] = A[i+i*ld]*s[i];
698: H[i+(i+1)*ld] = A[i+(i+1)*ld]*s[i];
699: }
700: s[ds->n-1] = PetscRealPart(B[ds->n-1+(ds->n-1)*ld]);
701: H[ds->n-1+(ds->n-2)*ld] = A[ds->n-1+(ds->n-2)*ld]*s[ds->n-1];
702: H[ds->n-1+(ds->n-1)*ld] = A[ds->n-1+(ds->n-1)*ld]*s[ds->n-1];
703: }
705: #if !defined(PETSC_USE_COMPLEX)
706: PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,wi+ds->l,NULL,&ld,ds->work,&lwork,&info));
707: #else
708: PetscStackCallBLAS("LAPACKhseqr",LAPACKhseqr_("E","N",&n1,&one,&n1,H+off,&ld,wr+ds->l,NULL,&ld,ds->work,&lwork,&info));
710: /* Sort to have consecutive conjugate pairs */
711: for (i=ds->l;i<ds->n;i++) {
712: j=i+1;
713: while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
714: if (j==ds->n) {
715: if (PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON) wr[i]=PetscRealPart(wr[i]);
716: else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
717: } else { /* complex eigenvalue */
718: wr[j] = wr[i+1];
719: if (PetscImaginaryPart(wr[i])<0) wr[i] = PetscConj(wr[i]);
720: wr[i+1] = PetscConj(wr[i]);
721: i++;
722: }
723: }
724: #endif
725: SlepcCheckLapackInfo("hseqr",info);
726: /* Compute Eigenvectors with Inverse Iteration */
727: DSGHIEPInverseIteration(ds,wr,wi);
729: /* Recover eigenvalues from diagonal */
730: DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
731: #if defined(PETSC_USE_COMPLEX)
732: if (wi) {
733: for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
734: }
735: #endif
736: return(0);
737: #endif
738: }
740: PetscErrorCode DSSolve_GHIEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
741: {
742: #if defined(SLEPC_MISSING_LAPACK_GEEVX)
744: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEEVX - Lapack routine is unavailable");
745: #else
747: PetscInt i,j,off,nwu=0,n,lw,lwr,nwru=0;
748: PetscBLASInt n_,ld,info,lwork,ilo,ihi;
749: PetscScalar *H,*A,*B,*Q,*X;
750: PetscReal *d,*s,*scale,nrm,*rcde,*rcdv;
751: #if defined(PETSC_USE_COMPLEX)
752: PetscInt k;
753: #endif
756: #if !defined(PETSC_USE_COMPLEX)
758: #endif
759: n = ds->n-ds->l;
760: PetscBLASIntCast(n,&n_);
761: PetscBLASIntCast(ds->ld,&ld);
762: off = ds->l + ds->l*ld;
763: A = ds->mat[DS_MAT_A];
764: B = ds->mat[DS_MAT_B];
765: Q = ds->mat[DS_MAT_Q];
766: d = ds->rmat[DS_MAT_T];
767: s = ds->rmat[DS_MAT_D];
768: lw = 14*ld+ld*ld;
769: lwr = 7*ld;
770: DSAllocateWork_Private(ds,lw,lwr,0);
771: scale = ds->rwork+nwru;
772: nwru += ld;
773: rcde = ds->rwork+nwru;
774: nwru += ld;
775: rcdv = ds->rwork+nwru;
776: nwru += ld;
777: /* Quick return if possible */
778: if (n_ == 1) {
779: Q[off] = 1.0;
780: if (!ds->compact) {
781: d[ds->l] = PetscRealPart(A[off]);
782: s[ds->l] = PetscRealPart(B[off]);
783: }
784: wr[ds->l] = d[ds->l]/s[ds->l];
785: if (wi) wi[ds->l] = 0.0;
786: return(0);
787: }
789: /* Form pseudo-symmetric matrix */
790: H = ds->work+nwu;
791: nwu += n*n;
792: PetscMemzero(H,n*n*sizeof(PetscScalar));
793: if (ds->compact) {
794: for (i=0;i<n-1;i++) {
795: H[i+i*n] = s[ds->l+i]*d[ds->l+i];
796: H[i+1+i*n] = s[ds->l+i+1]*d[ld+ds->l+i];
797: H[i+(i+1)*n] = s[ds->l+i]*d[ld+ds->l+i];
798: }
799: H[n-1+(n-1)*n] = s[ds->l+n-1]*d[ds->l+n-1];
800: for (i=0;i<ds->k-ds->l;i++) {
801: H[ds->k-ds->l+i*n] = s[ds->k]*d[2*ld+ds->l+i];
802: H[i+(ds->k-ds->l)*n] = s[i+ds->l]*d[2*ld+ds->l+i];
803: }
804: } else {
805: for (j=0;j<n;j++) {
806: for (i=0;i<n;i++) H[i+j*n] = B[off+i+i*ld]*A[off+i+j*ld];
807: }
808: }
810: /* Compute eigenpairs */
811: PetscBLASIntCast(lw-nwu,&lwork);
812: DSAllocateMat_Private(ds,DS_MAT_X);
813: X = ds->mat[DS_MAT_X];
814: #if !defined(PETSC_USE_COMPLEX)
815: PetscStackCallBLAS("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));
816: #else
817: PetscStackCallBLAS("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));
819: /* Sort to have consecutive conjugate pairs
820: Separate real and imaginary part of complex eigenvectors*/
821: for (i=ds->l;i<ds->n;i++) {
822: j=i+1;
823: while (j<ds->n && (PetscAbsScalar(wr[i]-PetscConj(wr[j]))>PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON)) j++;
824: if (j==ds->n) {
825: if (PetscAbsReal(PetscImaginaryPart(wr[i]))<PetscAbsScalar(wr[i])*PETSC_SQRT_MACHINE_EPSILON) {
826: wr[i]=PetscRealPart(wr[i]); /* real eigenvalue */
827: for (k=ds->l;k<ds->n;k++) {
828: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
829: }
830: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Found complex without conjugate pair");
831: } else { /* complex eigenvalue */
832: if (j!=i+1) {
833: wr[j] = wr[i+1];
834: PetscMemcpy(X+j*ds->ld,X+(i+1)*ds->ld,ds->ld*sizeof(PetscScalar));
835: }
836: if (PetscImaginaryPart(wr[i])<0) {
837: wr[i] = PetscConj(wr[i]);
838: for (k=ds->l;k<ds->n;k++) {
839: X[k+(i+1)*ds->ld] = -PetscImaginaryPart(X[k+i*ds->ld]);
840: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
841: }
842: } else {
843: for (k=ds->l;k<ds->n;k++) {
844: X[k+(i+1)*ds->ld] = PetscImaginaryPart(X[k+i*ds->ld]);
845: X[k+i*ds->ld] = PetscRealPart(X[k+i*ds->ld]);
846: }
847: }
848: wr[i+1] = PetscConj(wr[i]);
849: i++;
850: }
851: }
852: #endif
853: SlepcCheckLapackInfo("geevx",info);
855: /* Compute real s-orthonormal basis */
856: DSGHIEPOrthogEigenv(ds,DS_MAT_X,wr,wi,PETSC_FALSE);
858: /* Recover eigenvalues from diagonal */
859: DSGHIEPComplexEigs(ds,0,ds->l,wr,wi);
860: #if defined(PETSC_USE_COMPLEX)
861: if (wi) {
862: for (i=ds->l;i<ds->n;i++) wi[i] = 0.0;
863: }
864: #endif
865: return(0);
866: #endif
867: }
869: PetscErrorCode DSSynchronize_GHIEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
870: {
872: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
873: PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
876: if (ds->compact) kr = 4*ld;
877: else k = 2*(ds->n-l)*ld;
878: if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
879: if (eigr) k += (ds->n-l);
880: if (eigi) k += (ds->n-l);
881: DSAllocateWork_Private(ds,k+kr,0,0);
882: PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size);
883: PetscMPIIntCast(ds->n-l,&n);
884: PetscMPIIntCast(ld*(ds->n-l),&ldn);
885: PetscMPIIntCast(ld*3,&ld3);
886: PetscMPIIntCast(ld,&ld_);
887: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
888: if (!rank) {
889: if (ds->compact) {
890: MPI_Pack(ds->rmat[DS_MAT_T],ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
891: MPI_Pack(ds->rmat[DS_MAT_D],ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
892: } else {
893: MPI_Pack(ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
894: MPI_Pack(ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
895: }
896: if (ds->state>DS_STATE_RAW) {
897: MPI_Pack(ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
898: }
899: if (eigr) {
900: MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
901: }
902: if (eigi) {
903: MPI_Pack(eigi+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
904: }
905: }
906: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
907: if (rank) {
908: if (ds->compact) {
909: MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_T],ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds));
910: MPI_Unpack(ds->work,size,&off,ds->rmat[DS_MAT_D],ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds));
911: } else {
912: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_A]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
913: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_B]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
914: }
915: if (ds->state>DS_STATE_RAW) {
916: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Q]+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
917: }
918: if (eigr) {
919: MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
920: }
921: if (eigi) {
922: MPI_Unpack(ds->work,size,&off,eigi+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
923: }
924: }
925: return(0);
926: }
928: PetscErrorCode DSHermitian_GHIEP(DS ds,DSMatType m,PetscBool *flg)
929: {
931: if (m==DS_MAT_A || m==DS_MAT_B) *flg = PETSC_TRUE;
932: else *flg = PETSC_FALSE;
933: return(0);
934: }
936: SLEPC_EXTERN PetscErrorCode DSCreate_GHIEP(DS ds)
937: {
939: ds->ops->allocate = DSAllocate_GHIEP;
940: ds->ops->view = DSView_GHIEP;
941: ds->ops->vectors = DSVectors_GHIEP;
942: ds->ops->solve[0] = DSSolve_GHIEP_QR_II;
943: ds->ops->solve[1] = DSSolve_GHIEP_HZ;
944: ds->ops->solve[2] = DSSolve_GHIEP_QR;
945: ds->ops->sort = DSSort_GHIEP;
946: ds->ops->synchronize = DSSynchronize_GHIEP;
947: ds->ops->hermitian = DSHermitian_GHIEP;
948: return(0);
949: }