Actual source code: dshsvd.c
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,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: PetscCallLAPACKInfo("LAPACKsteqr",LAPACKsteqr_("V",&n1,dd+l,ee+l,V+off,&ld,ds->rwork+rwu,&info));
309: for (i=l;i<n;i++) wr[i] = PetscSqrtScalar(PetscAbs(dd[i]));
311: /* Build left singular vectors: U=A*V*Sigma^-1 */
312: PetscCall(PetscArrayzero(U+l*ld,n1*ld));
313: for (i=l;i<n-1;i++) {
314: scal = d[i];
315: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+i,&ld,U+l*ld+i,&ld));
316: j = (i<k)?k:i+1;
317: scal = e[i];
318: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+l*ld+j,&ld,U+l*ld+i,&ld));
319: }
320: scal = d[n-1];
321: PetscCallBLAS("BLASaxpy",BLASaxpy_(&n1,&scal,V+off+(n1-1),&ld,U+off+(n1-1),&ld));
322: /* Multiply by Sigma^-1 */
323: for (i=l;i<n;i++) {scal = 1.0/wr[i]; PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,U+i*ld+l,&one));}
325: } else { /* non-compact */
327: PetscCall(DSAllocateWork_Private(ds,(n+6)*ld,PetscDefined(USE_COMPLEX)?4*ld:ld,2*ld));
328: R = ds->work+swu;
329: swu += n*ld;
330: perm = ds->iwork+iwu;
331: iwu += n;
332: cmplx = ds->iwork+iwu;
333: dd = ds->rwork+rwu;
334: for (j=l;j<m;j++) {
335: for (i=0;i<n;i++) ds->work[i] = Omega[i]*A[i+j*ld];
336: PetscCallBLAS("BLASgemv",BLASgemv_("C",&n,&m,&sone,A,&ld,ds->work,&incx,&szero,V+j*ld,&incx));
337: }
339: /* compute eigenvalues */
340: lwork = (n+6)*ld;
341: #if PetscDefined(USE_COMPLEX)
342: rwu += ld;
343: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,ds->rwork+rwu,&info));
344: #else
345: PetscCallLAPACKInfo("LAPACKsyev",LAPACKsyev_("V","L",&m,V,&ld,dd,ds->work,&lwork,&info));
346: #endif
347: for (i=l;i<PetscMin(n,m);i++) d[i] = PetscSqrtReal(PetscAbsReal(dd[i]));
349: /* Build left singular vectors: U=A*V*Sigma^-1 */
350: for (j=l;j<PetscMin(n,m);j++) {
351: scal = 1.0/d[j];
352: PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&m,&scal,A,&ld,V+j*ld,&incx,&szero,U+j*ld,&incx));
353: }
354: }
356: if (ctx->reorth) { /* Reinforce orthogonality */
357: nv = n1;
358: for (i=0;i<n;i++) cmplx[i] = 0;
359: PetscCall(DSPseudoOrthog_HR(&nv,U+off,ld,Omega+l,R,ld,perm,cmplx,NULL,ds->work+swu));
360: } else { /* Update Omega */
361: for (i=l;i<PetscMin(n,m);i++) Omega[i] = PetscSign(dd[i]);
362: }
364: /* Update projected problem */
365: if (ds->compact) {
366: for (i=l;i<n;i++) d[i] = PetscRealPart(wr[i]);
367: PetscCall(PetscArrayzero(e,n-1));
368: } else {
369: for (i=l;i<m;i++) PetscCall(PetscArrayzero(A+l+i*ld,n-l));
370: for (i=l;i<n;i++) A[i+i*ld] = d[i];
371: }
372: for (i=l;i<PetscMin(n,m);i++) wr[i] = d[i];
373: if (wi) for (i=l;i<PetscMin(n,m);i++) wi[i] = 0.0;
375: if (ctx->reorth) { /* Update vectors V with R */
376: scal = -1.0;
377: for (i=0;i<nv;i++) {
378: if (PetscRealPart(R[i+i*ld]) < 0.0) PetscCallBLAS("BLASscal",BLASscal_(&n1,&scal,V+(i+l)*ld+l,&one));
379: }
380: }
382: if (!ds->compact) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
383: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_U],&U));
384: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_V],&V));
385: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&d));
386: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&Omega));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: #if !PetscDefined(HAVE_MPIUNI)
391: static PetscErrorCode DSSynchronize_HSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
392: {
393: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
394: PetscMPIInt n,rank,off=0,size,ldn,ld3,ld_;
395: PetscScalar *A,*U,*V;
396: PetscReal *T,*D;
398: PetscFunctionBegin;
399: if (ds->compact) kr = 3*ld;
400: else k = (ds->n-l)*ld;
401: kr += ld;
402: if (ds->state>DS_STATE_RAW) k += 2*(ds->n-l)*ld;
403: if (eigr) k += ds->n-l;
404: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
405: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
406: PetscCall(PetscMPIIntCast(ds->n-l,&n));
407: PetscCall(PetscMPIIntCast(ld*(ds->n-l),&ldn));
408: PetscCall(PetscMPIIntCast(3*ld,&ld3));
409: PetscCall(PetscMPIIntCast(ld,&ld_));
410: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
411: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
412: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
413: if (ds->state>DS_STATE_RAW) {
414: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
415: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
416: }
417: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
418: if (!rank) {
419: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
420: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
421: PetscCallMPI(MPI_Pack(D,ld_,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
422: if (ds->state>DS_STATE_RAW) {
423: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
424: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
425: }
426: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
427: }
428: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
429: if (rank) {
430: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
431: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
432: PetscCallMPI(MPI_Unpack(ds->work,size,&off,D,ld_,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
433: if (ds->state>DS_STATE_RAW) {
434: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
435: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
436: }
437: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
438: }
439: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
440: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
441: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
442: if (ds->state>DS_STATE_RAW) {
443: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
444: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
445: }
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
448: #endif
450: static PetscErrorCode DSMatGetSize_HSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
451: {
452: DS_HSVD *ctx = (DS_HSVD*)ds->data;
454: PetscFunctionBegin;
455: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the number of columns with DSHSVDSetDimensions()");
456: switch (t) {
457: case DS_MAT_A:
458: *rows = ds->n;
459: *cols = ds->extrarow? ctx->m+1: ctx->m;
460: break;
461: case DS_MAT_T:
462: *rows = ds->n;
463: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
464: break;
465: case DS_MAT_D:
466: *rows = ds->n;
467: *cols = 1;
468: break;
469: case DS_MAT_U:
470: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
471: *cols = ds->n;
472: break;
473: case DS_MAT_V:
474: *rows = ds->state==DS_STATE_TRUNCATED? ctx->t: ctx->m;
475: *cols = ctx->m;
476: break;
477: default:
478: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
479: }
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: static PetscErrorCode DSHSVDSetDimensions_HSVD(DS ds,PetscInt m)
484: {
485: DS_HSVD *ctx = (DS_HSVD*)ds->data;
487: PetscFunctionBegin;
488: DSCheckAlloc(ds,1);
489: if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
490: ctx->m = ds->ld;
491: } else {
492: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
493: ctx->m = m;
494: }
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*@
499: DSHSVDSetDimensions - Sets the number of columns for a `DSHSVD`.
501: Logically Collective
503: Input Parameters:
504: + ds - the direct solver context
505: - m - the number of columns
507: Notes:
508: This call is complementary to `DSSetDimensions()`, to provide a dimension
509: that is specific to this `DS` type.
511: Level: intermediate
513: .seealso: [](sec:ds), `DSHSVD`, `DSHSVDGetDimensions()`, `DSSetDimensions()`
514: @*/
515: PetscErrorCode DSHSVDSetDimensions(DS ds,PetscInt m)
516: {
517: PetscFunctionBegin;
520: PetscTryMethod(ds,"DSHSVDSetDimensions_C",(DS,PetscInt),(ds,m));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode DSHSVDGetDimensions_HSVD(DS ds,PetscInt *m)
525: {
526: DS_HSVD *ctx = (DS_HSVD*)ds->data;
528: PetscFunctionBegin;
529: *m = ctx->m;
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@
534: DSHSVDGetDimensions - Returns the number of columns for a `DSHSVD`.
536: Not Collective
538: Input Parameter:
539: . ds - the direct solver context
541: Output Parameter:
542: . m - the number of columns
544: Level: intermediate
546: .seealso: [](sec:ds), `DSHSVD`, `DSHSVDSetDimensions()`
547: @*/
548: PetscErrorCode DSHSVDGetDimensions(DS ds,PetscInt *m)
549: {
550: PetscFunctionBegin;
552: PetscAssertPointer(m,2);
553: PetscUseMethod(ds,"DSHSVDGetDimensions_C",(DS,PetscInt*),(ds,m));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode DSHSVDSetReorthogonalize_HSVD(DS ds,PetscBool reorth)
558: {
559: DS_HSVD *ctx = (DS_HSVD*)ds->data;
561: PetscFunctionBegin;
562: ctx->reorth = reorth;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*@
567: DSHSVDSetReorthogonalize - Sets the reorthogonalization of the left vectors in a `DSHSVD`.
569: Logically Collective
571: Input Parameters:
572: + ds - the direct solver context
573: - reorth - the reorthogonalization flag
575: Options Database Key:
576: . -ds_hsvd_reorthog (true|false) - sets the reorthogonalization flag
578: Note:
579: The computed left vectors (`U`) should be orthogonal with respect to the signature (`D`).
580: But it may be necessary to enforce this with a final reorthogonalization step (omitted
581: by default).
583: Level: intermediate
585: .seealso: [](sec:ds), `DSHSVD`, `DSHSVDGetReorthogonalize()`
586: @*/
587: PetscErrorCode DSHSVDSetReorthogonalize(DS ds,PetscBool reorth)
588: {
589: PetscFunctionBegin;
592: PetscTryMethod(ds,"DSHSVDSetReorthogonalize_C",(DS,PetscBool),(ds,reorth));
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: static PetscErrorCode DSHSVDGetReorthogonalize_HSVD(DS ds,PetscBool *reorth)
597: {
598: DS_HSVD *ctx = (DS_HSVD*)ds->data;
600: PetscFunctionBegin;
601: *reorth = ctx->reorth;
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@
606: DSHSVDGetReorthogonalize - Returns the reorthogonalization flag of a `DSHSVD`.
608: Not Collective
610: Input Parameter:
611: . ds - the direct solver context
613: Output Parameter:
614: . reorth - the reorthogonalization flag
616: Level: intermediate
618: .seealso: [](sec:ds), `DSHSVD`, `DSHSVDSetReorthogonalize()`
619: @*/
620: PetscErrorCode DSHSVDGetReorthogonalize(DS ds,PetscBool *reorth)
621: {
622: PetscFunctionBegin;
624: PetscAssertPointer(reorth,2);
625: PetscUseMethod(ds,"DSHSVDGetReorthogonalize_C",(DS,PetscBool*),(ds,reorth));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: static PetscErrorCode DSSetFromOptions_HSVD(DS ds,PetscOptionItems PetscOptionsObject)
630: {
631: PetscBool flg,reorth;
633: PetscFunctionBegin;
634: PetscOptionsHeadBegin(PetscOptionsObject,"DS HSVD Options");
636: PetscCall(PetscOptionsBool("-ds_hsvd_reorthog","Reorthogonalize U vectors","DSHSVDSetReorthogonalize",PETSC_FALSE,&reorth,&flg));
637: if (flg) PetscCall(DSHSVDSetReorthogonalize(ds,reorth));
639: PetscOptionsHeadEnd();
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode DSDestroy_HSVD(DS ds)
644: {
645: PetscFunctionBegin;
646: PetscCall(PetscFree(ds->data));
647: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",NULL));
648: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",NULL));
649: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",NULL));
650: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",NULL));
651: PetscFunctionReturn(PETSC_SUCCESS);
652: }
654: static PetscErrorCode DSSetCompact_HSVD(DS ds,PetscBool comp)
655: {
656: PetscFunctionBegin;
657: if (!comp) PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: static PetscErrorCode DSReallocate_HSVD(DS ds,PetscInt ld)
662: {
663: PetscInt i,*perm=ds->perm;
665: PetscFunctionBegin;
666: for (i=0;i<DS_NUM_MAT;i++) {
667: if (!ds->compact && i==DS_MAT_A) continue;
668: if (i!=DS_MAT_U && i!=DS_MAT_V && i!=DS_MAT_T && i!=DS_MAT_D) PetscCall(MatDestroy(&ds->omat[i]));
669: }
671: if (!ds->compact) PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
672: PetscCall(DSReallocateMat_Private(ds,DS_MAT_U,ld));
673: PetscCall(DSReallocateMat_Private(ds,DS_MAT_V,ld));
674: PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
675: PetscCall(DSReallocateMat_Private(ds,DS_MAT_D,ld));
677: PetscCall(PetscMalloc1(ld,&ds->perm));
678: PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
679: PetscCall(PetscFree(perm));
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: /*MC
684: DSHSVD - Dense Hyperbolic Singular Value Decomposition.
686: Notes:
687: The problem is expressed as $A = U\Sigma V^*$, where $A$ is rectangular in
688: general, with $n$ rows and $m$ columns. $U$ is orthogonal with respect to a
689: signature matrix $\Omega$, stored in $D$, $V$ is orthogonal, $\Sigma$ is a diagonal
690: matrix whose diagonal elements are the arguments of `DSSolve()`. After
691: solve, $A$ is overwritten with $\Sigma$, $D$ is overwritten with the new signature.
693: The matrices of left and right singular vectors, $U$ and $V$, have size $n$ and $m$,
694: respectively. The number of columns m must be specified via `DSHSVDSetDimensions()`.
696: If the `DS` object is in the intermediate state, $A$ is assumed to be in upper
697: bidiagonal form (possibly with an arrow) and is stored in compact format
698: on matrix $T$. The compact storage is implemented for the square case
699: only, $m=n$. The extra row should be interpreted in this case as an extra column.
701: Used DS matrices:
702: + `DS_MAT_A` - problem matrix (used only if `compact=PETSC_FALSE`)
703: . `DS_MAT_T` - upper bidiagonal matrix
704: . `DS_MAT_D` - diagonal matrix (signature)
705: . `DS_MAT_U` - left singular vectors
706: - `DS_MAT_V` - right singular vectors
708: Implemented methods:
709: . 0 - Cross product $A^*\Omega A$
711: Level: beginner
713: .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSHSVDSetDimensions()`, `DSSetCompact()`
714: M*/
715: SLEPC_EXTERN PetscErrorCode DSCreate_HSVD(DS ds)
716: {
717: DS_HSVD *ctx;
719: PetscFunctionBegin;
720: PetscCall(PetscNew(&ctx));
721: ds->data = (void*)ctx;
723: ds->ops->allocate = DSAllocate_HSVD;
724: ds->ops->setfromoptions = DSSetFromOptions_HSVD;
725: ds->ops->view = DSView_HSVD;
726: ds->ops->vectors = DSVectors_HSVD;
727: ds->ops->solve[0] = DSSolve_HSVD_CROSS;
728: ds->ops->sort = DSSort_HSVD;
729: ds->ops->truncate = DSTruncate_HSVD;
730: ds->ops->update = DSUpdateExtraRow_HSVD;
731: ds->ops->destroy = DSDestroy_HSVD;
732: ds->ops->matgetsize = DSMatGetSize_HSVD;
733: #if !PetscDefined(HAVE_MPIUNI)
734: ds->ops->synchronize = DSSynchronize_HSVD;
735: #endif
736: ds->ops->setcompact = DSSetCompact_HSVD;
737: ds->ops->reallocate = DSReallocate_HSVD;
738: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetDimensions_C",DSHSVDSetDimensions_HSVD));
739: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetDimensions_C",DSHSVDGetDimensions_HSVD));
740: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDSetReorthogonalize_C",DSHSVDSetReorthogonalize_HSVD));
741: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSHSVDGetReorthogonalize_C",DSHSVDGetReorthogonalize_HSVD));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }