Actual source code: dsgsvd.c

slepc-main 2024-11-09
Report Typos and Errors
  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,info,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 !defined (PETSC_USE_COMPLEX)
414:   PetscCallBLAS("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:   PetscCallBLAS("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 !defined (PETSC_USE_COMPLEX)
422:   PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
423:   alpha = ds->rwork;
424:   beta  = ds->rwork+ds->ld;
425:   PetscCallBLAS("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:   PetscCallBLAS("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
432:   SlepcCheckLapackInfo("ggsvd3",info);

434: #else  /* defined(SLEPC_MISSING_LAPACK_GGSVD3) */

436:   lwork = PetscMax(PetscMax(3*n,m),p)+n;
437: #if !defined (PETSC_USE_COMPLEX)
438:   PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
439:   alpha = ds->rwork;
440:   beta  = ds->rwork+ds->ld;
441:   PetscCallBLAS("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));
442: #else
443:   PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
444:   alpha = ds->rwork+2*ds->ld;
445:   beta  = ds->rwork+3*ds->ld;
446:   PetscCallBLAS("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));
447: #endif
448:   SlepcCheckLapackInfo("ggsvd",info);

450: #endif

452:   PetscCheck(k+l>=n1,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"The rank deficient case not supported yet");
453:   if (ds->compact) {
454:     PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
455:     PetscCall(DSGetArrayReal(ds,DS_MAT_D,&D));
456:     /* R is the identity matrix (except the sign) */
457:     for (i=lc;i<n;i++) {
458:       if (PetscRealPart(A[i+i*ld])<0.0) { /* scale column i */
459:         for (j=lc;j<n;j++) X[j+i*ld] = -X[j+i*ld];
460:       }
461:     }
462:     PetscCall(PetscArrayzero(T+ld,m-1));
463:     PetscCall(PetscArrayzero(T+2*ld,n-1));
464:     for (i=lc;i<n;i++) {
465:       T[i] = alpha[i-lc];
466:       D[i] = beta[i-lc];
467:       if (D[i]==0.0) wr[i] = PETSC_INFINITY;
468:       else wr[i] = T[i]/D[i];
469:     }
470:     ds->t = n;
471:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_D,&D));
472:     PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
473:   } else {
474:     /* X = X*inv(R) */
475:     q = PetscMin(m,n);
476:     PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&q,&sone,A,&ld,X,&ld));
477:     if (m<n) {
478:       r = n-m;
479:       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&r,&m,&sone,X,&ld,A,&ld,&smone,X+m*ld,&ld));
480:       PetscCallBLAS("BLAStrsm",BLAStrsm_("R","U","N","N",&n,&r,&sone,B+m*ld,&ld,X+m*ld,&ld));
481:     }
482:     if (k>0) {
483:       for (i=k;i<PetscMin(m,k+l);i++) {
484:         PetscCall(PetscArraycpy(X+(i-k)*ld,X+i*ld,ld));
485:         PetscCall(PetscArraycpy(U+(i-k)*ld,U+i*ld,ld));
486:       }
487:     }
488:     /* singular values */
489:     PetscCall(PetscArrayzero(A,ld*ld));
490:     PetscCall(PetscArrayzero(B,ld*ld));
491:     for (j=k;j<PetscMin(m,k+l);j++) {
492:       A[(j-k)*(1+ld)] = alpha[j];
493:       B[(j-k)*(1+ld)] = beta[j];
494:       wr[j-k] = alpha[j]/beta[j];
495:     }
496:     ds->t = PetscMin(m,k+l)-k; /* set number of computed values */
497:   }
498:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
499:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
500:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
501:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
502:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
503:   PetscFunctionReturn(PETSC_SUCCESS);
504: }

