Actual source code: test7.c
slepc-3.18.2 2023-01-26
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[] = "Test interface functions of spectrum-slicing STOAR.\n\n"
12: "This is based on ex38.c. 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: PetscBool showinertia=PETSC_TRUE,lock,detect,checket;
25: PetscInt n=100,Istart,Iend,i,*inertias,ns,nev,ncv,mpd;
26: PetscReal mu=1.0,tau=10.0,kappa=5.0,int0,int1,*shifts;
29: SlepcInitialize(&argc,&argv,(char*)0,help);
31: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
32: PetscOptionsGetBool(NULL,NULL,"-showinertia",&showinertia,NULL);
33: PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\n\n",n);
35: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
37: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
39: /* K is a tridiagonal */
40: MatCreate(PETSC_COMM_WORLD,&K);
41: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
42: MatSetFromOptions(K);
43: MatSetUp(K);
45: MatGetOwnershipRange(K,&Istart,&Iend);
46: for (i=Istart;i<Iend;i++) {
47: if (i>0) MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
48: MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
49: if (i<n-1) MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
50: }
52: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
53: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
55: /* C is a tridiagonal */
56: MatCreate(PETSC_COMM_WORLD,&C);
57: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
58: MatSetFromOptions(C);
59: MatSetUp(C);
61: MatGetOwnershipRange(C,&Istart,&Iend);
62: for (i=Istart;i<Iend;i++) {
63: if (i>0) MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
64: MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
65: if (i<n-1) MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
66: }
68: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
71: /* M is a diagonal matrix */
72: MatCreate(PETSC_COMM_WORLD,&M);
73: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
74: MatSetFromOptions(M);
75: MatSetUp(M);
76: MatGetOwnershipRange(M,&Istart,&Iend);
77: for (i=Istart;i<Iend;i++) MatSetValue(M,i,i,mu,INSERT_VALUES);
78: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Create the eigensolver and solve the problem
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PEPCreate(PETSC_COMM_WORLD,&pep);
86: A[0] = K; A[1] = C; A[2] = M;
87: PEPSetOperators(pep,3,A);
88: PEPSetProblemType(pep,PEP_HYPERBOLIC);
89: PEPSetType(pep,PEPSTOAR);
91: /*
92: Set interval and other settings for spectrum slicing
93: */
94: int0 = -11.3;
95: int1 = -9.5;
96: PEPSetInterval(pep,int0,int1);
97: PEPSetWhichEigenpairs(pep,PEP_ALL);
98: PEPGetST(pep,&st);
99: STSetType(st,STSINVERT);
100: STGetKSP(st,&ksp);
101: KSPSetType(ksp,KSPPREONLY);
102: KSPGetPC(ksp,&pc);
103: PCSetType(pc,PCCHOLESKY);
105: /*
106: Test interface functions of STOAR solver
107: */
108: PEPSTOARGetDetectZeros(pep,&detect);
109: PetscPrintf(PETSC_COMM_WORLD," Detect zeros before changing = %d",(int)detect);
110: PEPSTOARSetDetectZeros(pep,PETSC_TRUE);
111: PEPSTOARGetDetectZeros(pep,&detect);
112: PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)detect);
114: PEPSTOARGetLocking(pep,&lock);
115: PetscPrintf(PETSC_COMM_WORLD," Locking flag before changing = %d",(int)lock);
116: PEPSTOARSetLocking(pep,PETSC_TRUE);
117: PEPSTOARGetLocking(pep,&lock);
118: PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)lock);
120: PEPSTOARGetCheckEigenvalueType(pep,&checket);
121: PetscPrintf(PETSC_COMM_WORLD," Check eigenvalue type flag before changing = %d",(int)checket);
122: PEPSTOARSetCheckEigenvalueType(pep,PETSC_FALSE);
123: PEPSTOARGetCheckEigenvalueType(pep,&checket);
124: PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)checket);
126: PEPSTOARGetDimensions(pep,&nev,&ncv,&mpd);
127: PetscPrintf(PETSC_COMM_WORLD," Sub-solve dimensions before changing = [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]",nev,ncv,mpd);
128: PEPSTOARSetDimensions(pep,30,60,60);
129: PEPSTOARGetDimensions(pep,&nev,&ncv,&mpd);
130: PetscPrintf(PETSC_COMM_WORLD," ... changed to [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]\n",nev,ncv,mpd);
132: PEPSetFromOptions(pep);
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Compute all eigenvalues in interval and display info
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PEPSetUp(pep);
139: PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
140: PetscPrintf(PETSC_COMM_WORLD," Inertias (after setup):\n");
141: for (i=0;i<ns;i++) PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]);
142: PetscFree(shifts);
143: PetscFree(inertias);
145: PEPSolve(pep);
146: PEPGetDimensions(pep,&nev,NULL,NULL);
147: PEPGetInterval(pep,&int0,&int1);
148: PetscPrintf(PETSC_COMM_WORLD," Found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1);
150: if (showinertia) {
151: PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
152: PetscPrintf(PETSC_COMM_WORLD," Used %" PetscInt_FMT " shifts (inertia):\n",ns);
153: for (i=0;i<ns;i++) PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]);
154: PetscFree(shifts);
155: PetscFree(inertias);
156: }
158: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Clean up
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
163: PEPDestroy(&pep);
164: MatDestroy(&M);
165: MatDestroy(&C);
166: MatDestroy(&K);
167: SlepcFinalize();
168: return 0;
169: }
171: /*TEST
173: test:
174: requires: !single
175: args: -showinertia 0
177: TEST*/