Actual source code: ex23.c
slepc-3.9.2 2018-07-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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[] = "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\n"
14: " -m <m>, where <m> = number of grid subdivisions in each dimension.\n\n";
16: #include <slepcmfn.h>
18: /*
19: User-defined routines
20: */
21: PetscErrorCode MatMarkovModel(PetscInt m,Mat A);
23: int main(int argc,char **argv)
24: {
25: Mat A; /* problem matrix */
26: MFN mfn;
27: FN f;
28: PetscReal tol,norm;
29: PetscScalar t=2.0;
30: Vec v,y;
31: PetscInt N,m=15,ncv,maxit,its;
32: PetscErrorCode ierr;
33: PetscBool draw_sol;
34: MFNConvergedReason reason;
36: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
38: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
39: PetscOptionsGetScalar(NULL,NULL,"-t",&t,NULL);
40: N = m*(m+1)/2;
41: PetscPrintf(PETSC_COMM_WORLD,"\nMarkov y=exp(t*A)*e_1, N=%D (m=%D)\n\n",N,m);
43: PetscOptionsHasName(NULL,NULL,"-draw_sol",&draw_sol);
45: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: Compute the transition probability matrix, A
47: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49: MatCreate(PETSC_COMM_WORLD,&A);
50: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
51: MatSetFromOptions(A);
52: MatSetUp(A);
53: MatMarkovModel(m,A);
55: /* set v = e_1 */
56: MatCreateVecs(A,NULL,&y);
57: MatCreateVecs(A,NULL,&v);
58: VecSetValue(v,0,1.0,INSERT_VALUES);
59: VecAssemblyBegin(v);
60: VecAssemblyEnd(v);
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create the solver and set various options
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: /*
66: Create matrix function solver context
67: */
68: MFNCreate(PETSC_COMM_WORLD,&mfn);
70: /*
71: Set operator matrix, the function to compute, and other options
72: */
73: MFNSetOperator(mfn,A);
74: MFNGetFN(mfn,&f);
75: FNSetType(f,FNEXP);
76: FNSetScale(f,t,1.0);
77: MFNSetTolerances(mfn,1e-07,PETSC_DEFAULT);
79: /*
80: Set solver parameters at runtime
81: */
82: MFNSetFromOptions(mfn);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Solve the problem, y=exp(t*A)*v
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: MFNSolve(mfn,v,y);
89: MFNGetConvergedReason(mfn,&reason);
90: if (reason<0) SETERRQ(PETSC_COMM_WORLD,1,"Solver did not converge");
91: VecNorm(y,NORM_2,&norm);
93: /*
94: Optional: Get some information from the solver and display it
95: */
96: MFNGetIterationNumber(mfn,&its);
97: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
98: MFNGetDimensions(mfn,&ncv);
99: PetscPrintf(PETSC_COMM_WORLD," Subspace dimension: %D\n",ncv);
100: MFNGetTolerances(mfn,&tol,&maxit);
101: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Display solution and clean up
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscPrintf(PETSC_COMM_WORLD," Computed vector at time t=%.4g has norm %g\n\n",(double)PetscRealPart(t),(double)norm);
107: if (draw_sol) {
108: PetscViewerDrawSetPause(PETSC_VIEWER_DRAW_WORLD,-1);
109: VecView(y,PETSC_VIEWER_DRAW_WORLD);
110: }
112: /*
113: Free work space
114: */
115: MFNDestroy(&mfn);
116: MatDestroy(&A);
117: VecDestroy(&v);
118: VecDestroy(&y);
119: SlepcFinalize();
120: return ierr;
121: }
123: /*
124: Matrix generator for a Markov model of a random walk on a triangular grid.
125: See ex5.c for additional details.
126: */
127: PetscErrorCode MatMarkovModel(PetscInt m,Mat A)
128: {
129: const PetscReal cst = 0.5/(PetscReal)(m-1);
130: PetscReal pd,pu;
131: PetscInt Istart,Iend,i,j,jmax,ix=0;
132: PetscErrorCode ierr;
135: MatGetOwnershipRange(A,&Istart,&Iend);
136: for (i=1;i<=m;i++) {
137: jmax = m-i+1;
138: for (j=1;j<=jmax;j++) {
139: ix = ix + 1;
140: if (ix-1<Istart || ix>Iend) continue; /* compute only owned rows */
141: if (j!=jmax) {
142: pd = cst*(PetscReal)(i+j-1);
143: /* north */
144: if (i==1) {
145: MatSetValue(A,ix-1,ix,2*pd,INSERT_VALUES);
146: } else {
147: MatSetValue(A,ix-1,ix,pd,INSERT_VALUES);
148: }
149: /* east */
150: if (j==1) {
151: MatSetValue(A,ix-1,ix+jmax-1,2*pd,INSERT_VALUES);
152: } else {
153: MatSetValue(A,ix-1,ix+jmax-1,pd,INSERT_VALUES);
154: }
155: }
156: /* south */
157: pu = 0.5 - cst*(PetscReal)(i+j-3);
158: if (j>1) {
159: MatSetValue(A,ix-1,ix-2,pu,INSERT_VALUES);
160: }
161: /* west */
162: if (i>1) {
163: MatSetValue(A,ix-1,ix-jmax-2,pu,INSERT_VALUES);
164: }
165: }
166: }
167: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
169: return(0);
170: }