Actual source code: dspep.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: #include <slepc/private/dsimpl.h>       /*I "slepcds.h" I*/
 12: #include <slepcblaslapack.h>

 14: typedef struct {
 15:   PetscInt  d;              /* polynomial degree */
 16:   PetscReal *pbc;           /* polynomial basis coefficients */
 17: } DS_PEP;

 19: PetscErrorCode DSAllocate_PEP(DS ds,PetscInt ld)
 20: {
 22:   DS_PEP         *ctx = (DS_PEP*)ds->data;
 23:   PetscInt       i;

 26:   if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
 27:   DSAllocateMat_Private(ds,DS_MAT_X);
 28:   DSAllocateMat_Private(ds,DS_MAT_Y);
 29:   for (i=0;i<=ctx->d;i++) {
 30:     DSAllocateMat_Private(ds,DSMatExtra[i]);
 31:   }
 32:   PetscFree(ds->perm);
 33:   PetscMalloc1(ld*ctx->d,&ds->perm);
 34:   PetscLogObjectMemory((PetscObject)ds,ld*ctx->d*sizeof(PetscInt));
 35:   return(0);
 36: }

 38: PetscErrorCode DSView_PEP(DS ds,PetscViewer viewer)
 39: {
 40:   PetscErrorCode    ierr;
 41:   DS_PEP            *ctx = (DS_PEP*)ds->data;
 42:   PetscViewerFormat format;
 43:   PetscInt          i;

 46:   PetscViewerGetFormat(viewer,&format);
 47:   PetscViewerASCIIPrintf(viewer,"polynomial degree: %D\n",ctx->d);
 48:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
 49:   for (i=0;i<=ctx->d;i++) {
 50:     DSViewMat(ds,viewer,DSMatExtra[i]);
 51:   }
 52:   if (ds->state>DS_STATE_INTERMEDIATE) {
 53:     DSViewMat(ds,viewer,DS_MAT_X);
 54:   }
 55:   return(0);
 56: }

 58: PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
 59: {
 61:   if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
 62:   switch (mat) {
 63:     case DS_MAT_X:
 64:       break;
 65:     case DS_MAT_Y:
 66:       break;
 67:     default:
 68:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
 69:   }
 70:   return(0);
 71: }

 73: PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
 74: {
 76:   DS_PEP         *ctx = (DS_PEP*)ds->data;
 77:   PetscInt       n,i,j,k,p,*perm,told,ld;
 78:   PetscScalar    *A,*X,*Y,rtmp,rtmp2;

 81:   if (!ds->sc) return(0);
 82:   n = ds->n*ctx->d;
 83:   A  = ds->mat[DS_MAT_A];
 84:   perm = ds->perm;
 85:   for (i=0;i<n;i++) perm[i] = i;
 86:   told = ds->t;
 87:   ds->t = n;  /* force the sorting routines to consider d*n eigenvalues */
 88:   if (rr) {
 89:     DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
 90:   } else {
 91:     DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE);
 92:   }
 93:   ds->t = told;  /* restore value of t */
 94:   for (i=0;i<n;i++) A[i]  = wr[perm[i]];
 95:   for (i=0;i<n;i++) wr[i] = A[i];
 96:   for (i=0;i<n;i++) A[i]  = wi[perm[i]];
 97:   for (i=0;i<n;i++) wi[i] = A[i];
 98:   /* cannot use DSPermuteColumns_Private() since matrix is not square */
 99:   ld = ds->ld;
100:   X  = ds->mat[DS_MAT_X];
101:   Y  = ds->mat[DS_MAT_Y];
102:   for (i=0;i<n;i++) {
103:     p = perm[i];
104:     if (p != i) {
105:       j = i + 1;
106:       while (perm[j] != i) j++;
107:       perm[j] = p; perm[i] = i;
108:       /* swap columns i and j */
109:       for (k=0;k<ds->n;k++) {
110:         rtmp  = X[k+p*ld]; X[k+p*ld] = X[k+i*ld]; X[k+i*ld] = rtmp;
111:         rtmp2 = Y[k+p*ld]; Y[k+p*ld] = Y[k+i*ld]; Y[k+i*ld] = rtmp2;
112:       }
113:     }
114:   }
115:   return(0);
116: }

