Actual source code: ex38.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[] = "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;

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

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

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

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

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

 53:   MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
 54:   MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);

 56:   /* C is a tridiagonal */
 57:   MatCreate(PETSC_COMM_WORLD,&C);
 58:   MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
 59:   MatSetFromOptions(C);
 60:   MatSetUp(C);

 62:   MatGetOwnershipRange(C,&Istart,&Iend);
 63:   for (i=Istart;i<Iend;i++) {
 64:     if (i>0) MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
 65:     MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
 66:     if (i<n-1) MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
 67:   }

 69:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 70:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

 72:   /* M is a diagonal matrix */
 73:   MatCreate(PETSC_COMM_WORLD,&M);
 74:   MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
 75:   MatSetFromOptions(M);
 76:   MatSetUp(M);
 77:   MatGetOwnershipRange(M,&Istart,&Iend);
 78:   for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,mu,INSERT_VALUES);
 79:   MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:                 Create the eigensolver and solve the problem
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   /*
 87:      Create eigensolver context
 88:   */
 89:   PEPCreate(PETSC_COMM_WORLD,&pep);

 91:   /*
 92:      Set operators and set problem type
 93:   */
 94:   A[0] = K; A[1] = C; A[2] = M;
 95:   PEPSetOperators(pep,3,A);
 96:   PEPSetProblemType(pep,PEP_HYPERBOLIC);

 98:   /*
 99:      Set interval for spectrum slicing
100:   */
101:   inta = -11.3;
102:   intb = -9.5;
103:   PEPSetInterval(pep,inta,intb);
104:   PEPSetWhichEigenpairs(pep,PEP_ALL);

106:   /*
107:      Spectrum slicing requires STOAR
108:   */
109:   PEPSetType(pep,PEPSTOAR);

111:   /*
112:      Set shift-and-invert with Cholesky; select MUMPS if available
113:   */
114:   PEPGetST(pep,&st);
115:   STSetType(st,STSINVERT);

117:   STGetKSP(st,&ksp);
118:   KSPSetType(ksp,KSPPREONLY);
119:   KSPGetPC(ksp,&pc);
120:   PCSetType(pc,PCCHOLESKY);

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

136:      Note: depending on the interval, it may be necessary also to increase the workspace:
137:      '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
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:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:                       Solve the eigensystem
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150:   PEPSetUp(pep);
151:   if (show) {
152:     PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
153:     PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n");
154:     for (i=0;i<ns;i++) PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]);
155:     PetscPrintf(PETSC_COMM_WORLD,"\n");
156:     PetscFree(shifts);
157:     PetscFree(inertias);
158:   }
159:   PEPSolve(pep);
160:   if (show && !terse) {
161:     PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
162:     PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n");
163:     for (i=0;i<ns;i++) PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]);
164:     PetscPrintf(PETSC_COMM_WORLD,"\n");
165:     PetscFree(shifts);
166:     PetscFree(inertias);
167:   }

169:   /*
170:      Show eigenvalues in interval and print solution
171:   */
172:   PEPGetType(pep,&type);
173:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
174:   PEPGetDimensions(pep,&nev,NULL,NULL);
175:   PEPGetInterval(pep,&inta,&intb);
176:   PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb);

178:   /*
179:      Show detailed info unless -terse option is given by user
180:    */
181:   if (terse) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
182:   else {
183:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
184:     PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
185:     PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
186:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
187:   }

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:                     Clean up
191:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192:   PEPDestroy(&pep);
193:   MatDestroy(&M);
194:   MatDestroy(&C);
195:   MatDestroy(&K);
196:   SlepcFinalize();
197:   return 0;
198: }

200: /*TEST

202:    testset:
203:       requires: !single
204:       args: -show_inertias -terse
205:       output_file: output/ex38_1.out
206:       test:
207:          suffix: 1
208:          requires: !complex
209:       test:
210:          suffix: 1_complex
211:          requires: complex !mumps

213:    test:
214:       suffix: 2
215:       args: -pep_interval 0,2

217: TEST*/