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[] = "Test MFN interface functions.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n\n";
14 :
15 : #include <slepcmfn.h>
16 :
17 3 : int main(int argc,char **argv)
18 : {
19 3 : Mat A,B;
20 3 : MFN mfn;
21 3 : FN f;
22 3 : MFNConvergedReason reason;
23 3 : MFNType type;
24 3 : PetscReal norm,tol;
25 3 : Vec v,y;
26 3 : PetscInt N,n=4,Istart,Iend,i,j,II,ncv,its,maxit;
27 3 : PetscBool flg,testprefix=PETSC_FALSE;
28 3 : const char *prefix;
29 3 : PetscViewerAndFormat *vf;
30 :
31 3 : PetscFunctionBeginUser;
32 3 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
33 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
34 3 : N = n*n;
35 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root of Laplacian y=sqrt(A)*e_1, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,n));
36 3 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_prefix",&testprefix,NULL));
37 :
38 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39 : Compute the discrete 2-D Laplacian, A
40 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41 :
42 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
43 3 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
44 3 : PetscCall(MatSetFromOptions(A));
45 :
46 3 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
47 51 : for (II=Istart;II<Iend;II++) {
48 48 : i = II/n; j = II-i*n;
49 48 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
50 48 : if (i<n-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
51 48 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
52 48 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
53 48 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
54 : }
55 :
56 3 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
57 3 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
58 3 : PetscCall(MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE));
59 :
60 3 : PetscCall(MatCreateVecs(A,&v,&y));
61 3 : PetscCall(VecSetValue(v,0,1.0,INSERT_VALUES));
62 3 : PetscCall(VecAssemblyBegin(v));
63 3 : PetscCall(VecAssemblyEnd(v));
64 :
65 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 : Create the solver, set the matrix and the function
67 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68 3 : PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
69 3 : PetscCall(MFNSetOperator(mfn,A));
70 3 : PetscCall(MFNGetFN(mfn,&f));
71 3 : PetscCall(FNSetType(f,FNSQRT));
72 :
73 3 : PetscCall(MFNSetType(mfn,MFNKRYLOV));
74 3 : PetscCall(MFNGetType(mfn,&type));
75 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Type set to %s\n",type));
76 :
77 : /* test prefix usage */
78 3 : if (testprefix) {
79 1 : PetscCall(MFNSetOptionsPrefix(mfn,"check_"));
80 1 : PetscCall(MFNAppendOptionsPrefix(mfn,"myprefix_"));
81 1 : PetscCall(MFNGetOptionsPrefix(mfn,&prefix));
82 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," MFN prefix is currently: %s\n",prefix));
83 : }
84 :
85 : /* test some interface functions */
86 3 : PetscCall(MFNGetOperator(mfn,&B));
87 3 : PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD));
88 3 : PetscCall(MFNSetTolerances(mfn,1e-4,500));
89 3 : PetscCall(MFNSetDimensions(mfn,6));
90 3 : PetscCall(MFNSetErrorIfNotConverged(mfn,PETSC_TRUE));
91 : /* test monitors */
92 3 : PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf));
93 3 : PetscCall(MFNMonitorSet(mfn,(PetscErrorCode (*)(MFN,PetscInt,PetscReal,void*))MFNMonitorDefault,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy));
94 : /* PetscCall(MFNMonitorCancel(mfn)); */
95 3 : PetscCall(MFNSetFromOptions(mfn));
96 :
97 : /* query properties and print them */
98 3 : PetscCall(MFNGetTolerances(mfn,&tol,&maxit));
99 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Tolerance: %g, max iterations: %" PetscInt_FMT "\n",(double)tol,maxit));
100 3 : PetscCall(MFNGetDimensions(mfn,&ncv));
101 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv));
102 3 : PetscCall(MFNGetErrorIfNotConverged(mfn,&flg));
103 3 : if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Erroring out if convergence fails\n"));
104 :
105 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 : Solve y=sqrt(A)*v
107 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108 :
109 3 : PetscCall(MFNSolve(mfn,v,y));
110 3 : PetscCall(MFNGetConvergedReason(mfn,&reason));
111 3 : PetscCall(MFNGetIterationNumber(mfn,&its));
112 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Finished - converged reason = %d\n",(int)reason));
113 : /* PetscCall(PetscPrintf(PETSC_COMM_WORLD," its = %" PetscInt_FMT "\n",its)); */
114 3 : PetscCall(VecNorm(y,NORM_2,&norm));
115 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," sqrt(A)*v has norm %g\n",(double)norm));
116 :
117 : /*
118 : Free work space
119 : */
120 3 : PetscCall(MFNDestroy(&mfn));
121 3 : PetscCall(MatDestroy(&A));
122 3 : PetscCall(VecDestroy(&v));
123 3 : PetscCall(VecDestroy(&y));
124 3 : PetscCall(SlepcFinalize());
125 : return 0;
126 : }
127 :
128 : /*TEST
129 :
130 : testset:
131 : args: -mfn_monitor_cancel -mfn_converged_reason -mfn_view -log_exclude mfn,bv,fn
132 : output_file: output/test3_1.out
133 : test:
134 : suffix: 1
135 : test:
136 : suffix: 1_x
137 : args: -mfn_monitor draw::draw_lg -draw_virtual
138 : requires: x
139 :
140 : test:
141 : suffix: 2
142 : args: -test_prefix -check_myprefix_mfn_monitor
143 : filter: sed -e "s/estimate [0-9]\.[0-9]*e[+-]\([0-9]*\)/estimate (removed)/g" | sed -e "s/4.0[0-9]*e-10/4.03e-10/"
144 :
145 : TEST*/
|
