Actual source code: ex40.c

slepc-3.17.2 2022-08-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: */

 11: static char help[] = "Checking the definite property in quadratic symmetric eigenproblem.\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n> ... dimension of the matrices.\n"
 14:   "  -transform... whether to transform to a hyperbolic problem or not.\n"
 15:   "  -nonhyperbolic... to test with a modified (definite) problem that is not hyperbolic.\n\n";

 17: #include <slepcpep.h>

 19: /*
 20:   This example is based on spring.c, for fixed values mu=1,tau=10,kappa=5

 22:   The transformations are based on the method proposed in [Niendorf and Voss, LAA 2010].
 23: */

 25: PetscErrorCode QEPDefiniteTransformGetMatrices(PEP,PetscBool,PetscReal,PetscReal,Mat[3]);
 26: PetscErrorCode QEPDefiniteTransformMap(PetscBool,PetscReal,PetscReal,PetscInt,PetscScalar*,PetscBool);
 27: PetscErrorCode QEPDefiniteCheckError(Mat*,PEP,PetscBool,PetscReal,PetscReal);
 28: PetscErrorCode TransformMatricesMoebius(Mat[3],MatStructure,PetscReal,PetscReal,PetscReal,PetscReal,Mat[3]);

 30: int main(int argc,char **argv)
 31: {
 32:   Mat            M,C,K,*Op,A[3],At[3],B[3]; /* problem matrices */
 33:   PEP            pep;        /* polynomial eigenproblem solver context */
 34:   ST             st;         /* spectral transformation context */
 35:   KSP            ksp;
 36:   PC             pc;
 37:   PEPProblemType type;
 38:   PetscBool      terse,transform=PETSC_FALSE,nohyp=PETSC_FALSE;
 39:   PetscInt       n=100,Istart,Iend,i,def=0,hyp;
 40:   PetscReal      muu=1,tau=10,kappa=5,inta,intb;
 41:   PetscReal      alpha,beta,xi,mu,at[2]={0.0,0.0},c=.857,s;
 42:   PetscScalar    target,targett,ats[2];

 44:   SlepcInitialize(&argc,&argv,(char*)0,help);

 46:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 47:   PetscPrintf(PETSC_COMM_WORLD,"\nPEP example that checks definite property, n=%" PetscInt_FMT "\n\n",n);

 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 51:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 53:   /* K is a tridiagonal */
 54:   MatCreate(PETSC_COMM_WORLD,&K);
 55:   MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
 56:   MatSetFromOptions(K);
 57:   MatSetUp(K);

 59:   MatGetOwnershipRange(K,&Istart,&Iend);
 60:   for (i=Istart;i<Iend;i++) {
 61:     if (i>0) MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
 62:     MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
 63:     if (i<n-1) MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
 64:   }

 66:   MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
 67:   MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);

 69:   /* C is a tridiagonal */
 70:   MatCreate(PETSC_COMM_WORLD,&C);
 71:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
 72:   MatSetFromOptions(C);
 73:   MatSetUp(C);

 75:   MatGetOwnershipRange(C,&Istart,&Iend);
 76:   for (i=Istart;i<Iend;i++) {
 77:     if (i>0) MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
 78:     MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
 79:     if (i<n-1) MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
 80:   }

 82:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 83:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 85:   /* M is a diagonal matrix */
 86:   MatCreate(PETSC_COMM_WORLD,&M);
 87:   MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
 88:   MatSetFromOptions(M);
 89:   MatSetUp(M);
 90:   MatGetOwnershipRange(M,&Istart,&Iend);
 91:   for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,muu,INSERT_VALUES);
 92:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
 93:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

 95:   PetscOptionsGetBool(NULL,NULL,"-nonhyperbolic",&nohyp,NULL);
 96:   A[0] = K; A[1] = C; A[2] = M;
 97:   if (nohyp) {
 98:     s = c*.6;
 99:     TransformMatricesMoebius(A,UNKNOWN_NONZERO_PATTERN,c,s,-s,c,At);
100:     for (i=0;i<3;i++) MatDestroy(&A[i]);
101:     Op = At;
102:   } else Op = A;

104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105:                 Create the eigensolver and solve the problem
106:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

108:   /*
109:      Create eigensolver context
110:   */
111:   PEPCreate(PETSC_COMM_WORLD,&pep);
112:   PEPSetProblemType(pep,PEP_HERMITIAN);
113:   PEPSetType(pep,PEPSTOAR);
114:   /*
115:      Set operators and set problem type
116:   */
117:   PEPSetOperators(pep,3,Op);

119:   /*
120:      Set shift-and-invert with Cholesky; select MUMPS if available
121:   */
122:   PEPGetST(pep,&st);
123:   STGetKSP(st,&ksp);
124:   KSPSetType(ksp,KSPPREONLY);
125:   KSPGetPC(ksp,&pc);
126:   PCSetType(pc,PCCHOLESKY);

