Actual source code: ex58.c
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[] = "Linear Response eigenvalue problem.\n\n"
12: "The command line options are:\n"
13: " -n <n>, where <n> = dimension of the blocks.\n"
14: " -reduced <0/1>, to use the reduced form of the LREP.\n"
15: " -A <filename>, A matrix.\n"
16: " -B <filename>, B matrix.\n\n";
18: #include <slepceps.h>
20: /*
21: This example is similar to ex55.c, but with a real matrix
23: H = [ A B
24: -B -A ],
26: where A,B are real symmetric with the following Toeplitz structure:
28: A = pentadiag{a,b,c,b,a}
29: B = tridiag{b,d,b}
31: where a,b,c,d are real values.
33: Alternatively, use -A and -B command-line options to load matrices from file.
34: */
36: int main(int argc,char **argv)
37: {
38: Mat H,A,B,K,M; /* problem matrices */
39: EPS eps; /* eigenproblem solver context */
40: PetscScalar a,b,c,d;
41: PetscReal lev;
42: PetscInt n=24,Istart,Iend,i,nconv;
43: PetscBool flg,terse,checkorthog,nest=PETSC_FALSE,reduced=PETSC_FALSE;
44: Vec t,*x,*y;
45: PetscViewer viewer;
46: char filename[PETSC_MAX_PATH_LEN];
48: PetscFunctionBeginUser;
49: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
50: PetscCheck(!PetscDefined(USE_COMPLEX),PETSC_COMM_SELF,PETSC_ERR_SUP,"This example cannot be used with complex scalars");
52: PetscCall(PetscOptionsGetBool(NULL,NULL,"-reduced",&reduced,NULL));
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Compute (or load) the problem matrices A and B
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58: PetscCall(PetscOptionsHasName(NULL,NULL,"-A",&flg));
60: if (flg) {
62: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nLinear Response Eigenvalue Problem stored in file%s\n\n",reduced?" (reduced)":""));
63: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from a binary file...\n"));
65: PetscCall(PetscOptionsGetString(NULL,NULL,"-A",filename,sizeof(filename),&flg));
66: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
67: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
68: PetscCall(MatSetFromOptions(A));
69: PetscCall(MatLoad(A,viewer));
70: PetscCall(PetscViewerDestroy(&viewer));
72: PetscCall(PetscOptionsGetString(NULL,NULL,"-B",filename,sizeof(filename),&flg));
73: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate file for both A and B matrices");
74: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
75: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
76: PetscCall(MatSetFromOptions(B));
77: PetscCall(MatLoad(B,viewer));
78: PetscCall(PetscViewerDestroy(&viewer));
80: } else {
82: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
83: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nLinear Response Eigenvalue Problem, n=%" PetscInt_FMT "%s\n\n",n,reduced?" (reduced)":""));
85: a = -0.1;
86: b = 1.0;
87: c = 4.5;
88: d = 2.0;
90: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
91: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
92: PetscCall(MatSetFromOptions(A));
94: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
95: PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,n,n));
96: PetscCall(MatSetFromOptions(B));
98: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
99: for (i=Istart;i<Iend;i++) {
100: if (i>1) PetscCall(MatSetValue(A,i,i-2,a,INSERT_VALUES));
101: if (i>0) PetscCall(MatSetValue(A,i,i-1,b,INSERT_VALUES));
102: PetscCall(MatSetValue(A,i,i,c,INSERT_VALUES));
103: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,b,INSERT_VALUES));
104: if (i<n-2) PetscCall(MatSetValue(A,i,i+2,a,INSERT_VALUES));
105: }
107: PetscCall(MatGetOwnershipRange(B,&Istart,&Iend));
108: for (i=Istart;i<Iend;i++) {
109: if (i>0) PetscCall(MatSetValue(B,i,i-1,b,INSERT_VALUES));
110: PetscCall(MatSetValue(B,i,i,d,INSERT_VALUES));
111: if (i<n-1) PetscCall(MatSetValue(B,i,i+1,b,INSERT_VALUES));
112: }
114: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
115: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
116: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
117: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
118: }
120: if (reduced) {
122: PetscInt ma,na,Ma,Na;
123: const PetscScalar scal[] = { 1.0, -1.0 };
125: PetscCall(MatGetSize(A,&Ma,&Na));
126: PetscCall(MatGetLocalSize(A,&ma,&na));
128: /* K = A-B */
129: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&K));
130: PetscCall(MatSetSizes(K,ma,na,Ma,Na));
131: PetscCall(MatSetType(K,MATCOMPOSITE));
132: PetscCall(MatCompositeAddMat(K,A));
133: PetscCall(MatCompositeAddMat(K,B));
134: PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
135: PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
136: PetscCall(MatCompositeSetScalings(K,scal));
138: /* M = A+B */
139: PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&M));
140: PetscCall(MatSetSizes(M,ma,na,Ma,Na));
141: PetscCall(MatSetType(M,MATCOMPOSITE));
142: PetscCall(MatCompositeAddMat(M,A));
143: PetscCall(MatCompositeAddMat(M,B));
144: PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
145: PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
147: PetscCall(MatCreateLREP(K,M,reduced,&H));
148: PetscCall(MatDestroy(&K));
149: PetscCall(MatDestroy(&M));
151: } else PetscCall(MatCreateLREP(A,B,reduced,&H));
153: /* if you prefer, set the vector type so that MatCreateVecs() returns nested vectors */
154: PetscCall(PetscOptionsGetBool(NULL,NULL,"-nest",&nest,NULL));
155: if (nest) PetscCall(MatNestSetVecType(H,VECNEST));
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Create the eigensolver and set various options
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
162: PetscCall(EPSSetOperators(eps,H,NULL));
163: PetscCall(EPSSetProblemType(eps,EPS_LREP));
164: PetscCall(EPSSetFromOptions(eps));
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Solve the eigensystem and display solution
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: PetscCall(EPSSolve(eps));
172: /* show detailed info unless -terse option is given by user */
173: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
174: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
175: else {
176: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
177: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
178: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
179: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
180: }
182: /* check bi-orthogonality */
183: PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
184: PetscCall(EPSGetConverged(eps,&nconv));
185: if (checkorthog && nconv>0) {
186: PetscCall(MatCreateVecs(H,&t,NULL));
187: PetscCall(VecDuplicateVecs(t,nconv,&x));
188: PetscCall(VecDuplicateVecs(t,nconv,&y));
189: for (i=0;i<nconv;i++) {
190: PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
191: PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
192: }
193: PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
194: if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
195: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
196: PetscCall(VecDestroy(&t));
197: PetscCall(VecDestroyVecs(nconv,&x));
198: PetscCall(VecDestroyVecs(nconv,&y));
199: }
201: PetscCall(EPSDestroy(&eps));
202: PetscCall(MatDestroy(&A));
203: PetscCall(MatDestroy(&B));
204: PetscCall(MatDestroy(&H));
205: PetscCall(SlepcFinalize());
206: return 0;
207: }
209: /*TEST
211: build:
212: requires: !complex
214: testset:
215: args: -eps_nev 4 -eps_ncv 16 -terse -checkorthog -reduced {{0 1}} -nest {{0 1}}
216: nsize: {{1 2}}
217: filter: sed -e "s/ (reduced)//"
218: output_file: output/ex58_1.out
219: test:
220: suffix: 1
221: test:
222: suffix: 1_dense
223: args: -mat_type dense
224: test:
225: suffix: 1_cuda
226: args: -mat_type aijcusparse
227: requires: cuda
228: test:
229: suffix: 1_hip
230: args: -mat_type aijhipsparse
231: requires: hip
233: test:
234: args: -n 90 -eps_threshold_absolute 2.4 -eps_ncv 10 -terse -checkorthog -reduced {{0 1}}
235: filter: sed -e "s/ (reduced)//"
236: suffix: 2
238: TEST*/