Actual source code: ex38.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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[] = "Spectrum slicing on quadratic symmetric eigenproblem.\n\n"
 12:   "The command line options are:\n"
 13:   "  -n <n> ... dimension of the matrices.\n\n";

 15: #include <slepcpep.h>

 17: int main(int argc,char **argv)
 18: {
 19:   Mat            M,C,K,A[3]; /* problem matrices */
 20:   PEP            pep;        /* polynomial eigenproblem solver context */
 21:   ST             st;         /* spectral transformation context */
 22:   KSP            ksp;
 23:   PC             pc;
 24:   PEPType        type;
 25:   PetscBool      show=PETSC_FALSE,terse;
 26:   PetscInt       n=100,Istart,Iend,nev,i,*inertias,ns;
 27:   PetscReal      mu=1,tau=10,kappa=5,inta,intb,*shifts;

 30:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 32:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 33:   PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL);
 34:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
 35:   PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%D\n\n",n);

 37:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 38:      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 39:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 41:   /* K is a tridiagonal */
 42:   MatCreate(PETSC_COMM_WORLD,&K);
 43:   MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
 44:   MatSetFromOptions(K);
 45:   MatSetUp(K);

 47:   MatGetOwnershipRange(K,&Istart,&Iend);
 48:   for (i=Istart;i<Iend;i++) {
 49:     if (i>0) {
 50:       MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
 51:     }
 52:     MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
 53:     if (i<n-1) {
 54:       MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
 55:     }
 56:   }

 58:   MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);

 61:   /* C is a tridiagonal */
 62:   MatCreate(PETSC_COMM_WORLD,&C);
 63:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
 64:   MatSetFromOptions(C);
 65:   MatSetUp(C);

 67:   MatGetOwnershipRange(C,&Istart,&Iend);
 68:   for (i=Istart;i<Iend;i++) {
 69:     if (i>0) {
 70:       MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
 71:     }
 72:     MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
 73:     if (i<n-1) {
 74:       MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
 75:     }
 76:   }

 78:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 79:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

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

 93:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 94:                 Create the eigensolver and solve the problem
 95:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 97:   /*
 98:      Create eigensolver context
 99:   */
100:   PEPCreate(PETSC_COMM_WORLD,&pep);

102:   /*
103:      Set operators and set problem type
104:   */
105:   A[0] = K; A[1] = C; A[2] = M;
106:   PEPSetOperators(pep,3,A);
107:   PEPSetProblemType(pep,PEP_HYPERBOLIC);

109:   /*
110:      Set interval for spectrum slicing
111:   */
112:   inta = -11.3;
113:   intb = -9.5;
114:   PEPSetInterval(pep,inta,intb);
115:   PEPSetWhichEigenpairs(pep,PEP_ALL);

117:   /*
118:      Spectrum slicing requires STOAR
119:   */
120:   PEPSetType(pep,PEPSTOAR);

122:   /*
123:      Set shift-and-invert with Cholesky; select MUMPS if available
124:   */
125:   PEPGetST(pep,&st);
126:   STSetType(st,STSINVERT);

128:   STGetKSP(st,&ksp);
129:   KSPSetType(ksp,KSPPREONLY);
130:   KSPGetPC(ksp,&pc);
131:   PCSetType(pc,PCCHOLESKY);

133:   /*
134:      Use MUMPS if available.
135:      Note that in complex scalars we cannot use MUMPS for spectrum slicing,
136:      because MatGetInertia() is not available in that case.
137:   */
138: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
139:   PEPSTOARSetDetectZeros(pep,PETSC_TRUE);  /* enforce zero detection */
140:   PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
141:   /*
142:      Add several MUMPS options (see ex43.c for a better way of setting them in program):
143:      '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
144:      '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
145:      '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)

147:      Note: depending on the interval, it may be necessary also to increase the workspace:
148:      '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
149:   */
150:   PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12");
151: #endif

153:   /*
154:      Set solver parameters at runtime
155:   */
156:   PEPSetFromOptions(pep);

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:                       Solve the eigensystem
160:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   PEPSetUp(pep);
162:   if (show) {
163:     PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
164:     PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n");
165:     for (i=0;i<ns;i++) { PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %D \n",(double)shifts[i],inertias[i]); }
166:     PetscPrintf(PETSC_COMM_WORLD,"\n");
167:     PetscFree(shifts);
168:     PetscFree(inertias);
169:   }
170:   PEPSolve(pep);
171:   if (show && !terse) {
172:     PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
173:     PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n");
174:     for (i=0;i<ns;i++) { PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %D \n",(double)shifts[i],inertias[i]); }
175:     PetscPrintf(PETSC_COMM_WORLD,"\n");
176:     PetscFree(shifts);
177:     PetscFree(inertias);
178:   }

180:   /*
181:      Show eigenvalues in interval and print solution
182:   */
183:   PEPGetType(pep,&type);
184:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
185:   PEPGetDimensions(pep,&nev,NULL,NULL);
186:   PEPGetInterval(pep,&inta,&intb);
187:   PetscPrintf(PETSC_COMM_WORLD," %D eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb);

189:   /*
190:      Show detailed info unless -terse option is given by user
191:    */
192:   if (terse) {
193:     PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
194:   } else {
195:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
196:     PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
197:     PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
198:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
199:   }

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:                     Clean up
203:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204:   PEPDestroy(&pep);
205:   MatDestroy(&M);
206:   MatDestroy(&C);
207:   MatDestroy(&K);
208:   SlepcFinalize();
209:   return ierr;
210: }

212: /*TEST

214:    testset:
215:       requires: !single
216:       args: -show_inertias -terse
217:       output_file: output/ex38_1.out
218:       test:
219:          suffix: 1
220:          requires: !complex
221:       test:
222:          suffix: 1_complex
223:          requires: complex !mumps

225:    test:
226:       suffix: 2
227:       args: -pep_interval 0,2

229: TEST*/