128:   /*
129:      Use MUMPS if available.
130:      Note that in complex scalars we cannot use MUMPS for spectrum slicing,
131:      because MatGetInertia() is not available in that case.
132:   */
133: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
134:   PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
135:   /*
136:      Add several MUMPS options (see ex43.c for a better way of setting them in program):
137:      '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
138:   */
139:   PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12");
140: #endif

142:   /*
143:      Set solver parameters at runtime
144:   */
145:   PEPSetFromOptions(pep);

147:   PetscOptionsGetBool(NULL,NULL,"-transform",&transform,NULL);
148:   if (transform) {
149:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:                     Check if the problem is definite
151:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:     PEPCheckDefiniteQEP(pep,&xi,&mu,&def,&hyp);
153:     switch (def) {
154:       case 1:
155:         if (hyp==1) PetscPrintf(PETSC_COMM_WORLD,"Hyperbolic Problem xi=%g\n",(double)xi);
156:         else PetscPrintf(PETSC_COMM_WORLD,"Definite Problem xi=%g mu=%g\n",(double)xi,(double)mu);
157:         break;
158:       case -1:
159:         PetscPrintf(PETSC_COMM_WORLD,"Not Definite Problem\n");
160:         break;
161:       default:
162:         PetscPrintf(PETSC_COMM_WORLD,"Cannot determine definiteness\n");
163:         break;
164:     }

166:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:       Transform the QEP to have a definite inner product in the linearization
168:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:     if (def==1) {
170:       QEPDefiniteTransformGetMatrices(pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,B);
171:       PEPSetOperators(pep,3,B);
172:       PEPGetTarget(pep,&target);
173:       targett = target;
174:       QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,1,&targett,PETSC_FALSE);
175:       PEPSetTarget(pep,targett);
176:       PEPGetProblemType(pep,&type);
177:       PEPSetProblemType(pep,PEP_HYPERBOLIC);
178:       PEPSTOARGetLinearization(pep,&alpha,&beta);
179:       PEPSTOARSetLinearization(pep,1.0,0.0);
180:       PEPGetInterval(pep,&inta,&intb);
181:       if (inta!=intb) {
182:         ats[0] = inta; ats[1] = intb;
183:         QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,2,ats,PETSC_FALSE);
184:         at[0] = PetscRealPart(ats[0]); at[1] = PetscRealPart(ats[1]);
185:         if (at[0]<at[1]) PEPSetInterval(pep,at[0],at[1]);
186:         else PEPSetInterval(pep,PETSC_MIN_REAL,at[1]);
187:       }
188:     }
189:   }

191:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192:                       Solve the eigensystem
193:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194:   PEPSolve(pep);

196:   /* show detailed info unless -terse option is given by user */
197:   if (def!=1) {
198:     PetscOptionsHasName(NULL,NULL,"-terse",&terse);
199:     if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
200:     else {
201:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
202:       PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
203:       PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
204:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
205:     }
206:   } else {
207:     /* Map the solution */
208:     PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
209:     QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu);
210:     for (i=0;i<3;i++) MatDestroy(B+i);
211:   }
212:   if (at[0]>at[1]) {
213:     PEPSetInterval(pep,at[0],PETSC_MAX_REAL);
214:     PEPSolve(pep);
215:     PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
216:     /* Map the solution */
217:     QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu);
218:   }
219:   if (def==1) {
220:     PEPSetTarget(pep,target);
221:     PEPSetProblemType(pep,type);
222:     PEPSTOARSetLinearization(pep,alpha,beta);
223:     if (inta!=intb) PEPSetInterval(pep,inta,intb);
224:   }

226:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227:                     Clean up
228:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229:   PEPDestroy(&pep);
230:   for (i=0;i<3;i++) MatDestroy(Op+i);
231:   SlepcFinalize();
232:   return 0;
233: }