506: static PetscErrorCode DSCond_GSVD(DS ds,PetscReal *cond)
507: {
508:   DS_GSVD           *ctx = (DS_GSVD*)ds->data;
509:   PetscBLASInt      lwork,lrwork=0,info,m,n,p,ld;
510:   PetscScalar       *A,*work;
511:   const PetscScalar *M;
512:   PetscReal         *sigma,conda,condb;
513: #if defined(PETSC_USE_COMPLEX)
514:   PetscReal         *rwork;
515: #endif

517:   PetscFunctionBegin;
518:   PetscCall(PetscBLASIntCast(ds->n,&m));
519:   PetscCall(PetscBLASIntCast(ctx->m,&n));
520:   PetscCall(PetscBLASIntCast(ctx->p,&p));
521:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
522:   lwork = 5*n;
523: #if defined(PETSC_USE_COMPLEX)
524:   lrwork = 5*n;
525: #endif
526:   PetscCall(DSAllocateWork_Private(ds,ld*n+lwork,n+lrwork,0));
527:   A     = ds->work;
528:   work  = ds->work+ld*n;
529:   sigma = ds->rwork;
530: #if defined(PETSC_USE_COMPLEX)
531:   rwork = ds->rwork+n;
532: #endif
533:   if (ds->compact) PetscCall(DSSwitchFormat_GSVD(ds));

535:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_A],&M));
536:   PetscCall(PetscArraycpy(A,M,ld*n));
537:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_A],&M));
538: #if defined(PETSC_USE_COMPLEX)
539:   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
540: #else
541:   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&m,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
542: #endif
543:   SlepcCheckLapackInfo("gesvd",info);
544:   conda = sigma[0]/sigma[PetscMin(m,n)-1];

546:   PetscCall(MatDenseGetArrayRead(ds->omat[DS_MAT_B],&M));
547:   PetscCall(PetscArraycpy(A,M,ld*n));
548:   PetscCall(MatDenseRestoreArrayRead(ds->omat[DS_MAT_B],&M));
549: #if defined(PETSC_USE_COMPLEX)
550:   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,rwork,&info));
551: #else
552:   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("N","N",&p,&n,A,&ld,sigma,NULL,&ld,NULL,&ld,work,&lwork,&info));
553: #endif
554:   SlepcCheckLapackInfo("gesvd",info);
555:   condb = sigma[0]/sigma[PetscMin(p,n)-1];

557:   *cond = PetscMax(conda,condb);
558:   PetscFunctionReturn(PETSC_SUCCESS);
559: }

