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