118: PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
119: {
120: #if defined(SLEPC_MISSING_LAPACK_GGEV)
122:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GGEV - Lapack routine is unavailable");
123: #else
125:   DS_PEP         *ctx = (DS_PEP*)ds->data;
126:   PetscInt       i,j,k,off;
127:   PetscScalar    *A,*B,*W,*X,*U,*Y,*E,*work,*beta,norm;
128:   PetscReal      *ca,*cb,*cg;
129:   PetscBLASInt   info,n,ldd,nd,lrwork=0,lwork,one=1;
130: #if defined(PETSC_USE_COMPLEX)
131:   PetscReal      *rwork;
132: #else
133:   PetscScalar    norm0;
134: #endif

137:   if (!ds->mat[DS_MAT_A]) {
138:     DSAllocateMat_Private(ds,DS_MAT_A);
139:   }
140:   if (!ds->mat[DS_MAT_B]) {
141:     DSAllocateMat_Private(ds,DS_MAT_B);
142:   }
143:   if (!ds->mat[DS_MAT_W]) {
144:     DSAllocateMat_Private(ds,DS_MAT_W);
145:   }
146:   if (!ds->mat[DS_MAT_U]) {
147:     DSAllocateMat_Private(ds,DS_MAT_U);
148:   }
149:   PetscBLASIntCast(ds->n*ctx->d,&nd);
150:   PetscBLASIntCast(ds->n,&n);
151:   PetscBLASIntCast(ds->ld*ctx->d,&ldd);
152: #if defined(PETSC_USE_COMPLEX)
153:   PetscBLASIntCast(nd+2*nd,&lwork);
154:   PetscBLASIntCast(8*nd,&lrwork);
155: #else
156:   PetscBLASIntCast(nd+8*nd,&lwork);
157: #endif
158:   DSAllocateWork_Private(ds,lwork,lrwork,0);
159:   beta = ds->work;
160:   work = ds->work + nd;
161:   lwork -= nd;
162:   A = ds->mat[DS_MAT_A];
163:   B = ds->mat[DS_MAT_B];
164:   W = ds->mat[DS_MAT_W];
165:   U = ds->mat[DS_MAT_U];
166:   X = ds->mat[DS_MAT_X];
167:   Y = ds->mat[DS_MAT_Y];
168:   E = ds->mat[DSMatExtra[ctx->d]];

170:   /* build matrices A and B of the linearization */
171:   PetscMemzero(A,ldd*ldd*sizeof(PetscScalar));
172:   if (!ctx->pbc) { /* monomial basis */
173:     for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
174:     for (i=0;i<ctx->d;i++) {
175:       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
176:       for (j=0;j<ds->n;j++) {
177:         PetscMemcpy(A+off+j*ldd,ds->mat[DSMatExtra[i]]+j*ds->ld,ds->n*sizeof(PetscScalar));
178:       }
179:     }
180:   } else {
181:     ca = ctx->pbc;
182:     cb = ca+ctx->d+1;
183:     cg = cb+ctx->d+1;
184:     for (i=0;i<ds->n;i++) {
185:       A[i+(i+ds->n)*ldd] = ca[0];
186:       A[i+i*ldd] = cb[0];
187:     }
188:     for (;i<nd-ds->n;i++) {
189:       j = i/ds->n;
190:       A[i+(i+ds->n)*ldd] = ca[j];
191:       A[i+i*ldd] = cb[j];
192:       A[i+(i-ds->n)*ldd] = cg[j];
193:     }
194:     for (i=0;i<ctx->d-2;i++) {
195:       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
196:       for (j=0;j<ds->n;j++)
197:         for (k=0;k<ds->n;k++)
198:           *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1];
199:     }
200:     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
201:     for (j=0;j<ds->n;j++) 
202:       for (k=0;k<ds->n;k++)
203:         *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1]-E[j*ds->ld+k]*cg[ctx->d-1];
204:     off = (++i)*ds->n*ldd+(ctx->d-1)*ds->n;
205:     for (j=0;j<ds->n;j++) 
206:       for (k=0;k<ds->n;k++)
207:         *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1]-E[j*ds->ld+k]*cb[ctx->d-1];
208:   }
209:   PetscMemzero(B,ldd*ldd*sizeof(PetscScalar));
210:   for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
211:   off = (ctx->d-1)*ds->n*(ldd+1);
212:   for (j=0;j<ds->n;j++) {
213:     for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -E[i+j*ds->ld];
214:   }