235: /* ------------------------------------------------------------------- */
236: /*
237:   QEPDefiniteTransformMap_Initial - map a scalar value with a certain Moebius transform

239:                    a theta + b
240:          lambda = --------------
241:                    c theta + d

243:   Input:
244:     xi,mu: real values such that Q(xi)<0 and Q(mu)>0
245:     hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
246:   Input/Output:
247:     val (array of length n)
248:     if backtransform=true returns lambda from theta, else returns theta from lambda
249: */
250: static PetscErrorCode QEPDefiniteTransformMap_Initial(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
251: {
252:   PetscInt  i;
253:   PetscReal a,b,c,d,s;

255:   if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
256:   else { a = mu; b = mu*xi-1; c = 1.0; d = xi+mu; }
257:   if (!backtransform) { s = a; a = -d; d = -s; }
258:   for (i=0;i<n;i++) {
259:     if (PetscRealPart(val[i]) >= PETSC_MAX_REAL || PetscRealPart(val[i]) <= PETSC_MIN_REAL) val[i] = a/c;
260:     else if (val[i] == -d/c) val[i] = PETSC_MAX_REAL;
261:     else val[i] = (a*val[i]+b)/(c*val[i]+d);
262:   }
263:   PetscFunctionReturn(0);
264: }

266: /* ------------------------------------------------------------------- */
267: /*
268:   QEPDefiniteTransformMap - perform the mapping if the problem is hyperbolic, otherwise
269:   modify the value of xi in advance
270: */
271: PetscErrorCode QEPDefiniteTransformMap(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
272: {
273:   PetscReal      xit;
274:   PetscScalar    alpha;

276:   xit = xi;
277:   if (!hyperbolic) {
278:     alpha = xi;
279:     QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&alpha,PETSC_FALSE);
280:     xit = PetscRealPart(alpha);
281:   }
282:   QEPDefiniteTransformMap_Initial(hyperbolic,xit,mu,n,val,backtransform);
283:   PetscFunctionReturn(0);
284: }

286: /* ------------------------------------------------------------------- */
287: /*
288:   TransformMatricesMoebius - transform the coefficient matrices of a QEP

290:   Input:
291:     A: coefficient matrices of the original QEP
292:     a,b,c,d: parameters of the Moebius transform
293:     str: structure flag for MatAXPY operations
294:   Output:
295:     B: transformed matrices
296: */
297: PetscErrorCode TransformMatricesMoebius(Mat A[3],MatStructure str,PetscReal a,PetscReal b,PetscReal c,PetscReal d,Mat B[3])
298: {
299:   PetscInt       i,k;
300:   PetscReal      cf[9];

302:   for (i=0;i<3;i++) MatDuplicate(A[2],MAT_COPY_VALUES,&B[i]);
303:   /* Ct = b*b*A+b*d*B+d*d*C */
304:   cf[0] = d*d; cf[1] = b*d; cf[2] = b*b;
305:   /* Bt = 2*a*b*A+(b*c+a*d)*B+2*c*d*C*/
306:   cf[3] = 2*c*d; cf[4] = b*c+a*d; cf[5] = 2*a*b;
307:   /* At = a*a*A+a*c*B+c*c*C */
308:   cf[6] = c*c; cf[7] = a*c; cf[8] = a*a;
309:   for (k=0;k<3;k++) {
310:     MatScale(B[k],cf[k*3+2]);
311:     for (i=0;i<2;i++) MatAXPY(B[k],cf[3*k+i],A[i],str);
312:   }
313:   PetscFunctionReturn(0);
314: }

316: /* ------------------------------------------------------------------- */
317: /*
318:   QEPDefiniteTransformGetMatrices - given a PEP of degree 2, transform the three
319:   matrices with TransformMatricesMoebius

321:   Input:
322:     pep: polynomial eigenproblem to be transformed, with Q(.) being the quadratic polynomial
323:     xi,mu: real values such that Q(xi)<0 and Q(mu)>0
324:     hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
325:   Output:
326:     T: coefficient matrices of the transformed polynomial
327: */
328: PetscErrorCode QEPDefiniteTransformGetMatrices(PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu,Mat T[3])
329: {
330:   MatStructure   str;
331:   ST             st;
332:   PetscInt       i;
333:   PetscReal      a,b,c,d;
334:   PetscScalar    xit;
335:   Mat            A[3];

337:   for (i=2;i>=0;i--) PEPGetOperators(pep,i,&A[i]);
338:   if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
339:   else {
340:     xit = xi;
341:     QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&xit,PETSC_FALSE);
342:     a = mu; b = mu*PetscRealPart(xit)-1.0; c = 1.0; d = PetscRealPart(xit)+mu;
343:   }
344:   PEPGetST(pep,&st);
345:   STGetMatStructure(st,&str);
346:   TransformMatricesMoebius(A,str,a,b,c,d,T);
347:   PetscFunctionReturn(0);
348: }

