Actual source code: ex38.c

slepc-3.22.1 2024-10-28
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,NULL,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));

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

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

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

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

 68:   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
 69:   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));

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

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:                 Create the eigensolver and solve the problem
 82:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 84:   /*
 85:      Create eigensolver context
 86:   */
 87:   PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));

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

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

104:   /*
105:      Spectrum slicing requires STOAR
106:   */
107:   PetscCall(PEPSetType(pep,PEPSTOAR));

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

115:   PetscCall(STGetKSP(st,&ksp));
116:   PetscCall(KSPSetType(ksp,KSPPREONLY));
117:   PetscCall(KSPGetPC(ksp,&pc));
118:   PetscCall(PCSetType(pc,PCCHOLESKY));

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

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

140:   /*
141:      Set solver parameters at runtime
142:   */
143:   PetscCall(PEPSetFromOptions(pep));

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

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

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

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

198: /*TEST

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

211:    test:
212:       suffix: 2
213:       args: -pep_interval 0,2

215: TEST*/