Actual source code: dspriv.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: */
10: /*
11: Private DS routines
12: */
14: #include <slepc/private/dsimpl.h> /*I "slepcds.h" I*/
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;
25: PetscObjectTypeCompare((PetscObject)ds,DSPEP,&ispep);
26: if (ispep) {
27: DSPEPGetDegree(ds,&d);
28: }
29: if (ispep && (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;
30: else n = ds->ld;
31: switch (m) {
32: case DS_MAT_T:
33: nelem = 3*ds->ld;
34: break;
35: case DS_MAT_D:
36: nelem = ds->ld;
37: break;
38: case DS_MAT_X:
39: nelem = ds->ld*n;
40: break;
41: case DS_MAT_Y:
42: nelem = ds->ld*n;
43: break;
44: default:
45: nelem = n*n;
46: }
47: if (isreal) {
48: sz = nelem*sizeof(PetscReal);
49: if (ds->rmat[m]) {
50: PetscFree(ds->rmat[m]);
51: } else {
52: PetscLogObjectMemory((PetscObject)ds,sz);
53: }
54: PetscCalloc1(nelem,&ds->rmat[m]);
55: } else {
56: sz = nelem*sizeof(PetscScalar);
57: if (ds->mat[m]) {
58: PetscFree(ds->mat[m]);
59: } else {
60: PetscLogObjectMemory((PetscObject)ds,sz);
61: }
62: PetscCalloc1(nelem,&ds->mat[m]);
63: }
64: return(0);
65: }
67: PetscErrorCode DSAllocateWork_Private(DS ds,PetscInt s,PetscInt r,PetscInt i)
68: {
72: if (s>ds->lwork) {
73: PetscFree(ds->work);
74: PetscMalloc1(s,&ds->work);
75: PetscLogObjectMemory((PetscObject)ds,(s-ds->lwork)*sizeof(PetscScalar));
76: ds->lwork = s;
77: }
78: if (r>ds->lrwork) {
79: PetscFree(ds->rwork);
80: PetscMalloc1(r,&ds->rwork);
81: PetscLogObjectMemory((PetscObject)ds,(r-ds->lrwork)*sizeof(PetscReal));
82: ds->lrwork = r;
83: }
84: if (i>ds->liwork) {
85: PetscFree(ds->iwork);
86: PetscMalloc1(i,&ds->iwork);
87: PetscLogObjectMemory((PetscObject)ds,(i-ds->liwork)*sizeof(PetscBLASInt));
88: ds->liwork = i;
89: }
90: return(0);
91: }
93: /*@C
94: DSViewMat - Prints one of the internal DS matrices.
96: Collective on DS
98: Input Parameters:
99: + ds - the direct solver context
100: . viewer - visualization context
101: - m - matrix to display
103: Note:
104: Works only for ascii viewers. Set the viewer in Matlab format if
105: want to paste into Matlab.
107: Level: developer
109: .seealso: DSView()
110: @*/
111: PetscErrorCode DSViewMat(DS ds,PetscViewer viewer,DSMatType m)
112: {
113: PetscErrorCode ierr;
114: PetscInt i,j,rows,cols;
115: PetscScalar *v;
116: PetscViewerFormat format;
117: #if defined(PETSC_USE_COMPLEX)
118: PetscBool allreal = PETSC_TRUE;
119: #endif
124: DSCheckValidMat(ds,m,3);
125: if (!viewer) {
126: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ds),&viewer);
127: }
130: PetscViewerGetFormat(viewer,&format);
131: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
132: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
133: DSMatGetSize(ds,m,&rows,&cols);
134: #if defined(PETSC_USE_COMPLEX)
135: /* determine if matrix has all real values */
136: v = ds->mat[m];
137: for (i=0;i<rows;i++)
138: for (j=0;j<cols;j++)
139: if (PetscImaginaryPart(v[i+j*ds->ld])) { allreal = PETSC_FALSE; break; }
140: #endif
141: if (format == PETSC_VIEWER_ASCII_MATLAB) {
142: PetscViewerASCIIPrintf(viewer,"%% Size = %D %D\n",rows,cols);
143: PetscViewerASCIIPrintf(viewer,"%s = [\n",DSMatName[m]);
144: } else {
145: PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[m]);
146: }
148: for (i=0;i<rows;i++) {
149: v = ds->mat[m]+i;
150: for (j=0;j<cols;j++) {
151: #if defined(PETSC_USE_COMPLEX)
152: if (allreal) {
153: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)PetscRealPart(*v));
154: } else {
155: PetscViewerASCIIPrintf(viewer,"%18.16e%+18.16ei ",(double)PetscRealPart(*v),(double)PetscImaginaryPart(*v));
156: }
157: #else
158: PetscViewerASCIIPrintf(viewer,"%18.16e ",(double)*v);
159: #endif
160: v += ds->ld;
161: }
162: PetscViewerASCIIPrintf(viewer,"\n");
163: }
165: if (format == PETSC_VIEWER_ASCII_MATLAB) {
166: PetscViewerASCIIPrintf(viewer,"];\n");
167: }
168: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
169: PetscViewerFlush(viewer);
170: return(0);
171: }
173: PetscErrorCode DSSortEigenvalues_Private(DS ds,PetscScalar *wr,PetscScalar *wi,PetscInt *perm,PetscBool isghiep)
174: {
176: PetscScalar re,im,wi0;
177: PetscInt n,i,j,result,tmp1,tmp2=0,d=1;
180: n = ds->t; /* sort only first t pairs if truncated */
181: /* insertion sort */
182: i=ds->l+1;
183: #if !defined(PETSC_USE_COMPLEX)
184: if (wi && wi[perm[i-1]]!=0.0) i++; /* initial value is complex */
185: #else
186: if (isghiep && PetscImaginaryPart(wr[perm[i-1]])!=0.0) i++;
187: #endif
188: for (;i<n;i+=d) {
189: re = wr[perm[i]];
190: if (wi) im = wi[perm[i]];
191: else im = 0.0;
192: tmp1 = perm[i];
193: #if !defined(PETSC_USE_COMPLEX)
194: if (im!=0.0) { d = 2; tmp2 = perm[i+1]; }
195: else d = 1;
196: #else
197: if (isghiep && PetscImaginaryPart(re)!=0.0) { d = 2; tmp2 = perm[i+1]; }
198: else d = 1;
199: #endif
200: j = i-1;
201: if (wi) wi0 = wi[perm[j]];
202: else wi0 = 0.0;
203: SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
204: while (result<0 && j>=ds->l) {
205: perm[j+d] = perm[j];
206: j--;
207: #if !defined(PETSC_USE_COMPLEX)
208: if (wi && wi[perm[j+1]]!=0)
209: #else
210: if (isghiep && PetscImaginaryPart(wr[perm[j+1]])!=0)
211: #endif
212: { perm[j+d] = perm[j]; j--; }
214: if (j>=ds->l) {
215: if (wi) wi0 = wi[perm[j]];
216: else wi0 = 0.0;
217: SlepcSCCompare(ds->sc,re,im,wr[perm[j]],wi0,&result);
218: }
219: }
220: perm[j+1] = tmp1;
221: if (d==2) perm[j+2] = tmp2;
222: }
223: return(0);
224: }
226: PetscErrorCode DSSortEigenvaluesReal_Private(DS ds,PetscReal *eig,PetscInt *perm)
227: {
229: PetscScalar re;
230: PetscInt i,j,result,tmp,l,n;
233: n = ds->t; /* sort only first t pairs if truncated */
234: l = ds->l;
235: /* insertion sort */
236: for (i=l+1;i<n;i++) {
237: re = eig[perm[i]];
238: j = i-1;
239: SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
240: while (result<0 && j>=l) {
241: tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
242: if (j>=l) {
243: SlepcSCCompare(ds->sc,re,0.0,eig[perm[j]],0.0,&result);
244: }
245: }
246: }
247: return(0);
248: }
250: /*
251: DSCopyMatrix_Private - Copies the trailing block of a matrix (from
252: rows/columns l to n).
253: */
254: PetscErrorCode DSCopyMatrix_Private(DS ds,DSMatType dst,DSMatType src)
255: {
257: PetscInt j,m,off,ld;
258: PetscScalar *S,*D;
261: ld = ds->ld;
262: m = ds->n-ds->l;
263: off = ds->l+ds->l*ld;
264: S = ds->mat[src];
265: D = ds->mat[dst];
266: for (j=0;j<m;j++) {
267: PetscMemcpy(D+off+j*ld,S+off+j*ld,m*sizeof(PetscScalar));
268: }
269: return(0);
270: }
272: PetscErrorCode DSPermuteColumns_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat,PetscInt *perm)
273: {
274: PetscInt i,j,k,p,ld;
275: PetscScalar *Q,rtmp;
278: ld = ds->ld;
279: Q = ds->mat[mat];
280: for (i=l;i<n;i++) {
281: p = perm[i];
282: if (p != i) {
283: j = i + 1;
284: while (perm[j] != i) j++;
285: perm[j] = p; perm[i] = i;
286: /* swap columns i and j */
287: for (k=0;k<n;k++) {
288: rtmp = Q[k+p*ld]; Q[k+p*ld] = Q[k+i*ld]; Q[k+i*ld] = rtmp;
289: }
290: }
291: }
292: return(0);
293: }
295: PetscErrorCode DSPermuteRows_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat,PetscInt *perm)
296: {
297: PetscInt i,j,m=ds->m,k,p,ld;
298: PetscScalar *Q,rtmp;
301: if (!m) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"m was not set");
302: ld = ds->ld;
303: Q = ds->mat[mat];
304: for (i=l;i<n;i++) {
305: p = perm[i];
306: if (p != i) {
307: j = i + 1;
308: while (perm[j] != i) j++;
309: perm[j] = p; perm[i] = i;
310: /* swap rows i and j */
311: for (k=0;k<m;k++) {
312: rtmp = Q[p+k*ld]; Q[p+k*ld] = Q[i+k*ld]; Q[i+k*ld] = rtmp;
313: }
314: }
315: }
316: return(0);
317: }
319: PetscErrorCode DSPermuteBoth_Private(DS ds,PetscInt l,PetscInt n,DSMatType mat1,DSMatType mat2,PetscInt *perm)
320: {
321: PetscInt i,j,m=ds->m,k,p,ld;
322: PetscScalar *U,*VT,rtmp;
325: if (!m) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"m was not set");
326: ld = ds->ld;
327: U = ds->mat[mat1];
328: VT = ds->mat[mat2];
329: for (i=l;i<n;i++) {
330: p = perm[i];
331: if (p != i) {
332: j = i + 1;
333: while (perm[j] != i) j++;
334: perm[j] = p; perm[i] = i;
335: /* swap columns i and j of U */
336: for (k=0;k<n;k++) {
337: rtmp = U[k+p*ld]; U[k+p*ld] = U[k+i*ld]; U[k+i*ld] = rtmp;
338: }
339: /* swap rows i and j of VT */
340: for (k=0;k<m;k++) {
341: rtmp = VT[p+k*ld]; VT[p+k*ld] = VT[i+k*ld]; VT[i+k*ld] = rtmp;
342: }
343: }
344: }
345: return(0);
346: }
348: /*@
349: DSSetIdentity - Copy the identity (a diagonal matrix with ones) on the
350: active part of a matrix.
352: Logically Collective on DS
354: Input Parameters:
355: + ds - the direct solver context
356: - mat - the matrix to modify
358: Level: intermediate
359: @*/
360: PetscErrorCode DSSetIdentity(DS ds,DSMatType mat)
361: {
363: PetscScalar *x;
364: PetscInt i,ld,n,l;
369: DSCheckValidMat(ds,mat,2);
371: DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
372: DSGetLeadingDimension(ds,&ld);
373: PetscLogEventBegin(DS_Other,ds,0,0,0);
374: DSGetArray(ds,mat,&x);
375: PetscMemzero(&x[ld*l],ld*(n-l)*sizeof(PetscScalar));
376: for (i=l;i<n;i++) x[i+i*ld] = 1.0;
377: DSRestoreArray(ds,mat,&x);
378: PetscLogEventEnd(DS_Other,ds,0,0,0);
379: return(0);
380: }
382: /*@C
383: DSOrthogonalize - Orthogonalize the columns of a matrix.
385: Logically Collective on DS
387: Input Parameters:
388: + ds - the direct solver context
389: . mat - a matrix
390: - cols - number of columns to orthogonalize (starting from column zero)
392: Output Parameter:
393: . lindcols - (optional) number of linearly independent columns of the matrix
395: Level: developer
397: .seealso: DSPseudoOrthogonalize()
398: @*/
399: PetscErrorCode DSOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscInt *lindcols)
400: {
401: #if defined(PETSC_MISSING_LAPACK_GEQRF) || defined(PETSC_MISSING_LAPACK_ORGQR)
403: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GEQRF/ORGQR - Lapack routine is unavailable");
404: #else
406: PetscInt n,l,ld;
407: PetscBLASInt ld_,rA,cA,info,ltau,lw;
408: PetscScalar *A,*tau,*w,saux,dummy;
412: DSCheckAlloc(ds,1);
414: DSCheckValidMat(ds,mat,2);
417: DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
418: DSGetLeadingDimension(ds,&ld);
419: n = n - l;
420: if (cols > n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
421: if (n == 0 || cols == 0) return(0);
423: PetscLogEventBegin(DS_Other,ds,0,0,0);
424: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
425: DSGetArray(ds,mat,&A);
426: PetscBLASIntCast(PetscMin(cols,n),<au);
427: PetscBLASIntCast(ld,&ld_);
428: PetscBLASIntCast(n,&rA);
429: PetscBLASIntCast(cols,&cA);
430: lw = -1;
431: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,A,&ld_,&dummy,&saux,&lw,&info));
432: SlepcCheckLapackInfo("geqrf",info);
433: lw = (PetscBLASInt)PetscRealPart(saux);
434: DSAllocateWork_Private(ds,lw+ltau,0,0);
435: tau = ds->work;
436: w = &tau[ltau];
437: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&rA,&cA,&A[ld*l+l],&ld_,tau,w,&lw,&info));
438: SlepcCheckLapackInfo("geqrf",info);
439: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&rA,<au,<au,&A[ld*l+l],&ld_,tau,w,&lw,&info));
440: SlepcCheckLapackInfo("orgqr",info);
441: if (lindcols) *lindcols = ltau;
443: PetscFPTrapPop();
444: PetscLogEventEnd(DS_Other,ds,0,0,0);
445: DSRestoreArray(ds,mat,&A);
446: PetscObjectStateIncrease((PetscObject)ds);
447: return(0);
448: #endif
449: }
451: /*
452: Compute C <- a*A*B + b*C, where
453: ldC, the leading dimension of C,
454: ldA, the leading dimension of A,
455: rA, cA, rows and columns of A,
456: At, if true use the transpose of A instead,
457: ldB, the leading dimension of B,
458: rB, cB, rows and columns of B,
459: Bt, if true use the transpose of B instead
460: */
461: 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)
462: {
464: PetscInt tmp;
465: PetscBLASInt m, n, k, ldA = _ldA, ldB = _ldB, ldC = _ldC;
466: const char *N = "N", *T = "C", *qA = N, *qB = N;
469: if ((rA == 0) || (cB == 0)) return(0);
474: /* Transpose if needed */
475: if (At) tmp = rA, rA = cA, cA = tmp, qA = T;
476: if (Bt) tmp = rB, rB = cB, cB = tmp, qB = T;
478: /* Check size */
479: if (cA != rB) SETERRQ(PETSC_COMM_SELF,1,"Matrix dimensions do not match");
481: /* Do stub */
482: if ((rA == 1) && (cA == 1) && (cB == 1)) {
483: if (!At && !Bt) *C = *A * *B;
484: else if (At && !Bt) *C = PetscConj(*A) * *B;
485: else if (!At && Bt) *C = *A * PetscConj(*B);
486: else *C = PetscConj(*A) * PetscConj(*B);
487: m = n = k = 1;
488: } else {
489: m = rA; n = cB; k = cA;
490: PetscStackCallBLAS("BLASgemm",BLASgemm_(qA,qB,&m,&n,&k,&a,(PetscScalar*)A,&ldA,(PetscScalar*)B,&ldB,&b,C,&ldC));
491: }
493: PetscLogFlops(2.0*m*n*k);
494: return(0);
495: }
497: /*@C
498: DSPseudoOrthogonalize - Orthogonalize the columns of a matrix with Modified
499: Gram-Schmidt in an indefinite inner product space defined by a signature.
501: Logically Collective on DS
503: Input Parameters:
504: + ds - the direct solver context
505: . mat - the matrix
506: . cols - number of columns to orthogonalize (starting from column zero)
507: - s - the signature that defines the inner product
509: Output Parameters:
510: + lindcols - (optional) linearly independent columns of the matrix
511: - ns - (optional) the new signature of the vectors
513: Note:
514: After the call the matrix satisfies A'*s*A = ns.
516: Level: developer
518: .seealso: DSOrthogonalize()
519: @*/
520: PetscErrorCode DSPseudoOrthogonalize(DS ds,DSMatType mat,PetscInt cols,PetscReal *s,PetscInt *lindcols,PetscReal *ns)
521: {
523: PetscInt i,j,k,l,n,ld;
524: PetscBLASInt one=1,rA_;
525: PetscScalar alpha,*A,*A_,*m,*h,nr0;
526: PetscReal nr_o,nr,*ns_;
530: DSCheckAlloc(ds,1);
532: DSCheckValidMat(ds,mat,2);
535: DSGetDimensions(ds,&n,NULL,&l,NULL,NULL);
536: DSGetLeadingDimension(ds,&ld);
537: n = n - l;
538: if (cols > n) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONG,"Invalid number of columns");
539: if (n == 0 || cols == 0) return(0);
540: PetscBLASIntCast(n,&rA_);
541: DSGetArray(ds,mat,&A_);
542: A = &A_[ld*l+l];
543: DSAllocateWork_Private(ds,n+cols,ns?0:cols,0);
544: m = ds->work;
545: h = &m[n];
546: ns_ = ns ? ns : ds->rwork;
547: PetscLogEventBegin(DS_Other,ds,0,0,0);
548: for (i=0; i<cols; i++) {
549: /* m <- diag(s)*A[i] */
550: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
551: /* nr_o <- mynorm(A[i]'*m), mynorm(x) = sign(x)*sqrt(|x|) */
552: SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
553: nr = nr_o = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
554: for (j=0; j<3 && i>0; j++) {
555: /* h <- A[0:i-1]'*m */
556: SlepcMatDenseMult(h,i,0.0,1.0,A,ld,n,i,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
557: /* h <- diag(ns)*h */
558: for (k=0; k<i; k++) h[k] *= ns_[k];
559: /* A[i] <- A[i] - A[0:i-1]*h */
560: SlepcMatDenseMult(&A[ld*i],ld,1.0,-1.0,A,ld,n,i,PETSC_FALSE,h,i,i,1,PETSC_FALSE);
561: /* m <- diag(s)*A[i] */
562: for (k=0; k<n; k++) m[k] = s[k]*A[k+i*ld];
563: /* nr_o <- mynorm(A[i]'*m) */
564: SlepcMatDenseMult(&nr0,1,0.0,1.0,&A[ld*i],ld,n,1,PETSC_TRUE,m,n,n,1,PETSC_FALSE);
565: nr = PetscSign(PetscRealPart(nr0))*PetscSqrtReal(PetscAbsScalar(nr0));
566: if (PetscAbs(nr) < PETSC_MACHINE_EPSILON) SETERRQ(PETSC_COMM_SELF,1,"Linear dependency detected");
567: if (PetscAbs(nr) > 0.7*PetscAbs(nr_o)) break;
568: nr_o = nr;
569: }
570: ns_[i] = PetscSign(nr);
571: /* A[i] <- A[i]/|nr| */
572: alpha = 1.0/PetscAbs(nr);
573: PetscStackCallBLAS("BLASscal",BLASscal_(&rA_,&alpha,&A[i*ld],&one));
574: }
575: PetscLogEventEnd(DS_Other,ds,0,0,0);
576: DSRestoreArray(ds,mat,&A_);
577: PetscObjectStateIncrease((PetscObject)ds);
578: if (lindcols) *lindcols = cols;
579: return(0);
580: }