561: #if !defined(PETSC_HAVE_MPIUNI)
562: static PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
563: {
564:   DS_GSVD        *ctx = (DS_GSVD*)ds->data;
565:   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
566:   PetscMPIInt    m,rank,off=0,size,n,ldn,ld3;
567:   PetscScalar    *A,*U,*V,*X;
568:   PetscReal      *T;

570:   PetscFunctionBegin;
571:   PetscCall(PetscMPIIntCast(ctx->m,&m));
572:   if (ds->compact) kr = 3*ld;
573:   else k = 2*(m-l)*ld;
574:   if (ds->state>DS_STATE_RAW) k += 3*(m-l)*ld;
575:   if (eigr) k += m-l;
576:   PetscCall(DSAllocateWork_Private(ds,k+kr,0,0));
577:   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar)+kr*sizeof(PetscReal),&size));
578:   PetscCall(PetscMPIIntCast(m-l,&n));
579:   PetscCall(PetscMPIIntCast(ld*(m-l),&ldn));
580:   PetscCall(PetscMPIIntCast(3*ld,&ld3));
581:   if (ds->compact) PetscCall(DSGetArrayReal(ds,DS_MAT_T,&T));
582:   else PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
583:   if (ds->state>DS_STATE_RAW) {
584:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
585:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
586:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
587:   }
588:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
589:   if (!rank) {
590:     if (ds->compact) PetscCallMPI(MPI_Pack(T,ld3,MPIU_REAL,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
591:     else PetscCallMPI(MPI_Pack(A+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
592:     if (ds->state>DS_STATE_RAW) {
593:       PetscCallMPI(MPI_Pack(U+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
594:       PetscCallMPI(MPI_Pack(V+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
595:       PetscCallMPI(MPI_Pack(X+l*ld,ldn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
596:     }
597:     if (eigr) PetscCallMPI(MPI_Pack(eigr+l,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
598:   }
599:   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
600:   if (rank) {
601:     if (ds->compact) PetscCallMPI(MPI_Unpack(ds->work,size,&off,T,ld3,MPIU_REAL,PetscObjectComm((PetscObject)ds)));
602:     else PetscCallMPI(MPI_Unpack(ds->work,size,&off,A+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
603:     if (ds->state>DS_STATE_RAW) {
604:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,U+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
605:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,V+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
606:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,X+l*ld,ldn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
607:     }
608:     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr+l,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
609:   }
610:   if (ds->compact) PetscCall(DSRestoreArrayReal(ds,DS_MAT_T,&T));
611:   else PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
612:   if (ds->state>DS_STATE_RAW) {
613:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
614:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));
615:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
616:   }
617:   PetscFunctionReturn(PETSC_SUCCESS);
618: }
619: #endif

621: static PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
622: {
623:   DS_GSVD *ctx = (DS_GSVD*)ds->data;

625:   PetscFunctionBegin;
626:   PetscCheck(ctx->m,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"You should set the other dimensions with DSGSVDSetDimensions()");
627:   switch (t) {
628:     case DS_MAT_A:
629:       *rows = ds->n;
630:       *cols = ds->extrarow? ctx->m+1: ctx->m;
631:       break;
632:     case DS_MAT_B:
633:       *rows = ctx->p;
634:       *cols = ds->extrarow? ctx->m+1: ctx->m;
635:       break;
636:     case DS_MAT_T:
637:       *rows = ds->n;
638:       *cols = PetscDefined(USE_COMPLEX)? 2: 3;
639:       break;
640:     case DS_MAT_D:
641:       *rows = ctx->p;
642:       *cols = 1;
643:       break;
644:     case DS_MAT_U:
645:       *rows = ds->state==DS_STATE_TRUNCATED? ds->t: ds->n;
646:       *cols = ds->n;
647:       break;
648:     case DS_MAT_V:
649:       *rows = ds->state==DS_STATE_TRUNCATED? ctx->tp: ctx->p;
650:       *cols = ctx->p;
651:       break;
652:     case DS_MAT_X:
653:       *rows = ds->state==DS_STATE_TRUNCATED? ctx->tm: ctx->m;
654:       *cols = ctx->m;
655:       break;
656:     default:
657:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid t parameter");
658:   }
659:   PetscFunctionReturn(PETSC_SUCCESS);
660: }

662: static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
663: {
664:   DS_GSVD *ctx = (DS_GSVD*)ds->data;

666:   PetscFunctionBegin;
667:   DSCheckAlloc(ds,1);
668:   if (m == PETSC_DETERMINE) {
669:     ctx->m = ds->ld;
670:   } else if (m != PETSC_CURRENT) {
671:     PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
672:     ctx->m = m;
673:   }
674:   if (p == PETSC_DETERMINE) {
675:     ctx->p = ds->n;
676:   } else if (p != PETSC_CURRENT) {
677:     PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
678:     ctx->p = p;
679:   }
680:   PetscFunctionReturn(PETSC_SUCCESS);
681: }

683: /*@
684:    DSGSVDSetDimensions - Sets the number of columns and rows for a DSGSVD.

686:    Logically Collective

688:    Input Parameters:
689: +  ds - the direct solver context
690: .  m  - the number of columns
691: -  p  - the number of rows for the second matrix (B)

693:    Notes:
694:    This call is complementary to DSSetDimensions(), to provide two dimensions
695:    that are specific to this DS type. The number of rows for the first matrix (A)
696:    is set by DSSetDimensions().

698:    Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
699:    to set m to the leading dimension and p to the number of columns of B.

701:    Level: intermediate

703: .seealso: DSGSVDGetDimensions(), DSSetDimensions()
704: @*/
705: PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
706: {
707:   PetscFunctionBegin;
711:   PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
712:   PetscFunctionReturn(PETSC_SUCCESS);
713: }

715: static PetscErrorCode DSGSVDGetDimensions_GSVD(DS ds,PetscInt *m,PetscInt *p)
716: {
717:   DS_GSVD *ctx = (DS_GSVD*)ds->data;

719:   PetscFunctionBegin;
720:   if (m) *m = ctx->m;
721:   if (p) *p = ctx->p;
722:   PetscFunctionReturn(PETSC_SUCCESS);
723: }

725: /*@
726:    DSGSVDGetDimensions - Returns the number of columns and rows for a DSGSVD.

728:    Not Collective

730:    Input Parameter:
731: .  ds - the direct solver context

733:    Output Parameters:
734: +  m - the number of columns
735: -  p - the number of rows for the second matrix (B)

737:    Level: intermediate

739: .seealso: DSGSVDSetDimensions()
740: @*/
741: PetscErrorCode DSGSVDGetDimensions(DS ds,PetscInt *m,PetscInt *p)
742: {
743:   PetscFunctionBegin;
745:   PetscUseMethod(ds,"DSGSVDGetDimensions_C",(DS,PetscInt*,PetscInt*),(ds,m,p));
746:   PetscFunctionReturn(PETSC_SUCCESS);
747: }

749: static PetscErrorCode DSDestroy_GSVD(DS ds)
750: {
751:   PetscFunctionBegin;
752:   PetscCall(PetscFree(ds->data));
753:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",NULL));
754:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",NULL));
755:   PetscFunctionReturn(PETSC_SUCCESS);
756: }

758: /*MC
759:    DSGSVD - Dense Generalized Singular Value Decomposition.

761:    Level: beginner

763:    Notes:
764:    The problem is expressed as A*X = U*C, B*X = V*S, where A and B are
765:    matrices with the same number of columns, m, U and V are orthogonal
766:    (unitary), and X is an mxm invertible matrix. The DS object does not
767:    expose matrices C and S, instead the singular values sigma, which are
768:    the ratios c_i/s_i, are returned in the arguments of DSSolve().
769:    Note that the number of columns of the returned X, U, V may be smaller
770:    in the case that some c_i or s_i are zero.

772:    The number of rows of A (and U) is the value n passed with DSSetDimensions().
773:    The number of columns m and the number of rows of B (and V) must be
774:    set via DSGSVDSetDimensions().

776:    Internally, LAPACK's representation is used, U'*A*Q = C*[0 R], V'*B*Q = S*[0 R],
777:    where X = Q*inv(R) is computed at the end of DSSolve().

779:    If the compact storage format is selected, then a simplified problem is
780:    solved, where A and B are bidiagonal (possibly with an arrow), and [A;B]
781:    is assumed to have orthonormal columns. We consider two cases: (1) A and B
782:    are square mxm upper bidiagonal, and (2) A is lower bidiagonal with m+1
783:    rows and B is square upper bidiagonal. In these cases, R=I so it
784:    corresponds to the CS decomposition. The first matrix is stored in two
785:    diagonals of DS_MAT_T, while the second matrix is stored in DS_MAT_D
786:    and the remaining diagonal of DS_MAT_T.

788:    Allowed arguments of DSVectors() are DS_MAT_U, DS_MAT_V and DS_MAT_X.

790:    Used DS matrices:
791: +  DS_MAT_A - first problem matrix
792: .  DS_MAT_B - second problem matrix
793: .  DS_MAT_T - first upper bidiagonal matrix (if compact storage is selected)
794: .  DS_MAT_D - second upper bidiagonal matrix (if compact storage is selected)
795: .  DS_MAT_U - (upper) left generalized singular vectors
796: .  DS_MAT_V - (lower) left generalized singular vectors
797: -  DS_MAT_X - right generalized singular vectors

799:    Implemented methods:
800: .  0 - Lapack (_ggsvd3 if available, or _ggsvd)

802: .seealso: DSCreate(), DSSetType(), DSType, DSGSVDSetDimensions()
803: M*/
804: SLEPC_EXTERN PetscErrorCode DSCreate_GSVD(DS ds)
805: {
806:   DS_GSVD        *ctx;

808:   PetscFunctionBegin;
809:   PetscCall(PetscNew(&ctx));
810:   ds->data = (void*)ctx;

812:   ds->ops->allocate      = DSAllocate_GSVD;
813:   ds->ops->view          = DSView_GSVD;
814:   ds->ops->vectors       = DSVectors_GSVD;
815:   ds->ops->sort          = DSSort_GSVD;
816:   ds->ops->solve[0]      = DSSolve_GSVD;
817: #if !defined(PETSC_HAVE_MPIUNI)
818:   ds->ops->synchronize   = DSSynchronize_GSVD;
819: #endif
820:   ds->ops->truncate      = DSTruncate_GSVD;
821:   ds->ops->update        = DSUpdateExtraRow_GSVD;
822:   ds->ops->cond          = DSCond_GSVD;
823:   ds->ops->matgetsize    = DSMatGetSize_GSVD;
824:   ds->ops->destroy       = DSDestroy_GSVD;
825:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDSetDimensions_C",DSGSVDSetDimensions_GSVD));
826:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSGSVDGetDimensions_C",DSGSVDGetDimensions_GSVD));
827:   PetscFunctionReturn(PETSC_SUCCESS);
828: }