Actual source code: dshsvd.c
slepc-3.22.1 2024-10-28
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: typedef struct {
15: PetscInt m; /* number of columns */
16: PetscInt t; /* number of rows of V after truncating */
17: PetscBool reorth; /* reorthogonalize left vectors */
18: } DS_HSVD;
20: static PetscErrorCode DSAllocate_HSVD(DS ds,PetscInt ld)
21: {
22: PetscFunctionBegin;
23: if (!ds->compact) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24: PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
25: PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
26: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
27: PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
28: PetscCall(PetscFree(ds->perm));
29: PetscCall(PetscMalloc1(ld,&ds->perm));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: /* 0 l k m-1
34: -----------------------------------------
35: |* . . |
36: | * . . |
37: | * . . |
38: | * . . |
39: | o o |
40: | o o |
41: | o o |
42: | o o |
43: | o o |
44: | o o |
45: | o x |
46: | x x |
47: | x x |
48: | x x |
49: | x x |
50: | x x |
51: | x x |
52: | x x |
53: | x x|
54: n-1 | x|
55: -----------------------------------------
56: */
58: static PetscErrorCode DSView_HSVD(DS ds,PetscViewer viewer)
59: {
60: DS_HSVD *ctx = (DS_HSVD*)ds->data;
61: PetscViewerFormat format;
62: PetscInt i,j,r,c,m=ctx->m,rows,cols;
63: PetscReal *T,*S,value;
64: const char *methodname[] = {
65: "Cross product A'*Omega*A"
66: };
67: const int nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);
69: PetscFunctionBegin;
70: PetscCall(PetscViewerGetFormat(viewer,&format));
71: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
72: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
73: if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
74: if (ctx->reorth) PetscCall(PetscViewerASCIIPrintf(viewer,"reorthogonalizing left vectors\n"));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
77: PetscCheck(m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
78: if (ds->compact) {
79: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
80: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
81: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
82: rows = ds->n;
83: cols = ds->extrarow? m+1: m;
84: if (format == PETSC_VIEWER_ASCII_MATLAB) {
85: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rows,cols));
86: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
87: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
88: for (i=0;i<PetscMin(ds->n,m);i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
89: for (i=0;i<cols-1;i++) {
90: c = PetscMax(i+2,ds->k+1);
91: r = i+1;
92: value = i<ds->l? 0.0: T[i+ds->ld];
93: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",r,c,(double)value));
94: }
95: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
96: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",ds->n,ds->n));
97: PetscCall(PetscViewerASCIIPrintf(viewer,"omega = zeros(%" PetscInt_FMT ",3);\n",3*ds->n));
98: PetscCall(PetscViewerASCIIPrintf(viewer,"omega = [\n"));
99: for (i=0;i<ds->n;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
100: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(omega);\n",DSMatName[DS_MAT_D]));
101: } else {
102: PetscCall(PetscViewerASCIIPrintf(viewer,"T\n"));
103: for (i=0;i<rows;i++) {
104: for (j=0;j<cols;j++) {
105: if (i==j) value = T[i];
106: else if (i<ds->l) value = 0.0;
107: else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+ds->ld];
108: else if (i+1==j && i>=ds->k) value = T[i+ds->ld];
109: else value = 0.0;
110: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
111: }
112: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
113: }
114: PetscCall(PetscViewerASCIIPrintf(viewer,"omega\n"));
115: for (i=0;i<ds->n;i++) {
116: for (j=0;j<ds->n;j++) {
117: if (i==j) value = S[i];
118: else value = 0.0;
119: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
120: }
121: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
122: }
123: }
124: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
125: PetscCall(PetscViewerFlush(viewer));
126: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
127: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
128: } else {
129: PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
130: PetscCall(DSViewMat(ds,viewer,DS_MAT_D));
131: }
132: if (ds->state>DS_STATE_INTERMEDIATE) {
133: PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
134: PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
135: }
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode DSVectors_HSVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
140: {
141: PetscFunctionBegin;
142: switch (mat) {
143: case DS_MAT_U:
144: case DS_MAT_V:
145: if (rnorm) *rnorm = 0.0;
146: break;
147: default:
148: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
149: }
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode DSSort_HSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
154: {
155: DS_HSVD *ctx = (DS_HSVD*)ds->data;
156: PetscInt n,l,i,*perm,ld=ds->ld;
157: PetscScalar *A;
158: PetscReal *d,*s;
160: PetscFunctionBegin;
161: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
162: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
163: l = ds->l;
164: n = PetscMin(ds->n,ctx->m);
165: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
166: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&s));
167: PetscCall(DSAllocateWork_Private(ds,0,ds->ld,0));
168: perm = ds->perm;
169: if (!rr) PetscCall(DSSortEigenvaluesReal_Private(ds,d,perm));
170: else PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
171: PetscCall(PetscArraycpy(ds->rwork,s,n));
172: for (i=l;i<n;i++) s[i] = ds->rwork[perm[i]];
173: for (i=l;i<n;i++) wr[i] = d[perm[i]];
174: PetscCall(DSPermuteBoth_Private(ds,l,n,ds->n,ctx->m,DS_MAT_U,DS_MAT_V,perm));
175: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
176: if (!ds->compact) {
177: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
178: for (i=l;i<n;i++) A[i+i*ld] = wr[i];
179: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
180: }
181: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
182: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&s));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: static PetscErrorCode DSUpdateExtraRow_HSVD(DS ds)
187: {
188: DS_HSVD *ctx = (DS_HSVD*)ds->data;
189: PetscInt i;
190: PetscBLASInt n=0,m=0,ld,l;
191: const PetscScalar *U;
192: PetscReal *T,*e,*Omega,beta;
194: PetscFunctionBegin;
195: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
196: PetscCall(PetscBLASIntCast(ds->n,&n));
197: PetscCall(PetscBLASIntCast(ctx->m,&m));
198: PetscCall(PetscBLASIntCast(ds->ld,&ld));
199: PetscCall(PetscBLASIntCast(ds->l,&l));
200: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
201: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
202: PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
203: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
204: e = T+ld;
205: beta = PetscAbs(e[m-1]); /* in compact, we assume all entries are zero except the last one */
206: for (i=0;i<n;i++) e[i] = PetscRealPart(beta*U[n-1+i*ld]*Omega[i]);
207: ds->k = m;
208: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
209: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
210: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode DSTruncate_HSVD(DS ds,PetscInt n,PetscBool trim)
215: {
216: PetscInt i,ld=ds->ld,l=ds->l;
217: PetscScalar *A;
218: DS_HSVD *ctx = (DS_HSVD*)ds->data;
220: PetscFunctionBegin;
221: if (!ds->compact && ds->extrarow) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
222: if (trim) {
223: if (!ds->compact && ds->extrarow) { /* clean extra column */
224: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
225: }
226: ds->l = 0;
227: ds->k = 0;
228: ds->n = n;
229: ctx->m = n;
230: ds->t = ds->n; /* truncated length equal to the new dimension */
231: ctx->t = ctx->m; /* must also keep the previous dimension of V */
232: } else {
233: if (!ds->compact && ds->extrarow && ds->k==ds->n) {
234: /* copy entries of extra column to the new position, then clean last row */
235: for (i=l;i<n;i++) A[i+n*ld] = A[i+ctx->m*ld];
236: for (i=l;i<ds->n;i++) A[i+ctx->m*ld] = 0.0;
237: }
238: ds->k = (ds->extrarow)? n: 0;
239: ds->t = ds->n; /* truncated length equal to previous dimension */
240: ctx->t = ctx->m; /* must also keep the previous dimension of V */
241: ds->n = n;
242: ctx->m = n;
243: }
244: if (!ds->compact && ds->extrarow) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode DSSolve_HSVD_CROSS(DS ds,PetscScalar *wr,PetscScalar *wi)
249: {
250: DS_HSVD *ctx = (DS_HSVD*)ds->data;
251: PetscInt i,j,k=ds->k,rwu=0,iwu=0,swu=0,nv;
252: PetscBLASInt n1,n2,info,l=0,n=0,m=0,ld,off,one=1,*perm,*cmplx,incx=1,lwork;
253: PetscScalar *A,*U,*V,scal,*R,sone=1.0,szero=0.0;
254: PetscReal *d,*e,*dd,*ee,*Omega;
256: PetscFunctionBegin;
257: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
258: PetscCall(PetscBLASIntCast(ds->n,&n));
259: PetscCall(PetscBLASIntCast(ctx->m,&m));
260: PetscCheck(!ds->compact || n==m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-square matrices in compact storage");
261: PetscCheck(ds->compact || n>=m,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for the case of more columns than rows");
262: PetscCall(PetscBLASIntCast(ds->l,&l));
263: PetscCall(PetscBLASIntCast(ds->ld,&ld));
264: PetscCall(PetscBLASIntCast(PetscMax(0,ds->k-ds->l+1),&n2));
265: n1 = n-l; /* n1 = size of leading block, excl. locked + size of trailing block */
266: off = l+l*ld;
267: if (!ds->compact) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
268: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_U],&U));
269: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_V],&V));
270: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&d));
271: e = d+ld;
272: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&Omega));
273: PetscCall(PetscArrayzero(U,ld*ld));
274: for (i=0;i<l;i++) U[i+i*ld] = 1.0;
275: PetscCall(PetscArrayzero(V,ld*ld));
276: for (i=0;i<n;i++) V[i+i*ld] = 1.0;
277: for (i=0;i<l;i++) wr[i] = d[i];
278: if (wi) for (i=0;i<l;i++) wi[i] = 0.0;
280: if (ds->compact) {
281: /* Form the arrow tridiagonal cross product T=A'*Omega*A, where A is the arrow
282: bidiagonal matrix formed by d, e. T is stored in dd, ee */
283: PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,4*ld,2*ld));
284: R = ds->work+swu;
285: swu += n*ld;
286: perm = ds->iwork+iwu;
287: iwu += n;
288: cmplx = ds->iwork+iwu;
289: dd = ds->rwork+rwu;
290: rwu += ld;
291: ee = ds->rwork+rwu;
292: rwu += ld;
293: for (i=0;i<l;i++) {dd[i] = d[i]*d[i]*Omega[i]; ee[i] = 0.0;}
294: for (i=l;i<=ds->k;i++) {
295: dd[i] = Omega[i]*d[i]*d[i];
296: ee[i] = Omega[i]*d[i]*e[i];
297: }
298: for (i=l;i<k;i++) dd[k] += Omega[i]*e[i]*e[i];
299: for (i=k+1;i<n;i++) {
300: dd[i] = Omega[i]*d[i]*d[i]+Omega[i-1]*e[i-1]*e[i-1];
301: ee[i] = Omega[i]*d[i]*e[i];
302: }
304: /* Reduce T to tridiagonal form */
305: PetscCall(DSArrowTridiag(n2,dd+l,ee+l,V+off,ld));
307: /* Solve the tridiagonal eigenproblem corresponding to T */
308: PetscCallBLAS("LAPACKsteqr",LAPACKsteqr_("V",&n1,dd+l,ee+l,V+off,&ld,ds->rwork+rwu,&info));
309: SlepcCheckLapackInfo("steqr",info);
310: for (i=l;i<n;i++) wr[i] = PetscSqrtScalar(PetscAbs(dd[i]));
312: /* Build left singular vectors: U=A*V*Sigma^-1 */
313: PetscCall(PetscArrayzero(U+l*ld,n1*ld));
314: for (i=l;i<n-1;i++) {
315: scal = d[i];
316: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+i,&ld,U+l*ld+i,&ld));
317: j = (i<k)?k:i+1;
318: scal = e[i];
319: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+j,&ld,U+l*ld+i,&ld));
320: }
321: scal = d[n-1];
322: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+off+(n1-1),&ld,U+off+(n1-1),&ld));
323: /* Multiply by Sigma^-1 */
324: for (i=l;i<n;i++) {scal = 1.0/wr[i]; PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,U+i*ld+l,&one));}
326: } else { /* non-compact */
328: PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,PetscDefined(USE_COMPLEX)?4*ld:ld,2*ld));
329: R = ds->work+swu;
330: swu += n*ld;
331: perm = ds->iwork+iwu;
332: iwu += n;
333: cmplx = ds->iwork+iwu;
334: dd = ds->rwork+rwu;
335: for (j=l;j<m;j++) {
336: for (i=0;i<n;i++) ds->work[i] = Omega[i]*A[i+j*ld];
337: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&m,&sone,A,&ld,ds->work,&incx,&szero,V+j*ld,&incx));
338: }
340: /* compute eigenvalues */
341: lwork = (n+6)*ld;
342: #if defined(PETSC_USE_COMPLEX)
343: rwu += ld;
344: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,ds->rwork+rwu,&info));
345: #else
346: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,&info));
347: #endif
348: SlepcCheckLapackInfo("syev",info);
349: for (i=l;i<PetscMin(n,m);i++) d[i] = PetscSqrtReal(PetscAbsReal(dd[i]));
351: /* Build left singular vectors: U=A*V*Sigma^-1 */
352: for (j=l;j<PetscMin(n,m);j++) {
353: scal = 1.0/d[j];
354: PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&m,&scal,A,&ld,V+j*ld,&incx,&szero,U+j*ld,&incx));
355: }
356: }
358: if (ctx->reorth) { /* Reinforce orthogonality */
359: nv = n1;
360: for (i=0;i<n;i++) cmplx[i] = 0;
361: PetscCall(DSPseudoOrthog_HR(&nv,U+off,ld,Omega+l,R,ld,perm,cmplx,NULL,ds->work+swu));
362: } else { /* Update Omega */
363: for (i=l;i<PetscMin(n,m);i++) Omega[i] = PetscSign(dd[i]);
364: }
366: /* Update projected problem */
367: if (ds->compact) {
368: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
369: PetscCall(PetscArrayzero(e,n-1));
370: } else {
371: for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
372: for (i=l;i<n;i++) A[i+i*ld] = d[i];
373: }
374: for (i=l;i<PetscMin(n,m);i++) wr[i] = d[i];
375: if (wi) for (i=l;i<PetscMin(n,m);i++) wi[i] = 0.0;
377: if (ctx->reorth) { /* Update vectors V with R */
378: scal = -1.0;
379: for (i=0;i<nv;i++) {
380: if (PetscRealPart(R[i+i*ld]) < 0.0) PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,V+(i+l)*ld+l,&one));
381: }
382: }
384: if (!ds->compact) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
385: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
386: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
387: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
388: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: #if !defined(PETSC_HAVE_MPIUNI)
393: static PetscErrorCode DSSynchronize_HSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
394: {
395: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
396: PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
397: PetscScalar *A,*U,*V;
398: PetscReal *T,*D;
400: PetscFunctionBegin;
401: if (ds->compact) kr = 3*ld;
402: else k = (ds->n-l)*ld;
403: kr += ld;
404: if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
405: if (eigr) k += ds->n-l;
406: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
407: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
408: PetscCall(PetscMPIIntCast(ds->n-l,&n));
409: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
410: PetscCall(PetscMPIIntCast(3*ld,&ld3));
411: PetscCall(PetscMPIIntCast(ld,&ld_));
412: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
413: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
414: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
415: if (ds->state>DS_STATE_RAW) {
416: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
417: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
418: }
419: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
420: if (!rank) {
421: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
422: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
423: PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
424: if (ds->state>DS_STATE_RAW) {
425: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
426: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
427: }
428: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
429: }
430: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
431: if (rank) {
432: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
433: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
434: PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
435: if (ds->state>DS_STATE_RAW) {
436: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
437: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
438: }
439: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
440: }
441: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
442: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
443: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
444: if (ds->state>DS_STATE_RAW) {
445: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
446: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
447: }
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
450: #endif
452: static PetscErrorCode DSMatGetSize_HSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
453: {
454: DS_HSVD *ctx = (DS_HSVD*)ds->data;
456: PetscFunctionBegin;
457: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
458: switch (t) {
459: case DS_MAT_A:
460: *rows = ds->n;
461: *cols = ds->extrarow? ctx->m+1: ctx->m;
462: break;
463: case DS_MAT_T:
464: *rows = ds->n;
465: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
466: break;
467: case DS_MAT_D:
468: *rows = ds->n;
469: *cols = 1;
470: break;
471: case DS_MAT_U:
472: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
473: *cols = ds->n;
474: break;
475: case DS_MAT_V:
476: *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
477: *cols = ctx->m;
478: break;
479: default:
480: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
481: }
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: static PetscErrorCode DSHSVDSetDimensions_HSVD(DS ds,PetscInt m)
486: {
487: DS_HSVD *ctx = (DS_HSVD*)ds->data;
489: PetscFunctionBegin;
490: DSCheckAlloc(ds,1);
491: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
492: ctx->m = ds->ld;
493: } else {
494: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
495: ctx->m = m;
496: }
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: DSHSVDSetDimensions - Sets the number of columns for a DSHSVD.
503: Logically Collective
505: Input Parameters:
506: + ds - the direct solver context
507: - m - the number of columns
509: Notes:
510: This call is complementary to DSSetDimensions(), to provide a dimension
511: that is specific to this DS type.
513: Level: intermediate
515: .seealso: DSHSVDGetDimensions(), DSSetDimensions()
516: @*/
517: PetscErrorCode DSHSVDSetDimensions(DS ds,PetscInt m)
518: {
519: PetscFunctionBegin;
522: PetscTryMethod(ds,"DSHSVDSetDimensions_C",(DS,PetscInt),(ds,m));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: static PetscErrorCode DSHSVDGetDimensions_HSVD(DS ds,PetscInt *m)
527: {
528: DS_HSVD *ctx = (DS_HSVD*)ds->data;
530: PetscFunctionBegin;
531: *m = ctx->m;
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: DSHSVDGetDimensions - Returns the number of columns for a DSHSVD.
538: Not Collective
540: Input Parameter:
541: . ds - the direct solver context
543: Output Parameters:
544: . m - the number of columns
546: Level: intermediate
548: .seealso: DSHSVDSetDimensions()
549: @*/
550: PetscErrorCode DSHSVDGetDimensions(DS ds,PetscInt *m)
551: {
552: PetscFunctionBegin;
554: PetscAssertPointer(m,2);
555: PetscUseMethod(ds,"DSHSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: static PetscErrorCode DSHSVDSetReorthogonalize_HSVD(DS ds,PetscBool reorth)
560: {
561: DS_HSVD *ctx = (DS_HSVD*)ds->data;
563: PetscFunctionBegin;
564: ctx->reorth = reorth;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: /*@
569: DSHSVDSetReorthogonalize - Sets the reorthogonalization of the left vectors in a DSHSVD.
571: Logically Collective
573: Input Parameters:
574: + ds - the direct solver context
575: - reorth - the reorthogonalization flag
577: Options Database Key:
578: . -ds_hsvd_reorthog <bool> - sets the reorthogonalization flag
580: Note:
581: The computed left vectors (U) should be orthogonal with respect to the signature (D).
582: But it may be necessary to enforce this with a final reorthogonalization step (omitted
583: by default).
585: Level: intermediate
587: .seealso: DSHSVDGetReorthogonalize()
588: @*/
589: PetscErrorCode DSHSVDSetReorthogonalize(DS ds,PetscBool reorth)
590: {
591: PetscFunctionBegin;
594: PetscTryMethod(ds,"DSHSVDSetReorthogonalize_C",(DS,PetscBool),(ds,reorth));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode DSHSVDGetReorthogonalize_HSVD(DS ds,PetscBool *reorth)
599: {
600: DS_HSVD *ctx = (DS_HSVD*)ds->data;
602: PetscFunctionBegin;
603: *reorth = ctx->reorth;
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: /*@
608: DSHSVDGetReorthogonalize - Returns the reorthogonalization flag of a DSHSVD.
610: Not Collective
612: Input Parameter:
613: . ds - the direct solver context
615: Output Parameters:
616: . reorth - the reorthogonalization flag
618: Level: intermediate
620: .seealso: DSHSVDSetReorthogonalize()
621: @*/
622: PetscErrorCode DSHSVDGetReorthogonalize(DS ds,PetscBool *reorth)
623: {
624: PetscFunctionBegin;
626: PetscAssertPointer(reorth,2);
627: PetscUseMethod(ds,"DSHSVDGetReorthogonalize_C",(DS,PetscBool*),(ds,reorth));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: static PetscErrorCode DSSetFromOptions_HSVD(DS ds,PetscOptionItems *PetscOptionsObject)
632: {
633: PetscBool flg,reorth;
635: PetscFunctionBegin;
636: PetscOptionsHeadBegin(PetscOptionsObject,"DS HSVD Options");
638: PetscCall(PetscOptionsBool("-ds_hsvd_reorthog","Reorthogonalize U vectors","DSHSVDSetReorthogonalize",PETSC_FALSE,&reorth,&flg));
639: if (flg) PetscCall(DSHSVDSetReorthogonalize(ds,reorth));
641: PetscOptionsHeadEnd();
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: static PetscErrorCode DSDestroy_HSVD(DS ds)
646: {
647: PetscFunctionBegin;
648: PetscCall(PetscFree(ds->data));
649: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",NULL));
650: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",NULL));
651: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",NULL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",NULL));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: static PetscErrorCode DSSetCompact_HSVD(DS ds,PetscBool comp)
657: {
658: PetscFunctionBegin;
659: if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: /*MC
664: DSHSVD - Dense Hyperbolic Singular Value Decomposition.
666: Level: beginner
668: Notes:
669: The problem is expressed as A = U*Sigma*V', where A is rectangular in
670: general, with n rows and m columns. U is orthogonal with respect to a
671: signature matrix, stored in D. V is orthogonal. Sigma is a diagonal
672: matrix whose diagonal elements are the arguments of DSSolve(). After
673: solve, A is overwritten with Sigma, D is overwritten with the new signature.
675: The matrices of left and right singular vectors, U and V, have size n and m,
676: respectively. The number of columns m must be specified via DSHSVDSetDimensions().
678: If the DS object is in the intermediate state, A is assumed to be in upper
679: bidiagonal form (possibly with an arrow) and is stored in compact format
680: on matrix T. The compact storage is implemented for the square case
681: only, m=n. The extra row should be interpreted in this case as an extra column.
683: Used DS matrices:
684: + DS_MAT_A - problem matrix (used only if compact=false)
685: . DS_MAT_T - upper bidiagonal matrix
686: . DS_MAT_D - diagonal matrix (signature)
687: . DS_MAT_U - left singular vectors
688: - DS_MAT_V - right singular vectors
690: Implemented methods:
691: . 0 - Cross product A'*Omega*A
693: .seealso: DSCreate(), DSSetType(), DSType, DSHSVDSetDimensions()
694: M*/
695: SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS ds)
696: {
697: DS_HSVD *ctx;
699: PetscFunctionBegin;
700: PetscCall(PetscNew(&ctx));
701: ds->data = (void*)ctx;
703: ds->ops->allocate = DSAllocate_HSVD;
704: ds->ops->setfromoptions = DSSetFromOptions_HSVD;
705: ds->ops->view = DSView_HSVD;
706: ds->ops->vectors = DSVectors_HSVD;
707: ds->ops->solve[0] = DSSolve_HSVD_CROSS;
708: ds->ops->sort = DSSort_HSVD;
709: ds->ops->truncate = DSTruncate_HSVD;
710: ds->ops->update = DSUpdateExtraRow_HSVD;
711: ds->ops->destroy = DSDestroy_HSVD;
712: ds->ops->matgetsize = DSMatGetSize_HSVD;
713: #if !defined(PETSC_HAVE_MPIUNI)
714: ds->ops->synchronize = DSSynchronize_HSVD;
715: #endif
716: ds->ops->setcompact = DSSetCompact_HSVD;
717: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",DSHSVDSetDimensions_HSVD));
718: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",DSHSVDGetDimensions_HSVD));
719: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",DSHSVDSetReorthogonalize_HSVD));
720: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",DSHSVDGetReorthogonalize_HSVD));
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }