Line | Branch | Exec | Source |
---|---|---|---|
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[] = "Computes exp(t*A)*v for a matrix associated with a Markov model.\n\n" | ||
12 | "The command line options are:\n" | ||
13 | " -t <t>, where <t> = time parameter (multiplies the matrix).\n" | ||
14 | " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n" | ||
15 | "To draw the solution run with -mfn_view_solution draw -draw_pause -1\n\n"; | ||
16 | |||
17 | #include <slepcmfn.h> | ||
18 | |||
19 | /* | ||
20 | User-defined routines | ||
21 | */ | ||
22 | PetscErrorCode MatMarkovModel(PetscInt m,Mat A); | ||
23 | |||
24 | 10 | int main(int argc,char **argv) | |
25 | { | ||
26 | 10 | Mat A; /* problem matrix */ | |
27 | 10 | MFN mfn; | |
28 | 10 | FN f; | |
29 | 10 | PetscReal tol,norm; | |
30 | 10 | PetscScalar t=2.0; | |
31 | 10 | Vec v,y; | |
32 | 10 | PetscInt N,m=15,ncv,maxit,its; | |
33 | 10 | MFNConvergedReason reason; | |
34 | |||
35 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | PetscFunctionBeginUser; |
36 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(SlepcInitialize(&argc,&argv,NULL,help)); |
37 | |||
38 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL)); |
39 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL)); |
40 | 10 | N = m*(m+1)/2; | |
41 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nMarkov y=exp(t*A)*e_1, N=%" PetscInt_FMT " (m=%" PetscInt_FMT ")\n\n",N,m)); |
42 | |||
43 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
44 | Compute the transition probability matrix, A | ||
45 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
46 | |||
47 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); |
48 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N)); |
49 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatSetFromOptions(A)); |
50 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatMarkovModel(m,A)); |
51 | |||
52 | /* set v = e_1 */ | ||
53 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatCreateVecs(A,NULL,&y)); |
54 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatCreateVecs(A,NULL,&v)); |
55 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecSetValue(v,0,1.0,INSERT_VALUES)); |
56 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecAssemblyBegin(v)); |
57 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecAssemblyEnd(v)); |
58 | |||
59 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
60 | Create the solver and set various options | ||
61 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
62 | /* | ||
63 | Create matrix function solver context | ||
64 | */ | ||
65 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn)); |
66 | |||
67 | /* | ||
68 | Set operator matrix, the function to compute, and other options | ||
69 | */ | ||
70 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSetOperator(mfn,A)); |
71 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetFN(mfn,&f)); |
72 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(FNSetType(f,FNEXP)); |
73 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(FNSetScale(f,t,1.0)); |
74 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSetTolerances(mfn,1e-07,PETSC_CURRENT)); |
75 | |||
76 | /* | ||
77 | Set solver parameters at runtime | ||
78 | */ | ||
79 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSetFromOptions(mfn)); |
80 | |||
81 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
82 | Solve the problem, y=exp(t*A)*v | ||
83 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
84 | |||
85 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNSolve(mfn,v,y)); |
86 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetConvergedReason(mfn,&reason)); |
87 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
10 | PetscCheck(reason>=0,PETSC_COMM_WORLD,PETSC_ERR_CONV_FAILED,"Solver did not converge"); |
88 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecNorm(y,NORM_2,&norm)); |
89 | |||
90 | /* | ||
91 | Optional: Get some information from the solver and display it | ||
92 | */ | ||
93 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetIterationNumber(mfn,&its)); |
94 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its)); |
95 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetDimensions(mfn,&ncv)); |
96 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %" PetscInt_FMT "\n",ncv)); |
97 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNGetTolerances(mfn,&tol,&maxit)); |
98 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit)); |
99 | |||
100 | /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | ||
101 | Display solution and clean up | ||
102 | - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ | ||
103 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm)); |
104 | |||
105 | /* | ||
106 | Free work space | ||
107 | */ | ||
108 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MFNDestroy(&mfn)); |
109 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatDestroy(&A)); |
110 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecDestroy(&v)); |
111 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(VecDestroy(&y)); |
112 |
3/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
10 | PetscCall(SlepcFinalize()); |
113 | return 0; | ||
114 | } | ||
115 | |||
116 | /* | ||
117 | Matrix generator for a Markov model of a random walk on a triangular grid. | ||
118 | See ex5.c for additional details. | ||
119 | */ | ||
120 | 10 | PetscErrorCode MatMarkovModel(PetscInt m,Mat A) | |
121 | { | ||
122 | 10 | const PetscReal cst = 0.5/(PetscReal)(m-1); | |
123 | 10 | PetscReal pd,pu; | |
124 | 10 | PetscInt Istart,Iend,i,j,jmax,ix=0; | |
125 | |||
126 |
1/2✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
|
10 | PetscFunctionBeginUser; |
127 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatGetOwnershipRange(A,&Istart,&Iend)); |
128 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
160 | for (i=1;i<=m;i++) { |
129 | 150 | jmax = m-i+1; | |
130 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1350 | for (j=1;j<=jmax;j++) { |
131 | 1200 | ix = ix + 1; | |
132 |
2/4✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
|
1200 | if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */ |
133 |
2/2✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
|
1200 | if (j!=jmax) { |
134 | 1050 | pd = cst*(PetscReal)(i+j-1); | |
135 | /* north */ | ||
136 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1050 | if (i==1) PetscCall(MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES)); |
137 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
910 | else PetscCall(MatSetValue(A,ix-1,ix,pd,INSERT_VALUES)); |
138 | /* east */ | ||
139 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1050 | if (j==1) PetscCall(MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES)); |
140 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
910 | else PetscCall(MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES)); |
141 | } | ||
142 | /* south */ | ||
143 | 1200 | pu = 0.5 - cst*(PetscReal)(i+j-3); | |
144 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1200 | if (j>1) PetscCall(MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES)); |
145 | /* west */ | ||
146 |
6/8✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
|
1200 | if (i>1) PetscCall(MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES)); |
147 | } | ||
148 | } | ||
149 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); |
150 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
|
10 | PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); |
151 |
5/12✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
2 | PetscFunctionReturn(PETSC_SUCCESS); |
152 | } | ||
153 | |||
154 | /*TEST | ||
155 | |||
156 | test: | ||
157 | suffix: 1 | ||
158 | args: -mfn_ncv 6 | ||
159 | |||
160 | TEST*/ | ||
161 |