Actual source code: dsgsvd.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: */
10: #include <slepc/private/dsimpl.h>
11: #include <slepcblaslapack.h>
13: typedef struct {
14: PetscInt m; /* number of columns */
15: PetscInt p; /* number of rows of B */
16: PetscInt tm; /* number of rows of X after truncating */
17: PetscInt tp; /* number of rows of V after truncating */
18: } DS_GSVD;
20: static PetscErrorCode DSAllocate_GSVD(DS ds,PetscInt ld)
21: {
22: PetscFunctionBegin;
23: PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
24: PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
25: PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
26: PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
27: PetscCall(DSAllocateMat_Private(ds,DS_MAT_V));
28: PetscCall(DSAllocateMat_Private(ds,DS_MAT_T));
29: PetscCall(DSAllocateMat_Private(ds,DS_MAT_D));
30: PetscCall(PetscFree(ds->perm));
31: PetscCall(PetscMalloc1(ld,&ds->perm));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: /*
36: In compact form, A is either in form (a) or (b):
38: (a) (b)
39: lower bidiagonal with upper arrow (n=m+1) square upper bidiagonal with upper arrow (n=m)
40: 0 l k m-1
41: ----------------------------------------- 0 l k m-1
42: |* . | -----------------------------------------
43: | * . | |* . |
44: | * . | | * . |
45: | * . | | * . |
46: l |. . . . o o | l |. . . o o |
47: | o o | | o o |
48: | o o | | o o |
49: | o o | | o o |
50: | o o | | o o |
51: | o o | | o o |
52: k |. . . . . . . . . . o | k |. . . . . . . . . o x |
53: | x x | | x x |
54: | x x | | x x |
55: | x x | | x x |
56: | x x | | x x |
57: | x x | | x x |
58: | x x | | x x |
59: | x x | | x x |
60: | x x | | x x |
61: | x x| | x x|
62: n-1 | x| n-1 | x|
63: ----------------------------------------- -----------------------------------------
65: and B is square bidiagonal with upper arrow (p=m)
67: 0 l k m-1
68: -----------------------------------------
69: |* . |
70: | * . |
71: | * . |
72: | * . |
73: l |. . . . o o |
74: | o o |
75: | o o |
76: | o o |
77: | o o |
78: | o o |
79: k |. . . . . . . . . . o x |
80: | x x |
81: | x x |
82: | x x |
83: | x x |
84: | x x |
85: | x x |
86: | x x |
87: | x x|
88: p-1 | x|
89: ----------------------------------------
90: */
91: static PetscErrorCode DSView_GSVD(DS ds,PetscViewer viewer)
92: {
93: DS_GSVD *ctx = (DS_GSVD*)ds->data;
94: PetscViewerFormat format;
95: PetscInt i,j,r,k=ds->k,n=ds->n,m=ctx->m,p=ctx->p,rowsa,rowsb,colsa,colsb;
96: PetscReal *T,*S,value;
98: PetscFunctionBegin;
99: PetscCall(PetscViewerGetFormat(viewer,&format));
100: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
101: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
102: PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns: %" PetscInt_FMT "\n",m));
103: PetscCall(PetscViewerASCIIPrintf(viewer,"number of rows of B: %" PetscInt_FMT "\n",p));
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
106: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
107: if (ds->compact) {
108: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
109: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&S));
110: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE));
111: rowsa = n;
112: colsa = ds->extrarow? m+1: m;
113: rowsb = p;
114: colsb = ds->extrarow? m+1: m;
115: if (format == PETSC_VIEWER_ASCII_MATLAB) {
116: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsa,colsa));
117: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
118: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
119: for (i=0;i<PetscMin(rowsa,colsa);i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)T[i]));
120: for (i=0;i<k;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,k+1,(double)T[i+ds->ld]));
121: if (n>m) { /* A lower bidiagonal */
122: for (i=k;i<rowsa-1;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+2,i+1,(double)T[i+ds->ld]));
123: } else { /* A (square) upper bidiagonal */
124: for (i=k;i<colsa-1;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+2,(double)T[i+ds->ld]));
125: }
126: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_T]));
127: PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT "\n",rowsb,colsb));
128: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",2*ds->n));
129: PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n"));
130: for (i=0;i<rowsb;i++) PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,i+1,(double)S[i]));
131: for (i=0;i<colsb-1;i++) {
132: r = PetscMax(i+2,ds->k+1);
133: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,r,(double)T[i+2*ds->ld]));
134: }
135: PetscCall(PetscViewerASCIIPrintf(viewer,"];\n%s = spconvert(zzz);\n",DSMatName[DS_MAT_D]));
136: } else {
137: PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_T]));
138: for (i=0;i<rowsa;i++) {
139: for (j=0;j<colsa;j++) {
140: if (i==j) value = T[i];
141: else if (i<ds->k && j==ds->k) value = T[i+ds->ld];
142: else if (n>m && i==j+1 && i>ds->k) value = T[j+ds->ld];
143: else if (n<=m && i+1==j && i>=ds->k) value = T[i+ds->ld];
144: else value = 0.0;
145: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
146: }
147: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
148: }
149: PetscCall(PetscViewerASCIIPrintf(viewer,"Matrix %s =\n",DSMatName[DS_MAT_D]));
150: for (i=0;i<rowsb;i++) {
151: for (j=0;j<colsb;j++) {
152: if (i==j) value = S[i];
153: else if (i<ds->k && j==ds->k) value = T[PetscMin(i,j)+2*ds->ld];
154: else if (i+1==j && i>=ds->k) value = T[i+2*ds->ld];
155: else value = 0.0;
156: PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)value));
157: }
158: PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
159: }
160: }
161: PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE));
162: PetscCall(PetscViewerFlush(viewer));
163: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
164: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&S));
165: } else {
166: PetscCall(DSViewMat(ds,viewer,DS_MAT_A));
167: PetscCall(DSViewMat(ds,viewer,DS_MAT_B));
168: }
169: if (ds->state>DS_STATE_INTERMEDIATE) {
170: PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
171: PetscCall(DSViewMat(ds,viewer,DS_MAT_U));
172: PetscCall(DSViewMat(ds,viewer,DS_MAT_V));
173: }
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode DSVectors_GSVD(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
178: {
179: PetscFunctionBegin;
180: switch (mat) {
181: case DS_MAT_U:
182: case DS_MAT_V:
183: if (rnorm) *rnorm = 0.0;
184: break;
185: case DS_MAT_X:
186: break;
187: default:
188: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode DSSort_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *k)
194: {
195: DS_GSVD *ctx = (DS_GSVD*)ds->data;
196: PetscInt t,l,ld=ds->ld,i,*perm,*perm2;
197: PetscReal *T=NULL,*D=NULL,*eig;
198: PetscScalar *A=NULL,*B=NULL;
199: PetscBool compact=ds->compact;
201: PetscFunctionBegin;
202: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
203: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
204: l = ds->l;
205: t = ds->t;
206: perm = ds->perm;
207: PetscCall(PetscMalloc2(t,&eig,t,&perm2));
208: if (compact) {
209: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
210: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
211: for (i=0;i<t;i++) eig[i] = (D[i]==0)?PETSC_INFINITY:T[i]/D[i];
212: } else {
213: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
214: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
215: for (i=0;i<t;i++) eig[i] = (B[i+i*ld]==0)?PETSC_INFINITY:PetscRealPart(A[i+i*ld])/PetscRealPart(B[i*(1+ld)]);
216: }
217: PetscCall(DSSortEigenvaluesReal_Private(ds,eig,perm));
218: PetscCall(PetscArraycpy(perm2,perm,t));
219: for (i=l;i<t;i++) wr[i] = eig[perm[i]];
220: if (compact) {
221: PetscCall(PetscArraycpy(eig,T,t));
222: for (i=l;i<t;i++) T[i] = eig[perm[i]];
223: PetscCall(PetscArraycpy(eig,D,t));
224: for (i=l;i<t;i++) D[i] = eig[perm[i]];
225: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
226: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
227: } else {
228: for (i=l;i<t;i++) eig[i] = PetscRealPart(A[i*(1+ld)]);
229: for (i=l;i<t;i++) A[i*(1+ld)] = eig[perm[i]];
230: for (i=l;i<t;i++) eig[i] = PetscRealPart(B[i*(1+ld)]);
231: for (i=l;i<t;i++) B[i*(1+ld)] = eig[perm[i]];
232: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
233: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
234: }
235: PetscCall(DSPermuteColumns_Private(ds,l,t,ds->n,DS_MAT_U,perm2));
236: PetscCall(PetscArraycpy(perm2,perm,t));
237: PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->m,DS_MAT_X,perm2));
238: PetscCall(DSPermuteColumns_Private(ds,l,t,ctx->p,DS_MAT_V,perm));
239: PetscCall(PetscFree2(eig,perm2));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: static PetscErrorCode DSUpdateExtraRow_GSVD(DS ds)
244: {
245: DS_GSVD *ctx = (DS_GSVD*)ds->data;
246: PetscInt i;
247: PetscBLASInt n=0,m=0,ld=0;
248: const PetscScalar *U,*V;
249: PetscReal *T,*e,*f,alpha,beta,betah;
251: PetscFunctionBegin;
252: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
253: PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
254: PetscCall(PetscBLASIntCast(ds->n,&n));
255: PetscCall(PetscBLASIntCast(ctx->m,&m));
256: PetscCall(PetscBLASIntCast(ds->ld,&ld));
257: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
258: e = T+ld;
259: f = T+2*ld;
260: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_U],&U));
261: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_V],&V));
262: if (n<=m) { /* upper variant, A is square upper bidiagonal */
263: beta = e[m-1]; /* in compact, we assume all entries are zero except the last one */
264: betah = f[m-1];
265: for (i=0;i<m;i++) {
266: e[i] = PetscRealPart(beta*U[m-1+i*ld]);
267: f[i] = PetscRealPart(betah*V[m-1+i*ld]);
268: }
269: } else { /* lower variant, A is (m+1)xm lower bidiagonal */
270: alpha = T[m];
271: betah = f[m-1];
272: for (i=0;i<m;i++) {
273: e[i] = PetscRealPart(alpha*U[m+i*ld]);
274: f[i] = PetscRealPart(betah*V[m-1+i*ld]);
275: }
276: T[m] = PetscRealPart(alpha*U[m+m*ld]);
277: }
278: ds->k = m;
279: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_U],&U));
280: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_V],&V));
281: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
282: PetscFunctionReturn(PETSC_SUCCESS);
283: }
285: static PetscErrorCode DSTruncate_GSVD(DS ds,PetscInt n,PetscBool trim)
286: {
287: DS_GSVD *ctx = (DS_GSVD*)ds->data;
288: PetscScalar *U;
289: PetscReal *T;
290: PetscInt i,m=ctx->m,ld=ds->ld;
291: PetscBool lower=(ds->n>ctx->m)?PETSC_TRUE:PETSC_FALSE;
293: PetscFunctionBegin;
294: PetscCheck(ds->compact,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented for non-compact storage");
295: if (trim) {
296: ds->l = 0;
297: ds->k = 0;
298: ds->n = lower? n+1: n;
299: ctx->m = n;
300: ctx->p = n;
301: ds->t = ds->n; /* truncated length equal to the new dimension */
302: ctx->tm = ctx->m; /* must also keep the previous dimension of X */
303: ctx->tp = ctx->p; /* must also keep the previous dimension of V */
304: } else {
305: if (lower) {
306: /* move value of diagonal element of arrow (alpha) */
307: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
308: T[n] = T[m];
309: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
310: /* copy last column of U so that it updates the next initial vector of U1 */
311: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
312: for (i=0;i<=m;i++) U[i+n*ld] = U[i+m*ld];
313: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
314: }
315: ds->k = ds->extrarow? n: 0;
316: ds->t = ds->n; /* truncated length equal to previous dimension */
317: ctx->tm = ctx->m; /* must also keep the previous dimension of X */
318: ctx->tp = ctx->p; /* must also keep the previous dimension of V */
319: ds->n = lower? n+1: n;
320: ctx->m = n;
321: ctx->p = n;
322: }
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode DSSwitchFormat_GSVD(DS ds)
327: {
328: DS_GSVD *ctx = (DS_GSVD*)ds->data;
329: PetscReal *T,*D;
330: PetscScalar *A,*B;
331: PetscInt i,n=ds->n,k=ds->k,ld=ds->ld,m=ctx->m;
333: PetscFunctionBegin;
334: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
335: /* switch from compact (arrow) to dense storage */
336: /* bidiagonal associated to B is stored in D and T+2*ld */
337: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_A],&A));
338: PetscCall(MatDenseGetArrayWrite(ds->omat[DS_MAT_B],&B));
339: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
340: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
341: PetscCall(PetscArrayzero(A,ld*ld));
342: PetscCall(PetscArrayzero(B,ld*ld));
343: for (i=0;i<k;i++) {
344: A[i+i*ld] = T[i];
345: A[i+k*ld] = T[i+ld];
346: B[i+i*ld] = D[i];
347: B[i+k*ld] = T[i+2*ld];
348: }
349: /* B is upper bidiagonal */
350: B[k+k*ld] = D[k];
351: for (i=k+1;i<m;i++) {
352: B[i+i*ld] = D[i];
353: B[i-1+i*ld] = T[i-1+2*ld];
354: }
355: /* A can be upper (square) or lower bidiagonal */
356: for (i=k;i<m;i++) A[i+i*ld] = T[i];
357: if (n>m) for (i=k;i<m;i++) A[i+1+i*ld] = T[i+ld];
358: else for (i=k+1;i<m;i++) A[i-1+i*ld] = T[i-1+ld];
359: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_A],&A));
360: PetscCall(MatDenseRestoreArrayWrite(ds->omat[DS_MAT_B],&B));
361: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
362: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*
367: Compact format is used when [A;B] has orthonormal columns.
368: In this case R=I and the GSVD of (A,B) is the CS decomposition
369: */
370: static PetscErrorCode DSSolve_GSVD(DS ds,PetscScalar *wr,PetscScalar *wi)
371: {
372: DS_GSVD *ctx = (DS_GSVD*)ds->data;
373: PetscInt i,j;
374: PetscBLASInt n1,m1,lc = 0,n = 0,m = 0,p = 0,p1,l,k,q,ld,off,lwork,r;
375: PetscScalar *A,*B,*X,*U,*V,sone=1.0,smone=-1.0;
376: PetscReal *alpha,*beta,*T,*D;
377: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
378: PetscScalar a,dummy;
379: PetscReal rdummy;
380: PetscBLASInt idummy;
381: #endif
383: PetscFunctionBegin;
384: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
385: PetscCall(PetscBLASIntCast(ds->n,&m));
386: PetscCall(PetscBLASIntCast(ctx->m,&n));
387: PetscCall(PetscBLASIntCast(ctx->p,&p));
388: PetscCall(PetscBLASIntCast(ds->l,&lc));
389: PetscCheck(ds->compact || lc==0,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"DSGSVD with non-compact format does not support locking");
390: /* In compact storage B is always nxn and A can be either nxn or (n+1)xn */
391: PetscCheck(!ds->compact || (p==n && (m==p || m==p+1)),PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Dimensions not supported in compact format");
392: PetscCall(PetscBLASIntCast(ds->ld,&ld));
393: n1 = n-lc; /* n1 = size of leading block, excl. locked + size of trailing block */
394: m1 = m-lc;
395: p1 = p-lc;
396: off = lc+lc*ld;
397: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
398: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
399: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
400: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
401: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
402: PetscCall(PetscArrayzero(X,ld*ld));
403: for (i=0;i<lc;i++) X[i+i*ld] = 1.0;
404: PetscCall(PetscArrayzero(U,ld*ld));
405: for (i=0;i<lc;i++) U[i+i*ld] = 1.0;
406: PetscCall(PetscArrayzero(V,ld*ld));
407: for (i=0;i<lc;i++) V[i+i*ld] = 1.0;
408: if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
410: #if !defined(SLEPC_MISSING_LAPACK_GGSVD3)
411: /* workspace query and memory allocation */
412: lwork = -1;
413: #if !PetscDefined(USE_COMPLEX)
414: PetscCallLAPACKInfo("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,&dummy,&ld,&dummy,&ld,&rdummy,&rdummy,&dummy,&ld,&dummy,&ld,&dummy,&ld,&a,&lwork,&idummy,&info));
415: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
416: #else
417: PetscCallLAPACKInfo("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,&dummy,&ld,&dummy,&ld,&rdummy,&rdummy,&dummy,&ld,&dummy,&ld,&dummy,&ld,&a,&lwork,&rdummy,&idummy,&info));
418: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
419: #endif
421: #if !PetscDefined(USE_COMPLEX)
422: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
423: alpha = ds->rwork;
424: beta = ds->rwork+ds->ld;
425: PetscCallLAPACKInfo("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,&lwork,ds->iwork,&info));
426: #else
427: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
428: alpha = ds->rwork+2*ds->ld;
429: beta = ds->rwork+3*ds->ld;
430: PetscCallLAPACKInfo("LAPACKggsvd3",LAPACKggsvd3_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,&lwork,ds->rwork,ds->iwork,&info));
431: #endif
433: #else /* defined(SLEPC_MISSING_LAPACK_GGSVD3) */
435: lwork = PetscMax(PetscMax(3*n,m),p)+n;
436: #if !PetscDefined(USE_COMPLEX)
437: PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
438: alpha = ds->rwork;
439: beta = ds->rwork+ds->ld;
440: PetscCallLAPACKInfo("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,ds->iwork,&info));
441: #else
442: PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
443: alpha = ds->rwork+2*ds->ld;
444: beta = ds->rwork+3*ds->ld;
445: PetscCallLAPACKInfo("LAPACKggsvd",LAPACKggsvd_("U","V","Q",&m1,&n1,&p1,&k,&l,A+off,&ld,B+off,&ld,alpha,beta,U+off,&ld,V+off,&ld,X+off,&ld,ds->work,ds->rwork,ds->iwork,&info));
446: #endif
448: #endif
450: PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
451: if (ds->compact) {
452: PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
453: PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
454: /* R is the identity matrix (except the sign) */
455: for (i=lc;i<n;i++) {
456: if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
457: for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
458: }
459: }
460: PetscCall(PetscArrayzero(T+ld,m-1));
461: PetscCall(PetscArrayzero(T+2*ld,n-1));
462: for (i=lc;i<n;i++) {
463: T[i] = alpha[i-lc];
464: D[i] = beta[i-lc];
465: if (D[i]==0.0) wr[i] = PETSC_INFINITY;
466: else wr[i] = T[i]/D[i];
467: }
468: ds->t = n;
469: PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
470: PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
471: } else {
472: /* X = X*inv(R) */
473: q = PetscMin(m,n);
474: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&q,&sone,A,&ld,X,&ld));
475: if (m<n) {
476: r = n-m;
477: PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&r,&m,&sone,X,&ld,A,&ld,&smone,X+m*ld,&ld));
478: PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&r,&sone,B+m*ld,&ld,X+m*ld,&ld));
479: }
480: if (k>0) {
481: for (i=k;i<PetscMin(m,k+l);i++) {
482: PetscCall(PetscArraycpy(X+(i-k)*ld,X+i*ld,ld));
483: PetscCall(PetscArraycpy(U+(i-k)*ld,U+i*ld,ld));
484: }
485: }
486: /* singular values */
487: PetscCall(PetscArrayzero(A,ld*ld));
488: PetscCall(PetscArrayzero(B,ld*ld));
489: for (j=k;j<PetscMin(m,k+l);j++) {
490: A[(j-k)*(1+ld)] = alpha[j];
491: B[(j-k)*(1+ld)] = beta[j];
492: wr[j-k] = alpha[j]/beta[j];
493: }
494: ds->t = PetscMin(m,k+l)-k; /* set number of computed values */
495: }
496: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
497: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
498: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
499: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
500: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: static PetscErrorCode DSCond_GSVD(DS ds,PetscReal *cond)
505: {
506: DS_GSVD *ctx = (DS_GSVD*)ds->data;
507: PetscBLASInt lwork,lrwork=0,m,n,p,ld;
508: PetscScalar *A,*work;
509: const PetscScalar *M;
510: PetscReal *sigma,conda,condb;
511: #if PetscDefined(USE_COMPLEX)
512: PetscReal *rwork;
513: #endif
515: PetscFunctionBegin;
516: PetscCall(PetscBLASIntCast(ds->n,&m));
517: PetscCall(PetscBLASIntCast(ctx->m,&n));
518: PetscCall(PetscBLASIntCast(ctx->p,&p));
519: PetscCall(PetscBLASIntCast(ds->ld,&ld));
520: lwork = 5*n;
521: #if PetscDefined(USE_COMPLEX)
522: lrwork = 5*n;
523: #endif
524: PetscCall(DSAllocateWork_Private(ds,ld*n+lwork,n+lrwork,0));
525: A = ds->work;
526: work = ds->work+ld*n;
527: sigma = ds->rwork;
528: #if PetscDefined(USE_COMPLEX)
529: rwork = ds->rwork+n;
530: #endif
531: if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));
533: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&M));
534: PetscCall(PetscArraycpy(A,M,ld*n));
535: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&M));
536: #if PetscDefined(USE_COMPLEX)
537: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
538: #else
539: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
540: #endif
541: conda = sigma[0]/sigma[PetscMin(m,n)-1];
543: PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&M));
544: PetscCall(PetscArraycpy(A,M,ld*n));
545: PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&M));
546: #if PetscDefined(USE_COMPLEX)
547: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
548: #else
549: PetscCallLAPACKInfo("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
550: #endif
551: condb = sigma[0]/sigma[PetscMin(p,n)-1];
553: *cond = PetscMax(conda,condb);
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: #if !PetscDefined(HAVE_MPIUNI)
558: static PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
559: {
560: DS_GSVD *ctx = (DS_GSVD*)ds->data;
561: PetscInt ld=ds->ld,l=ds->l,k=0,kr=0;
562: PetscMPIInt m,rank,off=0,size,n,ldn,ld3;
563: PetscScalar *A,*U,*V,*X;
564: PetscReal *T;
566: PetscFunctionBegin;
567: PetscCall(PetscMPIIntCast(ctx->m,&m));
568: if (ds->compact) kr = 3*ld;
569: else k = 2*(m-l)*ld;
570: if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
571: if (eigr) k += m-l;
572: PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
573: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
574: PetscCall(PetscMPIIntCast(m-l,&n));
575: PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
576: PetscCall(PetscMPIIntCast(3*ld,&ld3));
577: if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
578: else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
579: if (ds->state>DS_STATE_RAW) {
580: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
581: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
582: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
583: }
584: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
585: if (!rank) {
586: if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
587: else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
588: if (ds->state>DS_STATE_RAW) {
589: PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
590: PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
591: PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
592: }
593: if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
594: }
595: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
596: if (rank) {
597: if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
598: else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
599: if (ds->state>DS_STATE_RAW) {
600: PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
601: PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
602: PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
603: }
604: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
605: }
606: if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
607: else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
608: if (ds->state>DS_STATE_RAW) {
609: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
610: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
611: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
612: }
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
615: #endif
617: static PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
618: {
619: DS_GSVD *ctx = (DS_GSVD*)ds->data;
621: PetscFunctionBegin;
622: PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
623: switch (t) {
624: case DS_MAT_A:
625: *rows = ds->n;
626: *cols = ds->extrarow? ctx->m+1: ctx->m;
627: break;
628: case DS_MAT_B:
629: *rows = ctx->p;
630: *cols = ds->extrarow? ctx->m+1: ctx->m;
631: break;
632: case DS_MAT_T:
633: *rows = ds->n;
634: *cols = PetscDefined(USE_COMPLEX)? 2: 3;
635: break;
636: case DS_MAT_D:
637: *rows = ctx->p;
638: *cols = 1;
639: break;
640: case DS_MAT_U:
641: *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
642: *cols = ds->n;
643: break;
644: case DS_MAT_V:
645: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
646: *cols = ctx->p;
647: break;
648: case DS_MAT_X:
649: *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
650: *cols = ctx->m;
651: break;
652: default:
653: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
654: }
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }
658: static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
659: {
660: DS_GSVD *ctx = (DS_GSVD*)ds->data;
662: PetscFunctionBegin;
663: DSCheckAlloc(ds,1);
664: if (m == PETSC_DETERMINE) {
665: ctx->m = ds->ld;
666: } else if (m != PETSC_CURRENT) {
667: PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
668: ctx->m = m;
669: }
670: if (p == PETSC_DETERMINE) {
671: ctx->p = ds->n;
672: } else if (p != PETSC_CURRENT) {
673: PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
674: ctx->p = p;
675: }
676: PetscFunctionReturn(PETSC_SUCCESS);
677: }
679: /*@
680: DSGSVDSetDimensions - Sets the number of columns and rows for a `DSGSVD`.
682: Logically Collective
684: Input Parameters:
685: + ds - the direct solver context
686: . m - the number of columns
687: - p - the number of rows for the second matrix (`DS_MAT_B`)
689: Notes:
690: This call is complementary to `DSSetDimensions()`, to provide two dimensions
691: that are specific to this `DS` type. The number of rows for the first matrix
692: (`DS_MAT_A`) is set by `DSSetDimensions()`.
694: Use `PETSC_CURRENT` to leave any of the values unchanged. Use `PETSC_DETERMINE`
695: to set `m` to the leading dimension and `p` to the number of columns of $B$.
697: Level: intermediate
699: .seealso: [](sec:ds), `DSGSVD`, `DSGSVDGetDimensions()`, `DSSetDimensions()`
700: @*/
701: PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
702: {
703: PetscFunctionBegin;
707: PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
712: {
713: DS_GSVD *ctx = (DS_GSVD*)ds->data;
715: PetscFunctionBegin;
716: if (m) *m = ctx->m;
717: if (p) *p = ctx->p;
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: DSGSVDGetDimensions - Returns the number of columns and rows for a `DSGSVD`.
724: Not Collective
726: Input Parameter:
727: . ds - the direct solver context
729: Output Parameters:
730: + m - the number of columns
731: - p - the number of rows for the second problem matrix (`DS_MAT_B`)
733: Level: intermediate
735: .seealso: [](sec:ds), `DSGSVD`, `DSGSVDSetDimensions()`
736: @*/
737: PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
738: {
739: PetscFunctionBegin;
741: PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode DSDestroy_GSVD(DS ds)
746: {
747: PetscFunctionBegin;
748: PetscCall(PetscFree(ds->data));
749: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
750: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: static PetscErrorCode DSReallocate_GSVD(DS ds,PetscInt ld)
755: {
756: PetscInt i,*perm=ds->perm;
758: PetscFunctionBegin;
759: for (i=0;i<DS_NUM_MAT;i++) {
760: if (i!=DS_MAT_A && i!=DS_MAT_B && i!=DS_MAT_X && i!=DS_MAT_U && i!=DS_MAT_V && i!=DS_MAT_T && i!=DS_MAT_D) PetscCall(MatDestroy(&ds->omat[i]));
761: }
763: PetscCall(DSReallocateMat_Private(ds,DS_MAT_A,ld));
764: PetscCall(DSReallocateMat_Private(ds,DS_MAT_B,ld));
765: PetscCall(DSReallocateMat_Private(ds,DS_MAT_X,ld));
766: PetscCall(DSReallocateMat_Private(ds,DS_MAT_U,ld));
767: PetscCall(DSReallocateMat_Private(ds,DS_MAT_V,ld));
768: PetscCall(DSReallocateMat_Private(ds,DS_MAT_T,ld));
769: PetscCall(DSReallocateMat_Private(ds,DS_MAT_D,ld));
771: PetscCall(PetscMalloc1(ld,&ds->perm));
772: PetscCall(PetscArraycpy(ds->perm,perm,ds->ld));
773: PetscCall(PetscFree(perm));
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: /*MC
778: DSGSVD - Dense Generalized Singular Value Decomposition.
780: Notes:
781: The problem is expressed as $AX = UC$, $BX = VS$, where $A$ and $B$ are
782: matrices with the same number of columns, $m$, $U$ and $V$ are orthogonal
783: (unitary), and $X$ is an $m\times m$ invertible matrix. The `DS` object does not
784: expose matrices $C$ and $S$, instead the singular values $\sigma_i$, which are
785: the ratios $c_i/s_i$, are returned in the arguments of `DSSolve()`.
786: Note that the number of columns of the returned $X$, $U$, $V$ may be smaller
787: in the case that some $c_i$ or $s_i$ are zero.
789: The number of rows of $A$ (and $U$) is the value $n$ passed with `DSSetDimensions()`.
790: The number of columns $m$ and the number of rows of $B$ (and $V$) must be
791: set via `DSGSVDSetDimensions()`.
793: Internally, LAPACK's representation is used, $U^*AQ = C[0,R]$, $V^*BQ = S[0, R]$,
794: where $X = QR^{-1}$ is computed at the end of `DSSolve()`.
796: If the compact storage format is selected, then a simplified problem is
797: solved, where $A$ and $B$ are bidiagonal (possibly with an arrow), and the
798: stacked matrix $\begin{bmatrix}A^* &B^*\end{bmatrix}^*$
799: is assumed to have orthonormal columns. We consider two cases,
801: 1. $A$ and $B$ are square $m\times m$ upper bidiagonal, or
802: 2. $A$ is lower bidiagonal with $m+1$ rows and $B$ is square upper bidiagonal.
804: In these cases, $R=I$ so it corresponds to the CS decomposition {cite:p}`Gol96`.
805: The first matrix is stored in two diagonals of `DS_MAT_T`, while the second matrix
806: is stored in `DS_MAT_D` and the remaining diagonal of `DS_MAT_T`.
808: Allowed arguments of `DSVectors()` are `DS_MAT_U`, `DS_MAT_V` and `DS_MAT_X`.
810: Used DS matrices:
811: + `DS_MAT_A` - first problem matrix
812: . `DS_MAT_B` - second problem matrix
813: . `DS_MAT_T` - first upper bidiagonal matrix (if compact storage is selected)
814: . `DS_MAT_D` - second upper bidiagonal matrix (if compact storage is selected)
815: . `DS_MAT_U` - (upper) left generalized singular vectors
816: . `DS_MAT_V` - (lower) left generalized singular vectors
817: - `DS_MAT_X` - right generalized singular vectors
819: Implemented methods:
820: . 0 - Lapack (`_ggsvd3` if available, or `_ggsvd`)
822: Level: beginner
824: .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSGSVDSetDimensions()`
825: M*/
826: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
827: {
828: DS_GSVD *ctx;
830: PetscFunctionBegin;
831: PetscCall(PetscNew(&ctx));
832: ds->data = (void*)ctx;
834: ds->ops->allocate = DSAllocate_GSVD;
835: ds->ops->view = DSView_GSVD;
836: ds->ops->vectors = DSVectors_GSVD;
837: ds->ops->sort = DSSort_GSVD;
838: ds->ops->solve[0] = DSSolve_GSVD;
839: #if !PetscDefined(HAVE_MPIUNI)
840: ds->ops->synchronize = DSSynchronize_GSVD;
841: #endif
842: ds->ops->truncate = DSTruncate_GSVD;
843: ds->ops->update = DSUpdateExtraRow_GSVD;
844: ds->ops->cond = DSCond_GSVD;
845: ds->ops->matgetsize = DSMatGetSize_GSVD;
846: ds->ops->destroy = DSDestroy_GSVD;
847: ds->ops->reallocate = DSReallocate_GSVD;
848: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
849: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }