Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : static char help[] = "Eigenproblem for the 1-D Laplacian with constraints. "
12 : "Based on ex1.\n"
13 : "The command line options are:\n"
14 : " -n <n>, where <n> = number of grid subdivisions = matrix dimension.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 1 : int main(int argc,char **argv)
19 : {
20 1 : Mat A;
21 1 : EPS eps;
22 1 : EPSType type;
23 1 : Vec *vi=NULL,*vc=NULL,t;
24 1 : PetscInt n=30,nev=4,i,j,Istart,Iend,nini=0,ncon=0,bs;
25 1 : PetscReal alpha,beta,restart;
26 1 : PetscBool flg,lock;
27 :
28 1 : PetscFunctionBeginUser;
29 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
30 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
31 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-nini",&nini,NULL));
32 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-ncon",&ncon,NULL));
33 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem, n=%" PetscInt_FMT " nini=%" PetscInt_FMT " ncon=%" PetscInt_FMT "\n\n",n,nini,ncon));
34 :
35 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36 : Compute the operator matrix that defines the eigensystem, Ax=kx
37 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38 :
39 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
40 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
41 1 : PetscCall(MatSetFromOptions(A));
42 :
43 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
44 31 : for (i=Istart;i<Iend;i++) {
45 30 : if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
46 30 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
47 30 : PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
48 : }
49 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
50 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
51 :
52 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53 : Create the eigensolver and set various options
54 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
56 1 : PetscCall(EPSSetOperators(eps,A,NULL));
57 1 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
58 1 : PetscCall(EPSSetType(eps,EPSLOBPCG));
59 1 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
60 1 : PetscCall(EPSSetConvergenceTest(eps,EPS_CONV_ABS));
61 1 : PetscCall(EPSSetDimensions(eps,nev,PETSC_DETERMINE,PETSC_DETERMINE));
62 1 : PetscCall(EPSLOBPCGSetBlockSize(eps,nev));
63 1 : PetscCall(EPSLOBPCGSetRestart(eps,0.7));
64 1 : PetscCall(EPSSetTolerances(eps,1e-8,1200));
65 1 : PetscCall(EPSSetFromOptions(eps));
66 :
67 1 : PetscCall(MatCreateVecs(A,&t,NULL));
68 1 : if (nini) {
69 0 : PetscCall(VecDuplicateVecs(t,nini,&vi));
70 0 : for (i=0;i<nini;i++) PetscCall(VecSetRandom(vi[i],NULL));
71 0 : PetscCall(EPSSetInitialSpace(eps,nini,vi));
72 : }
73 1 : if (ncon) { /* constraints are exact eigenvectors of lowest eigenvalues */
74 1 : alpha = PETSC_PI/(n+1);
75 1 : beta = PetscSqrtReal(2.0/(n+1));
76 1 : PetscCall(VecDuplicateVecs(t,ncon,&vc));
77 3 : for (i=0;i<ncon;i++) {
78 62 : for (j=0;j<n;j++) PetscCall(VecSetValue(vc[i],j,PetscSinReal(alpha*(j+1)*(i+1))*beta,INSERT_VALUES));
79 2 : PetscCall(VecAssemblyBegin(vc[i]));
80 2 : PetscCall(VecAssemblyEnd(vc[i]));
81 : }
82 1 : PetscCall(EPSSetDeflationSpace(eps,ncon,vc));
83 : }
84 :
85 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 : Solve the eigensystem
87 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 :
89 1 : PetscCall(EPSSolve(eps));
90 1 : PetscCall(EPSGetType(eps,&type));
91 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
92 1 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSLOBPCG,&flg));
93 1 : if (flg) {
94 1 : PetscCall(EPSLOBPCGGetLocking(eps,&lock));
95 1 : if (lock) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using soft locking\n"));
96 1 : PetscCall(EPSLOBPCGGetRestart(eps,&restart));
97 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," LOBPCG Restart parameter=%.4g\n",(double)restart));
98 1 : PetscCall(EPSLOBPCGGetBlockSize(eps,&bs));
99 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," LOBPCG Block size=%" PetscInt_FMT "\n",bs));
100 : }
101 :
102 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103 : Display solution and clean up
104 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105 :
106 1 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
107 1 : PetscCall(EPSDestroy(&eps));
108 1 : PetscCall(MatDestroy(&A));
109 1 : PetscCall(VecDestroyVecs(nini,&vi));
110 1 : PetscCall(VecDestroyVecs(ncon,&vc));
111 1 : PetscCall(VecDestroy(&t));
112 1 : PetscCall(SlepcFinalize());
113 : return 0;
114 : }
115 :
116 : /*TEST
117 :
118 : testset:
119 : args: -ncon 2
120 : output_file: output/test24_1.out
121 : test:
122 : suffix: 1
123 : requires: !single
124 : test:
125 : suffix: 1_cuda
126 : args: -mat_type aijcusparse
127 : requires: cuda !single
128 : test:
129 : suffix: 1_hip
130 : args: -mat_type aijhipsparse
131 : requires: hip !single
132 :
133 : TEST*/
|
