Actual source code: dsgsvd.c

slepc-3.19.2 2023-09-05
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: 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: 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: 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: 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: 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: 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: 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:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
412:   /* workspace query and memory allocation */
413:   lwork = -1;
414: #if !defined (PETSC_USE_COMPLEX)
415:   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));
416:   PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
417: #else
418:   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));
419:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
420: #endif

422: #if !defined (PETSC_USE_COMPLEX)
423:   PetscCall(DSAllocateWork_Private(ds,lwork,2*ds->ld,ds->ld));
424:   alpha = ds->rwork;
425:   beta  = ds->rwork+ds->ld;
426:   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));
427: #else
428:   PetscCall(DSAllocateWork_Private(ds,lwork,4*ds->ld,ds->ld));
429:   alpha = ds->rwork+2*ds->ld;
430:   beta  = ds->rwork+3*ds->ld;
431:   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));
432: #endif
433:   PetscCall(PetscFPTrapPop());
434:   SlepcCheckLapackInfo("ggsvd3",info);

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

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

454: #endif

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

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

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

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

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

563:   *cond = PetscMax(conda,condb);
564:   PetscFunctionReturn(PETSC_SUCCESS);
565: }

567: #if !defined(PETSC_HAVE_MPIUNI)
568: PetscErrorCode DSSynchronize_GSVD(DS ds,PetscScalar eigr[],PetscScalar eigi[])
569: {
570:   DS_GSVD        *ctx = (DS_GSVD*)ds->data;
571:   PetscInt       ld=ds->ld,l=ds->l,k=0,kr=0;
572:   PetscMPIInt    m=ctx->m,rank,off=0,size,n,ldn,ld3;
573:   PetscScalar    *A,*U,*V,*X;
574:   PetscReal      *T;

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

626: PetscErrorCode DSMatGetSize_GSVD(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
627: {
628:   DS_GSVD *ctx = (DS_GSVD*)ds->data;

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

667: static PetscErrorCode DSGSVDSetDimensions_GSVD(DS ds,PetscInt m,PetscInt p)
668: {
669:   DS_GSVD *ctx = (DS_GSVD*)ds->data;

671:   PetscFunctionBegin;
672:   DSCheckAlloc(ds,1);
673:   if (m==PETSC_DECIDE || m==PETSC_DEFAULT) {
674:     ctx->m = ds->ld;
675:   } else {
676:     PetscCheck(m>0 && m<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of m. Must be between 1 and ld");
677:     ctx->m = m;
678:   }
679:   if (p==PETSC_DECIDE || p==PETSC_DEFAULT) {
680:     ctx->p = ds->n;
681:   } else {
682:     PetscCheck(p>0 && p<=ds->ld,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of p. Must be between 1 and ld");
683:     ctx->p = p;
684:   }
685:   PetscFunctionReturn(PETSC_SUCCESS);
686: }

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

691:    Logically Collective

693:    Input Parameters:
694: +  ds - the direct solver context
695: .  m  - the number of columns
696: -  p  - the number of rows for the second matrix (B)

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

703:    Level: intermediate

705: .seealso: DSGSVDGetDimensions(), DSSetDimensions()
706: @*/
707: PetscErrorCode DSGSVDSetDimensions(DS ds,PetscInt m,PetscInt p)
708: {
709:   PetscFunctionBegin;
713:   PetscTryMethod(ds,"DSGSVDSetDimensions_C",(DS,PetscInt,PetscInt),(ds,m,p));
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

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

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

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

730:    Not Collective

732:    Input Parameter:
733: .  ds - the direct solver context

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

739:    Level: intermediate

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

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

760: /*MC
761:    DSGSVD - Dense Generalized Singular Value Decomposition.

763:    Level: beginner

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

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

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

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

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

792:    Used DS matrices:
793: +  DS_MAT_A - first problem matrix
794: .  DS_MAT_B - second problem matrix
795: .  DS_MAT_T - first upper bidiagonal matrix (if compact storage is selected)
796: -  DS_MAT_D - second upper bidiagonal matrix (if compact storage is selected)

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

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

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

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