350: /* ------------------------------------------------------------------- */
351: /*
352:   Auxiliary function to compute the residual norm of an eigenpair of a QEP defined
353:   by coefficient matrices A
354: */
355: static PetscErrorCode PEPResidualNorm(Mat *A,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
356: {
357:   PetscInt       i,nmat=3;
358:   PetscScalar    vals[3];
359:   Vec            u,w;
360: #if !defined(PETSC_USE_COMPLEX)
361:   Vec            ui,wi;
362:   PetscReal      ni;
363:   PetscBool      imag;
364:   PetscScalar    ivals[3];
365: #endif

367:   u = z[0]; w = z[1];
368:   VecSet(u,0.0);
369: #if !defined(PETSC_USE_COMPLEX)
370:   ui = z[2]; wi = z[3];
371: #endif
372:   vals[0] = 1.0;
373:   vals[1] = kr;
374:   vals[2] = kr*kr-ki*ki;
375: #if !defined(PETSC_USE_COMPLEX)
376:   ivals[0] = 0.0;
377:   ivals[1] = ki;
378:   ivals[2] = 2.0*kr*ki;
379:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
380:     imag = PETSC_FALSE;
381:   else {
382:     imag = PETSC_TRUE;
383:     VecSet(ui,0.0);
384:   }
385: #endif
386:   for (i=0;i<nmat;i++) {
387:     if (vals[i]!=0.0) {
388:       MatMult(A[i],xr,w);
389:       VecAXPY(u,vals[i],w);
390:     }
391: #if !defined(PETSC_USE_COMPLEX)
392:     if (imag) {
393:       if (ivals[i]!=0 || vals[i]!=0) {
394:         MatMult(A[i],xi,wi);
395:         if (vals[i]==0) MatMult(A[i],xr,w);
396:       }
397:       if (ivals[i]!=0) {
398:         VecAXPY(u,-ivals[i],wi);
399:         VecAXPY(ui,ivals[i],w);
400:       }
401:       if (vals[i]!=0) VecAXPY(ui,vals[i],wi);
402:     }
403: #endif
404:   }
405:   VecNorm(u,NORM_2,norm);
406: #if !defined(PETSC_USE_COMPLEX)
407:   if (imag) {
408:     VecNorm(ui,NORM_2,&ni);
409:     *norm = SlepcAbsEigenvalue(*norm,ni);
410:   }
411: #endif
412:   PetscFunctionReturn(0);
413: }

415: /* ------------------------------------------------------------------- */
416: /*
417:   QEPDefiniteCheckError - check and print the residual norm of a transformed PEP

419:   Input:
420:     A: coefficient matrices of the original problem
421:     pep: solver containing the computed solution of the transformed problem
422:     xi,mu,hyperbolic: parameters used in transformation
423: */
424: PetscErrorCode QEPDefiniteCheckError(Mat *A,PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu)
425: {
426:   PetscScalar    er,ei;
427:   PetscReal      re,im,error;
428:   Vec            vr,vi,w[4];
429:   PetscInt       i,nconv;
430:   BV             bv;
431:   char           ex[30],sep[]=" ---------------------- --------------------\n";

433:   PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||");
434:   PetscPrintf(PETSC_COMM_WORLD,"%s            k             %s\n%s",sep,ex,sep);
435:   PEPGetConverged(pep,&nconv);
436:   PEPGetBV(pep,&bv);
437:   BVCreateVec(bv,w);
438:   VecDuplicate(w[0],&vr);
439:   VecDuplicate(w[0],&vi);
440:   for (i=1;i<4;i++) VecDuplicate(w[0],w+i);
441:   for (i=0;i<nconv;i++) {
442:     PEPGetEigenpair(pep,i,&er,&ei,vr,vi);
443:     QEPDefiniteTransformMap(hyperbolic,xi,mu,1,&er,PETSC_TRUE);
444:     PEPResidualNorm(A,er,0.0,vr,vi,w,&error);
445:     error /= SlepcAbsEigenvalue(er,0.0);
446: #if defined(PETSC_USE_COMPLEX)
447:     re = PetscRealPart(er);
448:     im = PetscImaginaryPart(ei);
449: #else
450:     re = er;
451:     im = ei;
452: #endif
453:     if (im!=0.0) PetscPrintf(PETSC_COMM_WORLD,"  % 9f%+9fi      %12g\n",(double)re,(double)im,(double)error);
454:     else PetscPrintf(PETSC_COMM_WORLD,"    % 12f           %12g\n",(double)re,(double)error);
455:   }
456:   PetscPrintf(PETSC_COMM_WORLD,"%s",sep);
457:   for (i=0;i<4;i++) VecDestroy(w+i);
458:   VecDestroy(&vi);
459:   VecDestroy(&vr);
460:   PetscFunctionReturn(0);
461: }

463: /*TEST

465:    testset:
466:       requires: !single
467:       args: -pep_nev 3 -nonhyperbolic -pep_target 2
468:       output_file: output/ex40_1.out
469:       filter: grep -v "Definite" | sed -e "s/iterations [0-9]\([0-9]*\)/iterations xx/g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
470:       test:
471:          suffix: 1
472:          requires: !complex
473:       test:
474:          suffix: 1_complex
475:          requires: complex !mumps
476:       test:
477:          suffix: 1_transform
478:          requires: !complex
479:          args: -transform
480:       test:
481:          suffix: 1_transform_complex
482:          requires: complex !mumps
483:          args: -transform

485: TEST*/