Actual source code: ex38.c

slepc-3.20.1 2023-11-27
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:   PetscFunctionBeginUser;
 30:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

 32:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 33:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
 34:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
 35:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\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:   PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
 43:   PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
 44:   PetscCall(MatSetFromOptions(K));
 45:   PetscCall(MatSetUp(K));

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

 54:   PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
 55:   PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));

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

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

 70:   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
 71:   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));

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

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

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

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

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

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

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

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

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

137:      Note: depending on the interval, it may be necessary also to increase the workspace:
138:      '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
139:   */
140:   PetscCall(PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12"));
141: #endif

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

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:                       Solve the eigensystem
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151:   PetscCall(PEPSetUp(pep));
152:   if (show) {
153:     PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
154:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
155:     for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
156:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
157:     PetscCall(PetscFree(shifts));
158:     PetscCall(PetscFree(inertias));
159:   }
160:   PetscCall(PEPSolve(pep));
161:   if (show && !terse) {
162:     PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
163:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
164:     for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g  Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
165:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
166:     PetscCall(PetscFree(shifts));
167:     PetscCall(PetscFree(inertias));
168:   }

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

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

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

201: /*TEST

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

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

218: TEST*/