Actual source code: test21.c
slepc-3.22.1 2024-10-28
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[] = "Illustrates region filtering. "
12: "Based on ex5.\n"
13: "The command line options are:\n"
14: " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
16: #include <slepceps.h>
18: /*
19: User-defined routines
20: */
21: PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
23: int main(int argc,char **argv)
24: {
25: Mat A;
26: EPS eps;
27: ST st;
28: RG rg;
29: PetscReal radius,tol=PETSC_SMALL;
30: PetscScalar target=0.5,kr,ki;
31: PetscInt N,m=15,nev,i,nconv;
32: PetscBool checkfile;
33: char filename[PETSC_MAX_PATH_LEN];
34: PetscViewer viewer;
35: char str[50];
37: PetscFunctionBeginUser;
38: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39: PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
40: N = m*(m+1)/2;
41: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov Model, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n",N,m));
42: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-target",&target,NULL));
43: PetscCall(SlepcSNPrintfScalar(str,sizeof(str),target,PETSC_FALSE));
44: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Searching closest eigenvalues to the right of %s.\n\n",str));
46: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: Compute the operator matrix that defines the eigensystem, Ax=kx
48: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
51: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
52: PetscCall(MatSetFromOptions(A));
53: PetscCall(MatMarkovModel(m,A));
55: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56: Create a standalone spectral transformation
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: PetscCall(STCreate(PETSC_COMM_WORLD,&st));
60: PetscCall(STSetType(st,STSINVERT));
61: PetscCall(STSetShift(st,target));
63: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: Create a region for filtering
65: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
67: PetscCall(RGCreate(PETSC_COMM_WORLD,&rg));
68: PetscCall(RGSetType(rg,RGELLIPSE));
69: radius = (1.0-PetscRealPart(target))/2.0;
70: PetscCall(RGEllipseSetParameters(rg,target+radius,radius,0.1));
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create the eigensolver and set various options
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
77: PetscCall(EPSSetST(eps,st));
78: PetscCall(EPSSetRG(eps,rg));
79: PetscCall(EPSSetOperators(eps,A,NULL));
80: PetscCall(EPSSetProblemType(eps,EPS_NHEP));
81: PetscCall(EPSSetTolerances(eps,tol,PETSC_CURRENT));
82: PetscCall(EPSSetTarget(eps,target));
83: PetscCall(EPSSetFromOptions(eps));
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Solve the eigensystem and display solution
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: PetscCall(EPSSolve(eps));
90: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
91: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
92: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Check file containing the eigenvalues
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: PetscCall(PetscOptionsGetString(NULL,NULL,"-checkfile",filename,sizeof(filename),&checkfile));
98: if (checkfile) {
99: #if defined(PETSC_HAVE_COMPLEX)
100: PetscComplex *eigs,eval;
102: PetscCall(EPSGetConverged(eps,&nconv));
103: PetscCall(PetscMalloc1(nconv,&eigs));
104: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
105: PetscCall(PetscViewerBinaryRead(viewer,eigs,nconv,NULL,PETSC_COMPLEX));
106: PetscCall(PetscViewerDestroy(&viewer));
107: for (i=0;i<nconv;i++) {
108: PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,NULL,NULL));
109: #if defined(PETSC_USE_COMPLEX)
110: eval = kr;
111: #else
112: eval = PetscCMPLX(kr,ki);
113: #endif
114: PetscCheck(eval==eigs[i],PETSC_COMM_WORLD,PETSC_ERR_FILE_UNEXPECTED,"Eigenvalues in the file do not match");
115: }
116: PetscCall(PetscFree(eigs));
117: #else
118: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"The -checkfile option requires C99 complex numbers");
119: #endif
120: }
122: PetscCall(EPSDestroy(&eps));
123: PetscCall(STDestroy(&st));
124: PetscCall(RGDestroy(&rg));
125: PetscCall(MatDestroy(&A));
126: PetscCall(SlepcFinalize());
127: return 0;
128: }
130: PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
131: {
132: const PetscReal cst = 0.5/(PetscReal)(m-1);
133: PetscReal pd,pu;
134: PetscInt Istart,Iend,i,j,jmax,ix=0;
136: PetscFunctionBeginUser;
137: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
138: for (i=1;i<=m;i++) {
139: jmax = m-i+1;
140: for (j=1;j<=jmax;j++) {
141: ix = ix + 1;
142: if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
143: if (j!=jmax) {
144: pd = cst*(PetscReal)(i+j-1);
145: /* north */
146: if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES));
147: else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES));
148: /* east */
149: if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES));
150: else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES));
151: }
152: /* south */
153: pu = 0.5 - cst*(PetscReal)(i+j-3);
154: if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES));
155: /* west */
156: if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES));
157: }
158: }
159: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
160: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*TEST
166: test:
167: suffix: 1
168: args: -eps_nev 4 -eps_ncv 20 -eps_view_values binary:myvalues.bin -checkfile myvalues.bin
169: output_file: output/test11_1.out
170: requires: !single c99_complex
172: TEST*/