Actual source code: test5.c
slepc-3.21.1 2024-04-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 changing MFN type.\n\n";
13: #include <slepcmfn.h>
15: int main(int argc,char **argv)
16: {
17: Mat A; /* problem matrix */
18: MFN mfn;
19: FN f;
20: PetscReal norm;
21: PetscScalar t=0.3;
22: PetscInt N,n=25,m,Istart,Iend,II,i,j;
23: PetscBool flag;
24: Vec v,y;
26: PetscFunctionBeginUser;
27: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
29: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
30: PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
31: if (!flag) m=n;
32: N = n*m;
33: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL));
34: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMatrix exponential y=exp(t*A)*e, of the 2-D Laplacian, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
36: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
37: Build the 2-D Laplacian
38: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
41: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
42: PetscCall(MatSetFromOptions(A));
44: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
45: for (II=Istart;II<Iend;II++) {
46: i = II/n; j = II-i*n;
47: if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
48: if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
49: if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
50: if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
51: PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
52: }
54: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
55: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
57: /* set v = ones(n,1) */
58: PetscCall(MatCreateVecs(A,&v,&y));
59: PetscCall(VecSet(v,1.0));
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create the solver and set various options
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: PetscCall(FNCreate(PETSC_COMM_WORLD,&f));
66: PetscCall(FNSetType(f,FNEXP));
68: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
69: PetscCall(MFNSetOperator(mfn,A));
70: PetscCall(MFNSetType(mfn,MFNEXPOKIT));
71: PetscCall(MFNSetDimensions(mfn,24));
72: PetscCall(MFNSetTolerances(mfn,1e-5,1000));
73: PetscCall(MFNSetFN(mfn,f));
74: PetscCall(MFNSetErrorIfNotConverged(mfn,PETSC_TRUE));
75: PetscCall(MFNSetFromOptions(mfn));
76: PetscCall(MFNView(mfn,NULL));
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Change MFN type and solve the problem, y=exp(t*A)*v
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
82: PetscCall(MFNSetType(mfn,MFNKRYLOV));
83: PetscCall(FNSetScale(f,t,1.0));
84: PetscCall(MFNSolve(mfn,v,y));
85: PetscCall(VecNorm(y,NORM_2,&norm));
86: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm));
88: /*
89: Free work space
90: */
91: PetscCall(MFNDestroy(&mfn));
92: PetscCall(FNDestroy(&f));
93: PetscCall(MatDestroy(&A));
94: PetscCall(VecDestroy(&v));
95: PetscCall(VecDestroy(&y));
96: PetscCall(SlepcFinalize());
97: return 0;
98: }
100: /*TEST
102: test:
104: TEST*/