Actual source code: dshep.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: static PetscErrorCode DSAllocate_HEP(DS ds,PetscInt ld)
15: {
16: PetscFunctionBegin;
17: if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
18: PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q));
19: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
20: PetscCall(PetscFree(ds->perm));
21: PetscCall(PetscMalloc1(ld,&ds->perm));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: /* 0 l k n-1
26: -----------------------------------------
27: |* . . |
28: | * . . |
29: | * . . |
30: | * . . |
31: |. . . . o o |
32: | o o |
33: | o o |
34: | o o |
35: | o o |
36: | o o |
37: |. . . . o o o o o o o x |
38: | x x x |
39: | x x x |
40: | x x x |
41: | x x x |
42: | x x x |
43: | x x x |
44: | x x x |
45: | x x x|
46: | x x|
47: -----------------------------------------
48: */
50: static PetscErrorCode DSSwitchFormat_HEP(DS ds)
51: {
52: PetscReal *T;
53: PetscScalar *A;
54: PetscInt i,n=ds->n,k=ds->k,ld=ds->ld;
56: PetscFunctionBegin;
57: /* switch from compact (arrow) to dense storage */
58: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
59: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
60: PetscCall(PetscArrayzero(A,ld*ld));
61: for (i=0;i<k;i++) {
62: A[i+i*ld] = T[i];
63: A[k+i*ld] = T[i+ld];
64: A[i+k*ld] = T[i+ld];
65: }
66: A[k+k*ld] = T[k];
67: for (i=k+1;i<n;i++) {
68: A[i+i*ld] = T[i];
69: A[i-1+i*ld] = T[i-1+ld];
70: A[i+(i-1)*ld] = T[i-1+ld];
71: }
72: if (ds->extrarow) A[n+(n-1)*ld] = T[n-1+ld];
73: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
74: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode DSView_HEP(DS ds,PetscViewer viewer)
79: {
80: PetscViewerFormat format;
81: PetscInt i,j,r,c,rows;
82: PetscReal *T,value;
83: const char *methodname[] = {
84: "Implicit QR method (_steqr)",
85: "Relatively Robust Representations (_stevr)",
86: "Divide and Conquer method (_stedc)",
87: "Block Divide and Conquer method (dsbtdc)"
88: };
89: const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
91: PetscFunctionBegin;
92: PetscCall(PetscViewerGetFormat(viewer,&format));
93: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
94: if (ds->bs>1) PetscCall(PetscViewerASCIIPrintf(viewer,"block size: %" PetscInt_FMT "\n",ds->bs));
95: if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
98: if (ds->compact) {
99: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
100: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
101: rows = ds->extrarow? ds->n+1: ds->n;
102: if (format == PETSC_VIEWER_ASCII_MATLAB) {
103: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,ds->n));
104: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
105: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
106: for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
107: for (i=0;i<rows-1;i++) {
108: r = PetscMax(i+2,ds->k+1);
109: c = i+1;
110: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",r,c,(double)T[i+ds->ld]));
111: if (i<ds->n-1 && ds->k<ds->n) { /* do not print vertical arrow when k=n */
112: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",c,r,(double)T[i+ds->ld]));
113: }
114: }
115: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
116: } else {
117: for (i=0;i<rows;i++) {
118: for (j=0;j<ds->n;j++) {
119: if (i==j) value = T[i];
120: else if ((i<ds->k && j==ds->k) || (i==ds->k && j<ds->k)) value = T[PetscMin(i,j)+ds->ld];
121: else if (i==j+1 && i>ds->k) value = T[i-1+ds->ld];
122: else if (i+1==j && j>ds->k) value = T[j-1+ds->ld];
123: else value = 0.0;
124: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
125: }
126: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
127: }
128: }
129: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
130: PetscCall(PetscViewerFlush(viewer));
131: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
132: } else PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
133: if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_Q));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode DSVectors_HEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
138: {
139: PetscScalar *Z;
140: const PetscScalar *Q;
141: PetscInt ld = ds->ld;
143: PetscFunctionBegin;
144: switch (mat) {
145: case DS_MAT_X:
146: case DS_MAT_Y:
147: if (j) {
148: PetscCall(MatDenseGetArray(ds->omat[mat],&Z));
149: if (ds->state>=DS_STATE_CONDENSED) {
150: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
151: PetscCall(PetscArraycpy(Z+(*j)*ld,Q+(*j)*ld,ld));
152: if (rnorm) *rnorm = PetscAbsScalar(Q[ds->n-1+(*j)*ld]);
153: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
154: } else {
155: PetscCall(PetscArrayzero(Z+(*j)*ld,ld));
156: Z[(*j)+(*j)*ld] = 1.0;
157: if (rnorm) *rnorm = 0.0;
158: }
159: PetscCall(MatDenseRestoreArray(ds->omat[mat],&Z));
160: } else {
161: if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatCopy(ds->omat[DS_MAT_Q],ds->omat[mat],SAME_NONZERO_PATTERN));
162: else PetscCall(DSSetIdentity(ds,mat));
163: }
164: break;
165: case DS_MAT_U:
166: case DS_MAT_V:
167: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
168: default:
169: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
170: }
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*
175: ARROWTRIDIAG reduces a symmetric arrowhead matrix of the form
177: [ d 0 0 0 e ]
178: [ 0 d 0 0 e ]
179: A = [ 0 0 d 0 e ]
180: [ 0 0 0 d e ]
181: [ e e e e d ]
183: to tridiagonal form
185: [ d e 0 0 0 ]
186: [ e d e 0 0 ]
187: T = Q'*A*Q = [ 0 e d e 0 ],
188: [ 0 0 e d e ]
189: [ 0 0 0 e d ]
191: where Q is an orthogonal matrix. Rutishauser's algorithm is used to
192: perform the reduction, which requires O(n**2) flops. The accumulation
193: of the orthogonal factor Q, however, requires O(n**3) flops.
195: Arguments
196: =========
198: N (input) INTEGER
199: The order of the matrix A. N >= 0.
201: D (input/output) DOUBLE PRECISION array, dimension (N)
202: On entry, the diagonal entries of the matrix A to be
203: reduced.
204: On exit, the diagonal entries of the reduced matrix T.
206: E (input/output) DOUBLE PRECISION array, dimension (N-1)
207: On entry, the off-diagonal entries of the matrix A to be
208: reduced.
209: On exit, the subdiagonal entries of the reduced matrix T.
211: Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
212: On exit, the orthogonal matrix Q.
214: LDQ (input) INTEGER
215: The leading dimension of the array Q.
217: Note
218: ====
219: Based on Fortran code contributed by Daniel Kressner
220: */
221: PetscErrorCode DSArrowTridiag(PetscBLASInt n,PetscReal *d,PetscReal *e,PetscScalar *Q,PetscBLASInt ld)
222: {
223: PetscBLASInt i,j,j2,one=1;
224: PetscReal c,s,p,off,temp;
226: PetscFunctionBegin;
227: if (n<=2) PetscFunctionReturn(PETSC_SUCCESS);
229: for (j=0;j<n-2;j++) {
231: /* Eliminate entry e(j) by a rotation in the planes (j,j+1) */
232: temp = e[j+1];
233: PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp,&e[j],&c,&s,&e[j+1]));
234: s = -s;
236: /* Apply rotation to diagonal elements */
237: temp = d[j+1];
238: e[j] = c*s*(temp-d[j]);
239: d[j+1] = s*s*d[j] + c*c*temp;
240: d[j] = c*c*d[j] + s*s*temp;
242: /* Apply rotation to Q */
243: j2 = j+2;
244: PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+j*ld,&one,Q+(j+1)*ld,&one,&c,&s));
246: /* Chase newly introduced off-diagonal entry to the top left corner */
247: for (i=j-1;i>=0;i--) {
248: off = -s*e[i];
249: e[i] = c*e[i];
250: temp = e[i+1];
251: PetscCallBLAS("LAPACKlartg",LAPACKREALlartg_(&temp,&off,&c,&s,&e[i+1]));
252: s = -s;
253: temp = (d[i]-d[i+1])*s - 2.0*c*e[i];
254: p = s*temp;
255: d[i+1] += p;
256: d[i] -= p;
257: e[i] = -e[i] - c*temp;
258: PetscCallBLAS("BLASrot",BLASMIXEDrot_(&j2,Q+i*ld,&one,Q+(i+1)*ld,&one,&c,&s));
259: }
260: }
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*
265: Reduce to tridiagonal form by means of DSArrowTridiag.
266: */
267: static PetscErrorCode DSIntermediate_HEP(DS ds)
268: {
269: PetscInt i;
270: PetscBLASInt n1 = 0,n2,lwork,l = 0,n = 0,ld,off;
271: PetscScalar *Q,*work,*tau;
272: const PetscScalar *A;
273: PetscReal *d,*e;
274: Mat At,Qt; /* trailing submatrices */
276: PetscFunctionBegin;
277: PetscCall(PetscBLASIntCast(ds->n,&n));
278: PetscCall(PetscBLASIntCast(ds->l,&l));
279: PetscCall(PetscBLASIntCast(ds->ld,&ld));
280: PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-l+1),&n1)); /* size of leading block, excl. locked */
281: n2 = n-l; /* n2 = n1 + size of trailing block */
282: off = l+l*ld;
283: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
284: e = d+ld;
285: PetscCall(DSSetIdentity(ds,DS_MAT_Q));
286: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
288: if (ds->compact) {
290: if (ds->state<DS_STATE_INTERMEDIATE) PetscCall(DSArrowTridiag(n1,d+l,e+l,Q+off,ld));
292: } else {
294: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&A));
295: for (i=0;i<l;i++) { d[i] = PetscRealPart(A[i+i*ld]); e[i] = 0.0; }
297: if (ds->state<DS_STATE_INTERMEDIATE) {
298: PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_A],ds->l,ds->n,ds->l,ds->n,&At));
299: PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_Q],ds->l,ds->n,ds->l,ds->n,&Qt));
300: PetscCall(MatCopy(At,Qt,SAME_NONZERO_PATTERN));
301: PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_A],&At));
302: PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_Q],&Qt));
303: PetscCall(DSAllocateWork_Private(ds,ld+ld*ld,0,0));
304: tau = ds->work;
305: work = ds->work+ld;
306: lwork = ld*ld;
307: PetscCallLAPACKInfo("LAPACKsytrd",LAPACKsytrd_("L",&n2,Q+off,&ld,d+l,e+l,tau,work,&lwork,&info));
308: PetscCallLAPACKInfo("LAPACKorgtr",LAPACKorgtr_("L",&n2,Q+off,&ld,tau,work,&lwork,&info));
309: } else {
310: /* copy tridiagonal to d,e */
311: for (i=l;i<n;i++) d[i] = PetscRealPart(A[i+i*ld]);
312: for (i=l;i<n-1;i++) e[i] = PetscRealPart(A[(i+1)+i*ld]);
313: }
314: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&A));
315: }
316: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
317: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode DSSort_HEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
322: {
323: PetscInt n,l,i,*perm,ld=ds->ld;
324: PetscScalar *A;
325: PetscReal *d;
327: PetscFunctionBegin;
328: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
329: n = ds->n;
330: l = ds->l;
331: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
332: perm = ds->perm;
333: if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
334: else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
335: for (i=l;i<n;i++) wr[i] = d[perm[i]];
336: PetscCall(DSPermuteColumns_Private(ds,l,n,n,DS_MAT_Q,perm));
337: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
338: if (!ds->compact) {
339: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
340: for (i=l;i<n;i++) A[i+i*ld] = wr[i];
341: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
342: }
343: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: static PetscErrorCode DSUpdateExtraRow_HEP(DS ds)
348: {
349: PetscInt i;
350: PetscBLASInt n,ld,incx=1;
351: PetscScalar *A,*x,*y,one=1.0,zero=0.0;
352: PetscReal *T,*e,beta;
353: const PetscScalar *Q;
355: PetscFunctionBegin;
356: PetscCall(PetscBLASIntCast(ds->n,&n));
357: PetscCall(PetscBLASIntCast(ds->ld,&ld));
358: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_Q],&Q));
359: if (ds->compact) {
360: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
361: e = T+ld;
362: beta = e[n-1]; /* in compact, we assume all entries are zero except the last one */
363: for (i=0;i<n;i++) e[i] = PetscRealPart(beta*Q[n-1+i*ld]);
364: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
365: ds->k = n;
366: } else {
367: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
368: PetscCall(DSAllocateWork_Private(ds,2*ld,0,0));
369: x = ds->work;
370: y = ds->work+ld;
371: for (i=0;i<n;i++) x[i] = PetscConj(A[n+i*ld]);
372: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&n,&one,Q,&ld,x,&incx,&zero,y,&incx));
373: for (i=0;i<n;i++) A[n+i*ld] = PetscConj(y[i]);
374: ds->k = n;
375: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
376: }
377: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_Q],&Q));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: static PetscErrorCode DSSolve_HEP_QR(DS ds,PetscScalar *wr,PetscScalar *wi)
382: {
383: PetscInt i;
384: PetscBLASInt n1,l = 0,n = 0,ld,off;
385: PetscScalar *Q,*A;
386: PetscReal *d,*e;
388: PetscFunctionBegin;
389: PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
390: PetscCall(PetscBLASIntCast(ds->n,&n));
391: PetscCall(PetscBLASIntCast(ds->l,&l));
392: PetscCall(PetscBLASIntCast(ds->ld,&ld));
393: n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
394: off = l+l*ld;
395: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
396: e = d+ld;
398: /* Reduce to tridiagonal form */
399: PetscCall(DSIntermediate_HEP(ds));
401: /* Solve the tridiagonal eigenproblem */
402: for (i=0;i<l;i++) wr[i] = d[i];
404: PetscCall(DSAllocateWork_Private(ds,0,2*ld,0));
405: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
406: PetscCallLAPACKInfo("LAPACKsteqr",LAPACKsteqr_("V",&n1,d+l,e+l,Q+off,&ld,ds->rwork,&info));
407: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
408: for (i=l;i<n;i++) wr[i] = d[i];
410: /* Create diagonal matrix as a result */
411: if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
412: else {
413: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
414: for (i=l;i<n;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
415: for (i=l;i<n;i++) A[i+i*ld] = d[i];
416: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
417: }
418: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
420: /* Set zero wi */
421: if (wi) for (i=l;i<n;i++) wi[i] = 0.0;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: static PetscErrorCode DSSolve_HEP_MRRR(DS ds,PetscScalar *wr,PetscScalar *wi)
426: {
427: Mat At,Qt; /* trailing submatrices */
428: PetscInt i;
429: PetscBLASInt n1 = 0,n2 = 0,n3,lrwork,liwork,l = 0,n = 0,m = 0,ld,off,il,iu,*isuppz;
430: PetscScalar *A,*Q,*W=NULL,one=1.0,zero=0.0;
431: PetscReal *d,*e,abstol=0.0,vl,vu;
432: #if PetscDefined(USE_COMPLEX)
433: PetscInt j;
434: PetscReal *Qr,*ritz;
435: #endif
437: PetscFunctionBegin;
438: PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
439: PetscCall(PetscBLASIntCast(ds->n,&n));
440: PetscCall(PetscBLASIntCast(ds->l,&l));
441: PetscCall(PetscBLASIntCast(ds->ld,&ld));
442: PetscCall(PetscBLASIntCast(ds->k-l+1,&n1)); /* size of leading block, excl. locked */
443: PetscCall(PetscBLASIntCast(n-ds->k-1,&n2)); /* size of trailing block */
444: n3 = n1+n2;
445: off = l+l*ld;
446: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
447: e = d+ld;
449: /* Reduce to tridiagonal form */
450: PetscCall(DSIntermediate_HEP(ds));
452: /* Solve the tridiagonal eigenproblem */
453: for (i=0;i<l;i++) wr[i] = d[i];
455: if (ds->state<DS_STATE_INTERMEDIATE) { /* Q contains useful info */
456: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
457: PetscCall(MatCopy(ds->omat[DS_MAT_Q],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
458: }
459: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
460: lrwork = 20*ld;
461: liwork = 10*ld;
462: #if PetscDefined(USE_COMPLEX)
463: PetscCall(DSAllocateWork_Private(ds,0,lrwork+ld+ld*ld,liwork+2*ld));
464: #else
465: PetscCall(DSAllocateWork_Private(ds,0,lrwork+ld,liwork+2*ld));
466: #endif
467: isuppz = ds->iwork+liwork;
468: #if PetscDefined(USE_COMPLEX)
469: ritz = ds->rwork+lrwork;
470: Qr = ds->rwork+lrwork+ld;
471: PetscCallLAPACKInfo("LAPACKstevr",LAPACKstevr_("V","A",&n3,d+l,e+l,&vl,&vu,&il,&iu,&abstol,&m,ritz+l,Qr+off,&ld,isuppz,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
472: for (i=l;i<n;i++) wr[i] = ritz[i];
473: #else
474: PetscCallLAPACKInfo("LAPACKstevr",LAPACKstevr_("V","A",&n3,d+l,e+l,&vl,&vu,&il,&iu,&abstol,&m,wr+l,Q+off,&ld,isuppz,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
475: #endif
476: #if PetscDefined(USE_COMPLEX)
477: for (i=l;i<n;i++)
478: for (j=l;j<n;j++)
479: Q[i+j*ld] = Qr[i+j*ld];
480: #endif
481: if (ds->state<DS_STATE_INTERMEDIATE) { /* accumulate previous Q */
482: if (ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
483: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
484: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
485: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n3,&n3,&n3,&one,W+off,&ld,Q+off,&ld,&zero,A+off,&ld));
486: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
487: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
488: PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_A],ds->l,ds->n,ds->l,ds->n,&At));
489: PetscCall(MatDenseGetSubMatrix(ds->omat[DS_MAT_Q],ds->l,ds->n,ds->l,ds->n,&Qt));
490: PetscCall(MatCopy(At,Qt,SAME_NONZERO_PATTERN));
491: PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_A],&At));
492: PetscCall(MatDenseRestoreSubMatrix(ds->omat[DS_MAT_Q],&Qt));
493: }
494: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
495: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
497: /* Create diagonal matrix as a result */
498: if (ds->compact) PetscCall(PetscArrayzero(e,n-1));
499: else {
500: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
501: for (i=l;i<n;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
502: for (i=l;i<n;i++) A[i+i*ld] = d[i];
503: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
504: }
505: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
507: /* Set zero wi */
508: if (wi) for (i=l;i<n;i++) wi[i] = 0.0;
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode DSSolve_HEP_DC(DS ds,PetscScalar *wr,PetscScalar *wi)
513: {
514: PetscInt i;
515: PetscBLASInt n1,l = 0,ld,off,lrwork,liwork;
516: PetscScalar *Q,*A;
517: PetscReal *d,*e;
518: #if PetscDefined(USE_COMPLEX)
519: PetscBLASInt lwork;
520: PetscInt j;
521: #endif
523: PetscFunctionBegin;
524: PetscCheck(ds->bs==1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for bs>1");
525: PetscCall(PetscBLASIntCast(ds->l,&l));
526: PetscCall(PetscBLASIntCast(ds->ld,&ld));
527: PetscCall(PetscBLASIntCast(ds->n-ds->l,&n1));
528: off = l+l*ld;
529: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
530: e = d+ld;
532: /* Reduce to tridiagonal form */
533: PetscCall(DSIntermediate_HEP(ds));
535: /* Solve the tridiagonal eigenproblem */
536: for (i=0;i<l;i++) wr[i] = d[i];
538: lrwork = 5*n1*n1+3*n1+1;
539: liwork = 5*n1*n1+6*n1+6;
540: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
541: #if !PetscDefined(USE_COMPLEX)
542: PetscCall(DSAllocateWork_Private(ds,0,lrwork,liwork));
543: PetscCallLAPACKInfo("LAPACKstedc",LAPACKstedc_("V",&n1,d+l,e+l,Q+off,&ld,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
544: #else
545: lwork = ld*ld;
546: PetscCall(DSAllocateWork_Private(ds,lwork,lrwork,liwork));
547: PetscCallLAPACKInfo("LAPACKstedc",LAPACKstedc_("V",&n1,d+l,e+l,Q+off,&ld,ds->work,&lwork,ds->rwork,&lrwork,ds->iwork,&liwork,&info));
548: /* Fixing Lapack bug*/
549: for (j=ds->l;j<ds->n;j++)
550: for (i=0;i<ds->l;i++) Q[i+j*ld] = 0.0;
551: #endif
552: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
553: for (i=l;i<ds->n;i++) wr[i] = d[i];
555: /* Create diagonal matrix as a result */
556: if (ds->compact) PetscCall(PetscArrayzero(e,ds->n-1));
557: else {
558: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
559: for (i=l;i<ds->n;i++) PetscCall(PetscArrayzero(A+l+i*ld,ds->n-l));
560: for (i=l;i<ds->n;i++) A[i+i*ld] = d[i];
561: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
562: }
563: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
565: /* Set zero wi */
566: if (wi) for (i=l;i<ds->n;i++) wi[i] = 0.0;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: #if !PetscDefined(USE_COMPLEX)
571: static PetscErrorCode DSSolve_HEP_BDC(DS ds,PetscScalar *wr,PetscScalar *wi)
572: {
573: PetscBLASInt i,j,k,m,n = 0,info,nblks,bs = 0,ld = 0,lde,lrwork,liwork,*ksizes,*iwork,mingapi;
574: PetscScalar *Q,*A;
575: PetscReal *D,*E,*d,*e,tol=PETSC_MACHINE_EPSILON/2,tau1=1e-16,tau2=1e-18,*rwork,mingap;
577: PetscFunctionBegin;
578: PetscCheck(ds->l==0,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"This method is not prepared for l>1");
579: PetscCheck(!ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for compact storage");
580: PetscCall(PetscBLASIntCast(ds->ld,&ld));
581: PetscCall(PetscBLASIntCast(ds->bs,&bs));
582: PetscCall(PetscBLASIntCast(ds->n,&n));
583: nblks = n/bs;
584: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
585: e = d+ld;
586: lrwork = 4*n*n+60*n+1;
587: liwork = 5*n+5*nblks-1;
588: lde = 2*bs+1;
589: PetscCall(DSAllocateWork_Private(ds,bs*n+lde*lde*(nblks-1),lrwork,nblks+liwork));
590: D = ds->work;
591: E = ds->work+bs*n;
592: rwork = ds->rwork;
593: ksizes = ds->iwork;
594: iwork = ds->iwork+nblks;
595: PetscCall(PetscArrayzero(iwork,liwork));
597: /* Copy matrix to block tridiagonal format */
598: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
599: j=0;
600: for (i=0;i<nblks;i++) {
601: ksizes[i]=bs;
602: for (k=0;k<bs;k++)
603: for (m=0;m<bs;m++)
604: D[k+m*bs+i*bs*bs] = PetscRealPart(A[j+k+(j+m)*n]);
605: j = j + bs;
606: }
607: j=0;
608: for (i=0;i<nblks-1;i++) {
609: for (k=0;k<bs;k++)
610: for (m=0;m<bs;m++)
611: E[k+m*lde+i*lde*lde] = PetscRealPart(A[j+bs+k+(j+m)*n]);
612: j = j + bs;
613: }
614: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
616: /* Solve the block tridiagonal eigenproblem */
617: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
618: PetscCall(BDC_dsbtdc_("D","A",n,nblks,ksizes,D,bs,bs,E,lde,lde,tol,tau1,tau2,d,Q,n,rwork,lrwork,iwork,liwork,&mingap,&mingapi,&info,1,1));
619: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
620: for (i=0;i<ds->n;i++) wr[i] = d[i];
622: /* Create diagonal matrix as a result */
623: if (ds->compact) PetscCall(PetscArrayzero(e,ds->n-1));
624: else {
625: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
626: for (i=0;i<ds->n;i++) PetscCall(PetscArrayzero(A+i*ld,ds->n));
627: for (i=0;i<ds->n;i++) A[i+i*ld] = wr[i];
628: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
629: }
630: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
632: /* Set zero wi */
633: if (wi) for (i=0;i<ds->n;i++) wi[i] = 0.0;
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
636: #endif
638: static PetscErrorCode DSTruncate_HEP(DS ds,PetscInt n,PetscBool trim)
639: {
640: PetscInt i,ld=ds->ld,l=ds->l;
641: PetscScalar *A;
643: PetscFunctionBegin;
644: if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
645: if (trim) {
646: if (!ds->compact && ds->extrarow) { /* clean extra row */
647: for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
648: }
649: ds->l = 0;
650: ds->k = 0;
651: ds->n = n;
652: ds->t = ds->n; /* truncated length equal to the new dimension */
653: } else {
654: if (!ds->compact && ds->extrarow && ds->k==ds->n) {
655: /* copy entries of extra row to the new position, then clean last row */
656: for (i=l;i<n;i++) A[n+i*ld] = A[ds->n+i*ld];
657: for (i=l;i<ds->n;i++) A[ds->n+i*ld] = 0.0;
658: }
659: ds->k = ds->extrarow? n: 0;
660: ds->t = ds->n; /* truncated length equal to previous dimension */
661: ds->n = n;
662: }
663: if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: #if !PetscDefined(HAVE_MPIUNI)
668: static PetscErrorCode DSSynchronize_HEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
669: {
670: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
671: PetscMPIInt n,rank,off=0,size,ldn,ld3;
672: PetscScalar *A,*Q;
673: PetscReal *T;
675: PetscFunctionBegin;
676: if (ds->compact) kr = 3*ld;
677: else k = (ds->n-l)*ld;
678: if (ds->state>DS_STATE_RAW) k += (ds->n-l)*ld;
679: if (eigr) k += (ds->n-l);
680: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
681: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
682: PetscCall(PetscMPIIntCast(ds->n-l,&n));
683: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
684: PetscCall(PetscMPIIntCast(ld*3,&ld3));
685: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
686: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
687: if (ds->state>DS_STATE_RAW) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
688: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
689: if (!rank) {
690: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
691: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
692: if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Pack(Q+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
693: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
694: }
695: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
696: if (rank) {
697: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
698: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
699: if (ds->state>DS_STATE_RAW) PetscCallMPI(MPI_Unpack(ds->work,size,&off,Q+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
700: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
701: }
702: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
703: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
704: if (ds->state>DS_STATE_RAW) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
707: #endif
709: static PetscErrorCode DSCond_HEP(DS ds,PetscReal *cond)
710: {
711: PetscScalar *work;
712: PetscReal *rwork;
713: PetscBLASInt *ipiv;
714: PetscBLASInt lwork,n,ld;
715: PetscReal hn,hin;
716: PetscScalar *A;
718: PetscFunctionBegin;
719: PetscCall(PetscBLASIntCast(ds->n,&n));
720: PetscCall(PetscBLASIntCast(ds->ld,&ld));
721: lwork = 8*ld;
722: PetscCall(DSAllocateWork_Private(ds,lwork,ld,ld));
723: work = ds->work;
724: rwork = ds->rwork;
725: ipiv = ds->iwork;
726: if (ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
727: PetscCall(DSSwitchFormat_HEP(ds));
729: /* use workspace matrix W to avoid overwriting A */
730: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
731: PetscCall(MatCopy(ds->omat[DS_MAT_A],ds->omat[DS_MAT_W],SAME_NONZERO_PATTERN));
732: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&A));
734: /* norm of A */
735: hn = LAPACKlange_("I",&n,&n,A,&ld,rwork);
737: /* norm of inv(A) */
738: PetscCallLAPACKInfo("LAPACKgetrf",LAPACKgetrf_(&n,&n,A,&ld,ipiv,&info));
739: PetscCallLAPACKInfo("LAPACKgetri",LAPACKgetri_(&n,A,&ld,ipiv,work,&lwork,&info));
740: hin = LAPACKlange_("I",&n,&n,A,&ld,rwork);
741: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&A));
743: *cond = hn*hin;
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: static PetscErrorCode DSTranslateRKS_HEP(DS ds,PetscScalar alpha)
748: {
749: PetscInt i,j,k=ds->k;
750: PetscScalar *Q,*A,*R,*tau,*work;
751: PetscBLASInt ld,n1,n0,lwork;
753: PetscFunctionBegin;
754: PetscCall(PetscBLASIntCast(ds->ld,&ld));
755: PetscCall(DSAllocateWork_Private(ds,ld*ld,0,0));
756: tau = ds->work;
757: work = ds->work+ld;
758: PetscCall(PetscBLASIntCast(ld*(ld-1),&lwork));
759: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
760: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
761: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_Q],&Q));
762: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_W],&R));
764: /* copy I+alpha*A */
765: PetscCall(PetscArrayzero(Q,ld*ld));
766: PetscCall(PetscArrayzero(R,ld*ld));
767: for (i=0;i<k;i++) {
768: Q[i+i*ld] = 1.0 + alpha*A[i+i*ld];
769: Q[k+i*ld] = alpha*A[k+i*ld];
770: }
772: /* compute qr */
773: PetscCall(PetscBLASIntCast(k+1,&n1));
774: PetscCall(PetscBLASIntCast(k,&n0));
775: PetscCallLAPACKInfo("LAPACKgeqrf",LAPACKgeqrf_(&n1,&n0,Q,&ld,tau,work,&lwork,&info));
777: /* copy R from Q */
778: for (j=0;j<k;j++)
779: for (i=0;i<=j;i++)
780: R[i+j*ld] = Q[i+j*ld];
782: /* compute orthogonal matrix in Q */
783: PetscCallLAPACKInfo("LAPACKorgqr",LAPACKorgqr_(&n1,&n1,&n0,Q,&ld,tau,work,&lwork,&info));
785: /* compute the updated matrix of projected problem */
786: for (j=0;j<k;j++)
787: for (i=0;i<k+1;i++)
788: A[j*ld+i] = Q[i*ld+j];
789: alpha = -1.0/alpha;
790: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n1,&n0,&alpha,R,&ld,A,&ld));
791: for (i=0;i<k;i++)
792: A[ld*i+i] -= alpha;
794: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
795: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_Q],&Q));
796: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_W],&R));
797: PetscFunctionReturn(PETSC_SUCCESS);
798: }
800: static PetscErrorCode DSHermitian_HEP(DS ds,DSMatType m,PetscBool *flg)
801: {
802: PetscFunctionBegin;
803: if (m==DS_MAT_A && !ds->extrarow) *flg = PETSC_TRUE;
804: else *flg = PETSC_FALSE;
805: PetscFunctionReturn(PETSC_SUCCESS);
806: }
808: static PetscErrorCode DSSetCompact_HEP(DS ds,PetscBool comp)
809: {
810: PetscFunctionBegin;
811: if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
812: PetscFunctionReturn(PETSC_SUCCESS);
813: }
815: static PetscErrorCode DSReallocate_HEP(DS ds,PetscInt ld)
816: {
817: PetscInt i,*perm=ds->perm;
819: PetscFunctionBegin;
820: for (i=0;i<DS_NUM_MAT;i++) {
821: if (!ds->compact && i==DS_MAT_A) continue;
822: if (i!=DS_MAT_Q && i!=DS_MAT_T) PetscCall(MatDestroy(&ds->omat[i]));
823: }
825: if (!ds->compact) PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
826: PetscCall(DSReallocateMat_Private(ds,DS_MAT_Q,ld));
827: PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
829: PetscCall(PetscMalloc1(ld,&ds->perm));
830: PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
831: PetscCall(PetscFree(perm));
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: /*MC
836: DSHEP - Dense Hermitian Eigenvalue Problem.
838: Notes:
839: The problem is expressed as $AX = X\Lambda$, where $A$ is real symmetric
840: (or complex Hermitian). $\Lambda$ is a diagonal matrix whose diagonal
841: elements are the arguments of `DSSolve()`. After solve, $A$ is overwritten
842: with $\Lambda$.
844: In the intermediate state $A$ is reduced to tridiagonal form. In compact
845: storage format, the symmetric tridiagonal matrix is stored in $T$.
847: Used DS matrices:
848: + `DS_MAT_A` - problem matrix (used only if `compact=PETSC_FALSE`)
849: . `DS_MAT_T` - symmetric tridiagonal matrix
850: - `DS_MAT_Q` - orthogonal/unitary transformation that reduces to tridiagonal form
851: (intermediate step) or matrix of orthogonal eigenvectors, which is equal to $X$
853: Implemented methods:
854: + 0 - Implicit QR (`_steqr`)
855: . 1 - Multiple Relatively Robust Representations (`_stevr`)
856: . 2 - Divide and Conquer (`_stedc`)
857: - 3 - Block Divide and Conquer (real scalars only)
859: Level: beginner
861: .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSSetCompact()`
862: M*/
863: SLEPC_EXTERN PetscErrorCode DSCreate_HEP(DS ds)
864: {
865: PetscFunctionBegin;
866: ds->ops->allocate = DSAllocate_HEP;
867: ds->ops->view = DSView_HEP;
868: ds->ops->vectors = DSVectors_HEP;
869: ds->ops->solve[0] = DSSolve_HEP_QR;
870: ds->ops->solve[1] = DSSolve_HEP_MRRR;
871: ds->ops->solve[2] = DSSolve_HEP_DC;
872: #if !PetscDefined(USE_COMPLEX)
873: ds->ops->solve[3] = DSSolve_HEP_BDC;
874: #endif
875: ds->ops->sort = DSSort_HEP;
876: ds->ops->truncate = DSTruncate_HEP;
877: ds->ops->update = DSUpdateExtraRow_HEP;
878: ds->ops->cond = DSCond_HEP;
879: ds->ops->transrks = DSTranslateRKS_HEP;
880: ds->ops->hermitian = DSHermitian_HEP;
881: #if !PetscDefined(HAVE_MPIUNI)
882: ds->ops->synchronize = DSSynchronize_HEP;
883: #endif
884: ds->ops->setcompact = DSSetCompact_HEP;
885: ds->ops->reallocate = DSReallocate_HEP;
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }