Actual source code: dshsvd.c
slepc-3.21.1 2024-04-26
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: 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: 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: iwu += n;
290: dd = ds->rwork+rwu;
291: rwu += ld;
292: ee = ds->rwork+rwu;
293: rwu += ld;
294: for (i=0;i<l;i++) {dd[i] = d[i]*d[i]*Omega[i]; ee[i] = 0.0;}
295: for (i=l;i<=ds->k;i++) {
296: dd[i] = Omega[i]*d[i]*d[i];
297: ee[i] = Omega[i]*d[i]*e[i];
298: }
299: for (i=l;i<k;i++) dd[k] += Omega[i]*e[i]*e[i];
300: for (i=k+1;i<n;i++) {
301: dd[i] = Omega[i]*d[i]*d[i]+Omega[i-1]*e[i-1]*e[i-1];
302: ee[i] = Omega[i]*d[i]*e[i];
303: }
305: /* Reduce T to tridiagonal form */
306: PetscCall(DSArrowTridiag(n2,dd+l,ee+l,V+off,ld));
308: /* Solve the tridiagonal eigenproblem corresponding to T */
309: PetscCallBLAS("LAPACKsteqr",LAPACKsteqr_("V",&n1,dd+l,ee+l,V+off,&ld,ds->rwork+rwu,&info));
310: SlepcCheckLapackInfo("steqr",info);
311: for (i=l;i<n;i++) wr[i] = PetscSqrtScalar(PetscAbs(dd[i]));
313: /* Build left singular vectors: U=A*V*Sigma^-1 */
314: PetscCall(PetscArrayzero(U+l*ld,n1*ld));
315: for (i=l;i<n-1;i++) {
316: scal = d[i];
317: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+i,&ld,U+l*ld+i,&ld));
318: j = (i<k)?k:i+1;
319: scal = e[i];
320: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+j,&ld,U+l*ld+i,&ld));
321: }
322: scal = d[n-1];
323: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+off+(n1-1),&ld,U+off+(n1-1),&ld));
324: /* Multiply by Sigma^-1 */
325: for (i=l;i<n;i++) {scal = 1.0/wr[i]; PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,U+i*ld+l,&one));}
327: } else { /* non-compact */
329: PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,PetscDefined(USE_COMPLEX)?4*ld:ld,2*ld));
330: R = ds->work+swu;
331: swu += n*ld;
332: perm = ds->iwork+iwu;
333: iwu += n;
334: cmplx = ds->iwork+iwu;
335: iwu += n;
336: dd = ds->rwork+rwu;
337: rwu += ld;
338: for (j=l;j<m;j++) {
339: for (i=0;i<n;i++) ds->work[i] = Omega[i]*A[i+j*ld];
340: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&m,&sone,A,&ld,ds->work,&incx,&szero,V+j*ld,&incx));
341: }
343: /* compute eigenvalues */
344: lwork = (n+6)*ld;
345: #if defined(PETSC_USE_COMPLEX)
346: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,ds->rwork+rwu,&info));
347: #else
348: PetscCallBLAS("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,&info));
349: #endif
350: SlepcCheckLapackInfo("syev",info);
351: for (i=l;i<PetscMin(n,m);i++) d[i] = PetscSqrtReal(PetscAbsReal(dd[i]));
353: /* Build left singular vectors: U=A*V*Sigma^-1 */
354: for (j=l;j<PetscMin(n,m);j++) {
355: scal = 1.0/d[j];
356: PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&m,&scal,A,&ld,V+j*ld,&incx,&szero,U+j*ld,&incx));
357: }
358: }
360: if (ctx->reorth) { /* Reinforce orthogonality */
361: nv = n1;
362: for (i=0;i<n;i++) cmplx[i] = 0;
363: PetscCall(DSPseudoOrthog_HR(&nv,U+off,ld,Omega+l,R,ld,perm,cmplx,NULL,ds->work+swu));
364: } else { /* Update Omega */
365: for (i=l;i<PetscMin(n,m);i++) Omega[i] = PetscSign(dd[i]);
366: }
368: /* Update projected problem */
369: if (ds->compact) {
370: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
371: PetscCall(PetscArrayzero(e,n-1));
372: } else {
373: for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
374: for (i=l;i<n;i++) A[i+i*ld] = d[i];
375: }
376: for (i=l;i<PetscMin(n,m);i++) wr[i] = d[i];
377: if (wi) for (i=l;i<PetscMin(n,m);i++) wi[i] = 0.0;
379: if (ctx->reorth) { /* Update vectors V with R */
380: scal = -1.0;
381: for (i=0;i<nv;i++) {
382: if (PetscRealPart(R[i+i*ld]) < 0.0) PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,V+(i+l)*ld+l,&one));
383: }
384: }
386: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
387: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
388: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
389: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
390: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: #if !defined(PETSC_HAVE_MPIUNI)
395: static PetscErrorCode DSSynchronize_HSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
396: {
397: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
398: PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
399: PetscScalar *A,*U,*V;
400: PetscReal *T,*D;
402: PetscFunctionBegin;
403: if (ds->compact) kr = 3*ld;
404: else k = (ds->n-l)*ld;
405: kr += ld;
406: if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
407: if (eigr) k += ds->n-l;
408: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
409: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
410: PetscCall(PetscMPIIntCast(ds->n-l,&n));
411: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
412: PetscCall(PetscMPIIntCast(3*ld,&ld3));
413: PetscCall(PetscMPIIntCast(ld,&ld_));
414: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
415: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
416: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
417: if (ds->state>DS_STATE_RAW) {
418: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
419: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
420: }
421: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
422: if (!rank) {
423: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
424: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
425: PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
426: if (ds->state>DS_STATE_RAW) {
427: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
428: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
429: }
430: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
431: }
432: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
433: if (rank) {
434: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
435: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
436: PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
437: if (ds->state>DS_STATE_RAW) {
438: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
439: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
440: }
441: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
442: }
443: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
444: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
445: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
446: if (ds->state>DS_STATE_RAW) {
447: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
448: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
449: }
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
452: #endif
454: static PetscErrorCode DSMatGetSize_HSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
455: {
456: DS_HSVD *ctx = (DS_HSVD*)ds->data;
458: PetscFunctionBegin;
459: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
460: switch (t) {
461: case DS_MAT_A:
462: *rows = ds->n;
463: *cols = ds->extrarow? ctx->m+1: ctx->m;
464: break;
465: case DS_MAT_T:
466: *rows = ds->n;
467: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
468: break;
469: case DS_MAT_D:
470: *rows = ds->n;
471: *cols = 1;
472: break;
473: case DS_MAT_U:
474: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
475: *cols = ds->n;
476: break;
477: case DS_MAT_V:
478: *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
479: *cols = ctx->m;
480: break;
481: default:
482: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
483: }
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: static PetscErrorCode DSHSVDSetDimensions_HSVD(DS ds,PetscInt m)
488: {
489: DS_HSVD *ctx = (DS_HSVD*)ds->data;
491: PetscFunctionBegin;
492: DSCheckAlloc(ds,1);
493: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
494: ctx->m = ds->ld;
495: } else {
496: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
497: ctx->m = m;
498: }
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*@
503: DSHSVDSetDimensions - Sets the number of columns for a DSHSVD.
505: Logically Collective
507: Input Parameters:
508: + ds - the direct solver context
509: - m - the number of columns
511: Notes:
512: This call is complementary to DSSetDimensions(), to provide a dimension
513: that is specific to this DS type.
515: Level: intermediate
517: .seealso: DSHSVDGetDimensions(), DSSetDimensions()
518: @*/
519: PetscErrorCode DSHSVDSetDimensions(DS ds,PetscInt m)
520: {
521: PetscFunctionBegin;
524: PetscTryMethod(ds,"DSHSVDSetDimensions_C",(DS,PetscInt),(ds,m));
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode DSHSVDGetDimensions_HSVD(DS ds,PetscInt *m)
529: {
530: DS_HSVD *ctx = (DS_HSVD*)ds->data;
532: PetscFunctionBegin;
533: *m = ctx->m;
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*@
538: DSHSVDGetDimensions - Returns the number of columns for a DSHSVD.
540: Not Collective
542: Input Parameter:
543: . ds - the direct solver context
545: Output Parameters:
546: . m - the number of columns
548: Level: intermediate
550: .seealso: DSHSVDSetDimensions()
551: @*/
552: PetscErrorCode DSHSVDGetDimensions(DS ds,PetscInt *m)
553: {
554: PetscFunctionBegin;
556: PetscAssertPointer(m,2);
557: PetscUseMethod(ds,"DSHSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode DSHSVDSetReorthogonalize_HSVD(DS ds,PetscBool reorth)
562: {
563: DS_HSVD *ctx = (DS_HSVD*)ds->data;
565: PetscFunctionBegin;
566: ctx->reorth = reorth;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@
571: DSHSVDSetReorthogonalize - Sets the reorthogonalization of the left vectors in a DSHSVD.
573: Logically Collective
575: Input Parameters:
576: + ds - the direct solver context
577: - reorth - the reorthogonalization flag
579: Options Database Key:
580: . -ds_hsvd_reorthog <bool> - sets the reorthogonalization flag
582: Note:
583: The computed left vectors (U) should be orthogonal with respect to the signature (D).
584: But it may be necessary to enforce this with a final reorthogonalization step (omitted
585: by default).
587: Level: intermediate
589: .seealso: DSHSVDGetReorthogonalize()
590: @*/
591: PetscErrorCode DSHSVDSetReorthogonalize(DS ds,PetscBool reorth)
592: {
593: PetscFunctionBegin;
596: PetscTryMethod(ds,"DSHSVDSetReorthogonalize_C",(DS,PetscBool),(ds,reorth));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode DSHSVDGetReorthogonalize_HSVD(DS ds,PetscBool *reorth)
601: {
602: DS_HSVD *ctx = (DS_HSVD*)ds->data;
604: PetscFunctionBegin;
605: *reorth = ctx->reorth;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: DSHSVDGetReorthogonalize - Returns the reorthogonalization flag of a DSHSVD.
612: Not Collective
614: Input Parameter:
615: . ds - the direct solver context
617: Output Parameters:
618: . reorth - the reorthogonalization flag
620: Level: intermediate
622: .seealso: DSHSVDSetReorthogonalize()
623: @*/
624: PetscErrorCode DSHSVDGetReorthogonalize(DS ds,PetscBool *reorth)
625: {
626: PetscFunctionBegin;
628: PetscAssertPointer(reorth,2);
629: PetscUseMethod(ds,"DSHSVDGetReorthogonalize_C",(DS,PetscBool*),(ds,reorth));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: static PetscErrorCode DSSetFromOptions_HSVD(DS ds,PetscOptionItems *PetscOptionsObject)
634: {
635: PetscBool flg,reorth;
637: PetscFunctionBegin;
638: PetscOptionsHeadBegin(PetscOptionsObject,"DS HSVD Options");
640: PetscCall(PetscOptionsBool("-ds_hsvd_reorthog","Reorthogonalize U vectors","DSHSVDSetReorthogonalize",PETSC_FALSE,&reorth,&flg));
641: if (flg) PetscCall(DSHSVDSetReorthogonalize(ds,reorth));
643: PetscOptionsHeadEnd();
644: PetscFunctionReturn(PETSC_SUCCESS);
645: }
647: static PetscErrorCode DSDestroy_HSVD(DS ds)
648: {
649: PetscFunctionBegin;
650: PetscCall(PetscFree(ds->data));
651: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",NULL));
652: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",NULL));
653: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",NULL));
654: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",NULL));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: /*MC
659: DSHSVD - Dense Hyperbolic Singular Value Decomposition.
661: Level: beginner
663: Notes:
664: The problem is expressed as A = U*Sigma*V', where A is rectangular in
665: general, with n rows and m columns. U is orthogonal with respect to a
666: signature matrix, stored in D. V is orthogonal. Sigma is a diagonal
667: matrix whose diagonal elements are the arguments of DSSolve(). After
668: solve, A is overwritten with Sigma, D is overwritten with the new signature.
670: The matrices of left and right singular vectors, U and V, have size n and m,
671: respectively. The number of columns m must be specified via DSHSVDSetDimensions().
673: If the DS object is in the intermediate state, A is assumed to be in upper
674: bidiagonal form (possibly with an arrow) and is stored in compact format
675: on matrix T. The compact storage is implemented for the square case
676: only, m=n. The extra row should be interpreted in this case as an extra column.
678: Used DS matrices:
679: + DS_MAT_A - problem matrix
680: . DS_MAT_T - upper bidiagonal matrix
681: - DS_MAT_D - diagonal matrix (signature)
683: Implemented methods:
684: . 0 - Cross product A'*Omega*A
686: .seealso: DSCreate(), DSSetType(), DSType, DSHSVDSetDimensions()
687: M*/
688: SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS ds)
689: {
690: DS_HSVD *ctx;
692: PetscFunctionBegin;
693: PetscCall(PetscNew(&ctx));
694: ds->data = (void*)ctx;
696: ds->ops->allocate = DSAllocate_HSVD;
697: ds->ops->setfromoptions = DSSetFromOptions_HSVD;
698: ds->ops->view = DSView_HSVD;
699: ds->ops->vectors = DSVectors_HSVD;
700: ds->ops->solve[0] = DSSolve_HSVD_CROSS;
701: ds->ops->sort = DSSort_HSVD;
702: #if !defined(PETSC_HAVE_MPIUNI)
703: ds->ops->synchronize = DSSynchronize_HSVD;
704: #endif
705: ds->ops->truncate = DSTruncate_HSVD;
706: ds->ops->update = DSUpdateExtraRow_HSVD;
707: ds->ops->destroy = DSDestroy_HSVD;
708: ds->ops->matgetsize = DSMatGetSize_HSVD;
709: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",DSHSVDSetDimensions_HSVD));
710: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",DSHSVDGetDimensions_HSVD));
711: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",DSHSVDSetReorthogonalize_HSVD));
712: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",DSHSVDGetReorthogonalize_HSVD));
713: PetscFunctionReturn(PETSC_SUCCESS);
714: }