Actual source code: dspriv.c
slepc-3.17.0 2022-03-31
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: */
10: /*
11: Private DS routines
12: */
14: #include <slepc/private/dsimpl.h>
15: #include <slepcblaslapack.h>
17: PetscErrorCode DSAllocateMatrix_Private(DS ds,DSMatType m,PetscBool isreal)
18: {
19: size_t sz;
20: PetscInt n,d,nelem;
21: PetscBool ispep,isnep;
23: PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep);
24: PetscObjectTypeCompare((PetscObject)ds,DSNEP,&isnep);
25: if (ispep) DSPEPGetDegree(ds,&d);
26: if (isnep) DSNEPGetMinimality(ds,&d);
27: if ((ispep || isnep) && (m==DS_MAT_A || m==DS_MAT_B || m==DS_MAT_W || m==DS_MAT_U || m==DS_MAT_X || m==DS_MAT_Y)) n = d*ds->ld;
28: else n = ds->ld;
30: switch (m) {
31: case DS_MAT_T:
32: nelem = 3*ds->ld;
33: break;
34: case DS_MAT_D:
35: nelem = ds->ld;
36: break;
37: case DS_MAT_X:
38: nelem = ds->ld*n;
39: break;
40: case DS_MAT_Y:
41: nelem = ds->ld*n;
42: break;
43: default:
44: nelem = n*n;
45: }
46: if (isreal) {
47: sz = nelem*sizeof(PetscReal);
48: if (ds->rmat[m]) PetscFree(ds->rmat[m]);
49: else PetscLogObjectMemory((PetscObject)ds,sz);
50: PetscCalloc1(nelem,&ds->rmat[m]);
51: } else {
52: sz = nelem*sizeof(PetscScalar);
53: if (ds->mat[m]) PetscFree(ds->mat[m]);
54: else PetscLogObjectMemory((PetscObject)ds,sz);
55: PetscCalloc1(nelem,&ds->mat[m]);
56: }
57: PetscFunctionReturn(0);
58: }
60: PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
61: {
62: if (s>ds->lwork) {
63: PetscFree(ds->work);
64: PetscMalloc1(s,&ds->work);
65: PetscLogObjectMemory((PetscObject)ds,(s-ds->lwork)*sizeof(PetscScalar));
66: ds->lwork = s;
67: }
68: if (r>ds->lrwork) {
69: PetscFree(ds->rwork);
70: PetscMalloc1(r,&ds->rwork);
71: PetscLogObjectMemory((PetscObject)ds,(r-ds->lrwork)*sizeof(PetscReal));
72: ds->lrwork = r;
73: }
74: if (i>ds->liwork) {
75: PetscFree(ds->iwork);
76: PetscMalloc1(i,&ds->iwork);
77: PetscLogObjectMemory((PetscObject)ds,(i-ds->liwork)*sizeof(PetscBLASInt));
78: ds->liwork = i;
79: }
80: PetscFunctionReturn(0);
81: }
83: /*@C
84: DSViewMat - Prints one of the internal DS matrices.
86: Collective on ds
88: Input Parameters:
89: + ds - the direct solver context
90: . viewer - visualization context
91: - m - matrix to display
93: Note:
94: Works only for ascii viewers. Set the viewer in Matlab format if
95: want to paste into Matlab.
97: Level: developer
99: .seealso: DSView()
100: @*/
101: PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
102: {
103: PetscInt i,j,rows,cols;
104: PetscScalar *v;
105: PetscViewerFormat format;
106: #if defined(PETSC_USE_COMPLEX)
107: PetscBool allreal = PETSC_TRUE;
108: #endif
112: DSCheckValidMat(ds,m,3);
113: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer);
116: PetscViewerGetFormat(viewer,&format);
117: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0);
118: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
119: DSMatGetSize(ds,m,&rows,&cols);
120: #if defined(PETSC_USE_COMPLEX)
121: /* determine if matrix has all real values */
122: v = ds->mat[m];
123: for (i=0;i<rows;i++)
124: for (j=0;j<cols;j++)
125: if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
126: #endif
127: if (format == PETSC_VIEWER_ASCII_MATLAB) {
128: PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols);
129: PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]);
130: } else PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]);
132: for (i=0;i<rows;i++) {
133: v = ds->mat[m]+i;
134: for (j=0;j<cols;j++) {
135: #if defined(PETSC_USE_COMPLEX)
136: if (allreal) PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
137: else PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
138: #else
139: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
140: #endif
141: v += ds->ld;
142: }
143: PetscViewerASCIIPrintf(viewer,"\n");
144: }
146: if (format == PETSC_VIEWER_ASCII_MATLAB) PetscViewerASCIIPrintf(viewer,"];\n");
147: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
148: PetscViewerFlush(viewer);
149: PetscFunctionReturn(0);
150: }
152: PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
153: {
154: PetscScalar re,im,wi0;
155: PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
157: n = ds->t; /* sort only first t pairs if truncated */
158: /* insertion sort */
159: i=ds->l+1;
160: #if !defined(PETSC_USE_COMPLEX)
161: if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
162: #else
163: if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
164: #endif
165: for (;i<n;i+=d) {
166: re = wr[perm[i]];
167: if (wi) im = wi[perm[i]];
168: else im = 0.0;
169: tmp1 = perm[i];
170: #if !defined(PETSC_USE_COMPLEX)
171: if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
172: else d = 1;
173: #else
174: if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
175: else d = 1;
176: #endif
177: j = i-1;
178: if (wi) wi0 = wi[perm[j]];
179: else wi0 = 0.0;
180: SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
181: while (result<0 && j>=ds->l) {
182: perm[j+d] = perm[j];
183: j--;
184: #if !defined(PETSC_USE_COMPLEX)
185: if (wi && wi[perm[j+1]]!=0)
186: #else
187: if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
188: #endif
189: { perm[j+d] = perm[j]; j--; }
191: if (j>=ds->l) {
192: if (wi) wi0 = wi[perm[j]];
193: else wi0 = 0.0;
194: SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
195: }
196: }
197: perm[j+1] = tmp1;
198: if (d==2) perm[j+2] = tmp2;
199: }
200: PetscFunctionReturn(0);
201: }
203: PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
204: {
205: PetscScalar re;
206: PetscInt i,j,result,tmp,l,n;
208: n = ds->t; /* sort only first t pairs if truncated */
209: l = ds->l;
210: /* insertion sort */
211: for (i=l+1;i<n;i++) {
212: re = eig[perm[i]];
213: j = i-1;
214: SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
215: while (result<0 && j>=l) {
216: tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
217: if (j>=l) SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
218: }
219: }
220: PetscFunctionReturn(0);
221: }
223: /*
224: DSCopyMatrix_Private - Copies the trailing block of a matrix (from
225: rows/columns l to n).
226: */
227: PetscErrorCode DSCopyMatrix_Private(DS ds,DSMatType dst,DSMatType src)
228: {
229: PetscInt j,m,off,ld;
230: PetscScalar *S,*D;
232: ld = ds->ld;
233: m = ds->n-ds->l;
234: off = ds->l+ds->l*ld;
235: S = ds->mat[src];
236: D = ds->mat[dst];
237: for (j=0;j<m;j++) PetscArraycpy(D+off+j*ld,S+off+j*ld,m);
238: PetscFunctionReturn(0);
239: }
241: /*
242: Permute comumns [istart..iend-1] of [mat] according to perm. Columns have length n
243: */
244: PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat,PetscInt *perm)
245: {
246: PetscInt i,j,k,p,ld;
247: PetscScalar *Q,rtmp;
249: ld = ds->ld;
250: Q = ds->mat[mat];
251: for (i=istart;i<iend;i++) {
252: p = perm[i];
253: if (p != i) {
254: j = i + 1;
255: while (perm[j] != i) j++;
256: perm[j] = p; perm[i] = i;
257: /* swap columns i and j */
258: for (k=0;k<n;k++) {
259: rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
260: }
261: }
262: }
263: PetscFunctionReturn(0);
264: }
266: /*
267: The same as DSPermuteColumns_Private but for two matrices [mat1] and [mat2]
268: */
269: PetscErrorCode DSPermuteColumnsTwo_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
270: {
271: PetscInt i,j,k,p,ld;
272: PetscScalar *Q,*Z,rtmp,rtmp2;
274: ld = ds->ld;
275: Q = ds->mat[mat1];
276: Z = ds->mat[mat2];
277: for (i=istart;i<iend;i++) {
278: p = perm[i];
279: if (p != i) {
280: j = i + 1;
281: while (perm[j] != i) j++;
282: perm[j] = p; perm[i] = i;
283: /* swap columns i and j */
284: for (k=0;k<n;k++) {
285: rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
286: rtmp2 = Z[k+p*ld]; Z[k+p*ld] = Z[k+i*ld]; Z[k+i*ld] = rtmp2;
287: }
288: }
289: }
290: PetscFunctionReturn(0);
291: }
293: /*
294: Permute rows [istart..iend-1] of [mat] according to perm. Rows have length m
295: */
296: PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt m,DSMatType mat,PetscInt *perm)
297: {
298: PetscInt i,j,k,p,ld;
299: PetscScalar *Q,rtmp;
301: ld = ds->ld;
302: Q = ds->mat[mat];
303: for (i=istart;i<iend;i++) {
304: p = perm[i];
305: if (p != i) {
306: j = i + 1;
307: while (perm[j] != i) j++;
308: perm[j] = p; perm[i] = i;
309: /* swap rows i and j */
310: for (k=0;k<m;k++) {
311: rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
312: }
313: }
314: }
315: PetscFunctionReturn(0);
316: }
318: /*
319: Permute columns [istart..iend-1] of [mat1] and [mat2] according to perm.
320: Columns of [mat1] have length n, columns of [mat2] have length m
321: */
322: PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt istart,PetscInt iend,PetscInt n,PetscInt m,DSMatType mat1,DSMatType mat2,PetscInt *perm)
323: {
324: PetscInt i,j,k,p,ld;
325: PetscScalar *U,*V,rtmp;
327: ld = ds->ld;
328: U = ds->mat[mat1];
329: V = ds->mat[mat2];
330: for (i=istart;i<iend;i++) {
331: p = perm[i];
332: if (p != i) {
333: j = i + 1;
334: while (perm[j] != i) j++;
335: perm[j] = p; perm[i] = i;
336: /* swap columns i and j of U */
337: for (k=0;k<n;k++) {
338: rtmp = U[k+p*ld]; U[k+p*ld] = U[k+i*ld]; U[k+i*ld] = rtmp;
339: }
340: /* swap columns i and j of V */
341: for (k=0;k<m;k++) {
342: rtmp = V[k+p*ld]; V[k+p*ld] = V[k+i*ld]; V[k+i*ld] = rtmp;
343: }
344: }
345: }
346: PetscFunctionReturn(0);
347: }
349: /*@
350: DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
351: active part of a matrix.
353: Logically Collective on ds
355: Input Parameters:
356: + ds - the direct solver context
357: - mat - the matrix to modify
359: Level: intermediate
361: .seealso: DSGetMat()
362: @*/
363: PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
364: {
365: PetscScalar *x;
366: PetscInt i,ld,n,l;
370: DSCheckValidMat(ds,mat,2);
372: DSGetDimensions(ds,&n,&l,NULL,NULL);
373: DSGetLeadingDimension(ds,&ld);
374: PetscLogEventBegin(DS_Other,ds,0,0,0);
375: DSGetArray(ds,mat,&x);
376: PetscArrayzero(&x[ld*l],ld*(n-l));
377: for (i=l;i<n;i++) x[i+i*ld] = 1.0;
378: DSRestoreArray(ds,mat,&x);
379: PetscLogEventEnd(DS_Other,ds,0,0,0);
380: PetscFunctionReturn(0);
381: }
383: /*@C
384: DSOrthogonalize - Orthogonalize the columns of a matrix.
386: Logically Collective on ds
388: Input Parameters:
389: + ds - the direct solver context
390: . mat - a matrix
391: - cols - number of columns to orthogonalize (starting from column zero)
393: Output Parameter:
394: . lindcols - (optional) number of linearly independent columns of the matrix
396: Level: developer
398: .seealso: DSPseudoOrthogonalize()
399: @*/
400: PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
401: {
402: PetscInt n,l,ld;
403: PetscBLASInt ld_,rA,cA,info,ltau,lw;
404: PetscScalar *A,*tau,*w,saux,dummy;
407: DSCheckAlloc(ds,1);
409: DSCheckValidMat(ds,mat,2);
412: DSGetDimensions(ds,&n,&l,NULL,NULL);
413: DSGetLeadingDimension(ds,&ld);
414: n = n - l;
416: if (n == 0 || cols == 0) PetscFunctionReturn(0);
418: PetscLogEventBegin(DS_Other,ds,0,0,0);
419: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
420: DSGetArray(ds,mat,&A);
421: PetscBLASIntCast(PetscMin(cols,n),<au);
422: PetscBLASIntCast(ld,&ld_);
423: PetscBLASIntCast(n,&rA);
424: PetscBLASIntCast(cols,&cA);
425: lw = -1;
426: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
427: SlepcCheckLapackInfo("geqrf",info);
428: lw = (PetscBLASInt)PetscRealPart(saux);
429: DSAllocateWork_Private(ds,lw+ltau,0,0);
430: tau = ds->work;
431: w = &tau[ltau];
432: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
433: SlepcCheckLapackInfo("geqrf",info);
434: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
435: SlepcCheckLapackInfo("orgqr",info);
436: if (lindcols) *lindcols = ltau;
438: PetscFPTrapPop();
439: PetscLogEventEnd(DS_Other,ds,0,0,0);
440: DSRestoreArray(ds,mat,&A);
441: PetscObjectStateIncrease((PetscObject)ds);
442: PetscFunctionReturn(0);
443: }
445: /*
446: Compute C <- a*A*B + b*C, where
447: ldC, the leading dimension of C,
448: ldA, the leading dimension of A,
449: rA, cA, rows and columns of A,
450: At, if true use the transpose of A instead,
451: ldB, the leading dimension of B,
452: rB, cB, rows and columns of B,
453: Bt, if true use the transpose of B instead
454: */
455: static PetscErrorCode SlepcMatDenseMult(PetscScalar *C,PetscInt _ldC,PetscScalar b,PetscScalar a,const PetscScalar *A,PetscInt _ldA,PetscInt rA,PetscInt cA,PetscBool At,const PetscScalar *B,PetscInt _ldB,PetscInt rB,PetscInt cB,PetscBool Bt)
456: {
457: PetscInt tmp;
458: PetscBLASInt m, n, k, ldA = _ldA, ldB = _ldB, ldC = _ldC;
459: const char *N = "N", *T = "C", *qA = N, *qB = N;
461: if ((rA == 0) || (cB == 0)) PetscFunctionReturn(0);
466: /* Transpose if needed */
467: if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
468: if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;
470: /* Check size */
473: /* Do stub */
474: if ((rA == 1) && (cA == 1) && (cB == 1)) {
475: if (!At && !Bt) *C = *A * *B;
476: else if (At && !Bt) *C = PetscConj(*A) * *B;
477: else if (!At && Bt) *C = *A * PetscConj(*B);
478: else *C = PetscConj(*A) * PetscConj(*B);
479: m = n = k = 1;
480: } else {
481: m = rA; n = cB; k = cA;
482: PetscStackCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
483: }
485: PetscLogFlops(2.0*m*n*k);
486: PetscFunctionReturn(0);
487: }
489: /*@C
490: DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
491: Gram-Schmidt in an indefinite inner product space defined by a signature.
493: Logically Collective on ds
495: Input Parameters:
496: + ds - the direct solver context
497: . mat - the matrix
498: . cols - number of columns to orthogonalize (starting from column zero)
499: - s - the signature that defines the inner product
501: Output Parameters:
502: + lindcols - (optional) linearly independent columns of the matrix
503: - ns - (optional) the new signature of the vectors
505: Note:
506: After the call the matrix satisfies A'*s*A = ns.
508: Level: developer
510: .seealso: DSOrthogonalize()
511: @*/
512: PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal *s,PetscInt *lindcols,PetscReal *ns)
513: {
514: PetscInt i,j,k,l,n,ld;
515: PetscBLASInt info,one=1,zero=0,rA_,ld_;
516: PetscScalar *A,*A_,*m,*h,nr0;
517: PetscReal nr_o,nr,nr_abs,*ns_,done=1.0;
520: DSCheckAlloc(ds,1);
522: DSCheckValidMat(ds,mat,2);
525: DSGetDimensions(ds,&n,&l,NULL,NULL);
526: DSGetLeadingDimension(ds,&ld);
527: n = n - l;
529: if (n == 0 || cols == 0) PetscFunctionReturn(0);
530: PetscBLASIntCast(n,&rA_);
531: PetscBLASIntCast(ld,&ld_);
532: DSGetArray(ds,mat,&A_);
533: A = &A_[ld*l+l];
534: DSAllocateWork_Private(ds,n+cols,ns?0:cols,0);
535: m = ds->work;
536: h = &m[n];
537: ns_ = ns ? ns : ds->rwork;
538: PetscLogEventBegin(DS_Other,ds,0,0,0);
539: for (i=0; i<cols; i++) {
540: /* m <- diag(s)*A[i] */
541: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
542: /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
543: SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
544: nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
545: for (j=0; j<3 && i>0; j++) {
546: /* h <- A[0:i-1]'*m */
547: SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
548: /* h <- diag(ns)*h */
549: for (k=0; k<i; k++) h[k] *= ns_[k];
550: /* A[i] <- A[i] - A[0:i-1]*h */
551: SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE);
552: /* m <- diag(s)*A[i] */
553: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
554: /* nr_o <- mynorm(A[i]'*m) */
555: SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
556: nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
558: if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
559: nr_o = nr;
560: }
561: ns_[i] = PetscSign(nr);
562: /* A[i] <- A[i]/|nr| */
563: nr_abs = PetscAbs(nr);
564: PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&nr_abs,&done,&rA_,&one,A+i*ld,&ld_,&info));
565: SlepcCheckLapackInfo("lascl",info);
566: }
567: PetscLogEventEnd(DS_Other,ds,0,0,0);
568: DSRestoreArray(ds,mat,&A_);
569: PetscObjectStateIncrease((PetscObject)ds);
570: if (lindcols) *lindcols = cols;
571: PetscFunctionReturn(0);
572: }