Actual source code: ex38.c

slepc-3.18.2 2023-01-26
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;

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

 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=%" 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:   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) MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
 50:     MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
 51:     if (i<n-1) MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
 52:   }

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

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

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

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

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

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

 87:   /*
 88:      Create eigensolver context
 89:   */
 90:   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:   PEPSetOperators(pep,3,A);
 97:   PEPSetProblemType(pep,PEP_HYPERBOLIC);

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

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

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

118:   STGetKSP(st,&ksp);
119:   KSPSetType(ksp,KSPPREONLY);
120:   KSPGetPC(ksp,&pc);
121:   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:   PEPSTOARSetDetectZeros(pep,PETSC_TRUE);  /* enforce zero detection */
130:   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:   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:   PEPSetFromOptions(pep);

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

170:   /*
171:      Show eigenvalues in interval and print solution
172:   */
173:   PEPGetType(pep,&type);
174:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
175:   PEPGetDimensions(pep,&nev,NULL,NULL);
176:   PEPGetInterval(pep,&inta,&intb);
177:   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) PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
183:   else {
184:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
185:     PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
186:     PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
187:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
188:   }

190:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191:                     Clean up
192:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193:   PEPDestroy(&pep);
194:   MatDestroy(&M);
195:   MatDestroy(&C);
196:   MatDestroy(&K);
197:   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*/