216:   /* solve generalized eigenproblem */
217: #if defined(PETSC_USE_COMPLEX)
218:   rwork = ds->rwork;
219:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,rwork,&info));
220: #else
221:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
222: #endif
223:   SlepcCheckLapackInfo("ggev",info);

225:   /* copy eigenvalues */
226:   for (i=0;i<nd;i++) {
227:     if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
228:     else wr[i] /= beta[i];
229: #if !defined(PETSC_USE_COMPLEX)
230:     if (beta[i]==0.0) wi[i] = 0.0;
231:     else wi[i] /= beta[i];
232: #else
233:     if (wi) wi[i] = 0.0;
234: #endif
235:   }

237:   /* copy and normalize eigenvectors */
238:   for (j=0;j<nd;j++) {
239:     PetscMemcpy(X+j*ds->ld,W+j*ldd,ds->n*sizeof(PetscScalar));
240:     PetscMemcpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n*sizeof(PetscScalar));
241:   }
242:   for (j=0;j<nd;j++) {
243: #if !defined(PETSC_USE_COMPLEX)
244:     if (wi[j] != 0.0) {
245:       norm = BLASnrm2_(&n,X+j*ds->ld,&one);
246:       norm0 = BLASnrm2_(&n,X+(j+1)*ds->ld,&one);
247:       norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
248:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
249:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+(j+1)*ds->ld,&one));
250:       norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
251:       norm0 = BLASnrm2_(&n,Y+(j+1)*ds->ld,&one);
252:       norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
253:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
254:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+(j+1)*ds->ld,&one));
255:       j++;
256:     } else
257: #endif
258:     {
259:       norm = 1.0/BLASnrm2_(&n,X+j*ds->ld,&one);
260:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
261:       norm = 1.0/BLASnrm2_(&n,Y+j*ds->ld,&one);
262:       PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
263:     }
264:   }
265:   return(0);
266: #endif
267: }

269: PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
270: {
272:   DS_PEP         *ctx = (DS_PEP*)ds->data;  
273:   PetscInt       ld=ds->ld,k=0;
274:   PetscMPIInt    ldnd,rank,off=0,size,dn;

277:   if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
278:   if (eigr) k += ctx->d*ds->n; 
279:   if (eigi) k += ctx->d*ds->n; 
280:   DSAllocateWork_Private(ds,k,0,0);
281:   PetscMPIIntCast(k*sizeof(PetscScalar),&size);
282:   PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd);
283:   PetscMPIIntCast(ctx->d*ds->n,&dn);
284:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
285:   if (!rank) {
286:     if (ds->state>=DS_STATE_CONDENSED) {
287:       MPI_Pack(ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
288:       MPI_Pack(ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
289:     }
290:     if (eigr) {
291:       MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
292:     }
293:     if (eigi) {
294:       MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
295:     }
296:   }
297:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
298:   if (rank) {
299:     if (ds->state>=DS_STATE_CONDENSED) {
300:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
301:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
302:     }
303:     if (eigr) {
304:       MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
305:     }
306:     if (eigi) {
307:       MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
308:     }
309:   }
310:   return(0);
311: }

313: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
314: {
315:   DS_PEP *ctx = (DS_PEP*)ds->data;

318:   if (d<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
319:   if (d>=DS_NUM_EXTRA) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Only implemented for polynomials of degree at most %D",DS_NUM_EXTRA-1);
320:   ctx->d = d;
321:   return(0);
322: }

324: /*@
325:    DSPEPSetDegree - Sets the polynomial degree for a DSPEP.

327:    Logically Collective on DS

329:    Input Parameters:
330: +  ds - the direct solver context
331: -  d  - the degree

333:    Level: intermediate

335: .seealso: DSPEPGetDegree()
336: @*/
337: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
338: {

344:   PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
345:   return(0);
346: }

348: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
349: {
350:   DS_PEP *ctx = (DS_PEP*)ds->data;

353:   *d = ctx->d;
354:   return(0);
355: }

357: /*@
358:    DSPEPGetDegree - Returns the polynomial degree for a DSPEP.

360:    Not collective

362:    Input Parameter:
363: .  ds - the direct solver context

365:    Output Parameters:
366: .  d - the degree

368:    Level: intermediate

370: .seealso: DSPEPSetDegree()
371: @*/
372: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
373: {

379:   PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
380:   return(0);
381: }

383: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
384: {
386:   DS_PEP         *ctx = (DS_PEP*)ds->data;
387:   PetscInt       i;

390:   if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
391:   if (ctx->pbc) { PetscFree(ctx->pbc); }
392:   PetscMalloc1(3*(ctx->d+1),&ctx->pbc);
393:   for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
394:   ds->state = DS_STATE_RAW;
395:   return(0);
396: }

398: /*@C
399:    DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.

401:    Logically Collective on DS

403:    Input Parameters:
404: +  ds  - the direct solver context
405: -  pbc - the polynomial basis coefficients

407:    Notes:
408:    This function is required only in the case of a polynomial specified in a
409:    non-monomial basis, to provide the coefficients that will be used
410:    during the linearization, multiplying the identity blocks on the three main
411:    diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...) 
412:    the coefficients must be different.

414:    There must be a total of 3*(d+1) coefficients, where d is the degree of the
415:    polynomial. The coefficients are arranged in three groups: alpha, beta, and
416:    gamma, according to the definition of the three-term recurrence. In the case
417:    of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
418:    necessary to invoke this function.

420:    Level: advanced

422: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
423: @*/
424: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
425: {

430:   PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
431:   return(0);
432: }

434: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
435: {
437:   DS_PEP         *ctx = (DS_PEP*)ds->data;
438:   PetscInt       i;

441:   if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
442:   PetscCalloc1(3*(ctx->d+1),pbc);
443:   if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
444:   else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
445:   return(0);
446: }

448: /*@C
449:    DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.

451:    Not collective

453:    Input Parameter:
454: .  ds - the direct solver context

456:    Output Parameters:
457: .  pbc - the polynomial basis coefficients

459:    Note:
460:    The returned array has length 3*(d+1) and should be freed by the user.

462:    Fortran Note:
463:    The calling sequence from Fortran is
464: .vb
465:    DSPEPGetCoefficients(eps,pbc,ierr)
466:    double precision pbc(d+1) output
467: .ve

469:    Level: advanced

471: .seealso: DSPEPSetCoefficients()
472: @*/
473: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
474: {

480:   PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
481:   return(0);
482: }

484: PetscErrorCode DSDestroy_PEP(DS ds)
485: {
487:   DS_PEP         *ctx = (DS_PEP*)ds->data;

490:   if (ctx->pbc) { PetscFree(ctx->pbc); }
491:   PetscFree(ds->data);
492:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL);
493:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL);
494:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL);
495:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL);
496:   return(0);
497: }

499: PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
500: {
501:   DS_PEP *ctx = (DS_PEP*)ds->data;

504:   if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
505:   *rows = ds->n;
506:   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
507:   *cols = ds->n;
508:   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U || t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->d;
509:   return(0);
510: }

512: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
513: {
514:   DS_PEP         *ctx;

518:   PetscNewLog(ds,&ctx);
519:   ds->data = (void*)ctx;

521:   ds->ops->allocate      = DSAllocate_PEP;
522:   ds->ops->view          = DSView_PEP;
523:   ds->ops->vectors       = DSVectors_PEP;
524:   ds->ops->solve[0]      = DSSolve_PEP_QZ;
525:   ds->ops->sort          = DSSort_PEP;
526:   ds->ops->synchronize   = DSSynchronize_PEP;
527:   ds->ops->destroy       = DSDestroy_PEP;
528:   ds->ops->matgetsize    = DSMatGetSize_PEP;
529:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP);
530:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP);
531:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP);
532:   PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP);
533:   return(0);
534: }