Actual source code: mfnexpokit.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: */
10: /*
11: SLEPc matrix function solver: "expokit"
13: Method: Arnoldi method tailored for the matrix exponential
15: Algorithm:
17: Uses Arnoldi relations to compute exp(t_step*A)*v_last for
18: several time steps.
20: References:
22: [1] R. Sidje, "Expokit: a software package for computing matrix
23: exponentials", ACM Trans. Math. Softw. 24(1):130-156, 1998.
24: */
26: #include <slepc/private/mfnimpl.h>
28: PetscErrorCode MFNSetUp_Expokit(MFN mfn)
29: {
31: PetscInt N;
32: PetscBool isexp;
35: MatGetSize(mfn->A,&N,NULL);
36: if (!mfn->ncv) mfn->ncv = PetscMin(30,N);
37: if (!mfn->max_it) mfn->max_it = 100;
38: MFNAllocateSolution(mfn,2);
40: PetscObjectTypeCompare((PetscObject)mfn->fn,FNEXP,&isexp);
41: if (!isexp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This solver only supports the exponential function");
42: return(0);
43: }
45: PetscErrorCode MFNSolve_Expokit(MFN mfn,Vec b,Vec x)
46: {
48: PetscInt mxstep,mxrej,m,mb,ld,i,j,ireject,mx,k1;
49: Vec v,r;
50: Mat M=NULL,K=NULL;
51: FN fn;
52: PetscScalar *H,*B,*F,*betaF,t,sgn,sfactor;
53: PetscReal anorm,avnorm,tol,err_loc,rndoff;
54: PetscReal t_out,t_new,t_now,t_step;
55: PetscReal xm,fact,s,p1,p2;
56: PetscReal beta,beta2,gamma,delta;
57: PetscBool breakdown;
60: m = mfn->ncv;
61: tol = mfn->tol;
62: mxstep = mfn->max_it;
63: mxrej = 10;
64: gamma = 0.9;
65: delta = 1.2;
66: mb = m;
67: FNGetScale(mfn->fn,&t,&sfactor);
68: FNDuplicate(mfn->fn,PetscObjectComm((PetscObject)mfn->fn),&fn);
69: FNSetScale(fn,1.0,1.0);
70: t_out = PetscAbsScalar(t);
71: t_now = 0.0;
72: MatNorm(mfn->A,NORM_INFINITY,&anorm);
73: rndoff = anorm*PETSC_MACHINE_EPSILON;
75: k1 = 2;
76: xm = 1.0/(PetscReal)m;
77: beta = mfn->bnorm;
78: fact = PetscPowRealInt((m+1)/2.72,m+1)*PetscSqrtReal(2*PETSC_PI*(m+1));
79: t_new = (1.0/anorm)*PetscPowReal((fact*tol)/(4.0*beta*anorm),xm);
80: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
81: t_new = PetscCeilReal(t_new/s)*s;
82: sgn = t/PetscAbsScalar(t);
84: VecCopy(b,x);
85: ld = m+2;
86: PetscCalloc3(m+1,&betaF,ld*ld,&H,ld*ld,&B);
88: while (mfn->reason == MFN_CONVERGED_ITERATING) {
89: mfn->its++;
90: if (PetscIsInfOrNanReal(t_new)) t_new = PETSC_MAX_REAL;
91: t_step = PetscMin(t_out-t_now,t_new);
92: BVInsertVec(mfn->V,0,x);
93: BVScaleColumn(mfn->V,0,1.0/beta);
94: MFNBasicArnoldi(mfn,H,ld,0,&mb,&beta2,&breakdown);
95: if (breakdown) {
96: k1 = 0;
97: t_step = t_out-t_now;
98: }
99: if (k1!=0) {
100: H[m+1+ld*m] = 1.0;
101: BVGetColumn(mfn->V,m,&v);
102: BVGetColumn(mfn->V,m+1,&r);
103: MatMult(mfn->A,v,r);
104: BVRestoreColumn(mfn->V,m,&v);
105: BVRestoreColumn(mfn->V,m+1,&r);
106: BVNormColumn(mfn->V,m+1,NORM_2,&avnorm);
107: }
108: PetscMemcpy(B,H,ld*ld*sizeof(PetscScalar));
110: ireject = 0;
111: while (ireject <= mxrej) {
112: mx = mb + k1;
113: for (i=0;i<mx;i++) {
114: for (j=0;j<mx;j++) {
115: H[i+j*ld] = sgn*B[i+j*ld]*t_step;
116: }
117: }
118: MFN_CreateDenseMat(mx,&M);
119: MFN_CreateDenseMat(mx,&K);
120: MatDenseGetArray(M,&F);
121: for (j=0;j<mx;j++) {
122: PetscMemcpy(F+j*mx,H+j*ld,mx*sizeof(PetscScalar));
123: }
124: MatDenseRestoreArray(M,&F);
125: FNEvaluateFunctionMat(fn,M,K);
127: if (k1==0) {
128: err_loc = tol;
129: break;
130: } else {
131: MatDenseGetArray(K,&F);
132: p1 = PetscAbsScalar(beta*F[m]);
133: p2 = PetscAbsScalar(beta*F[m+1]*avnorm);
134: MatDenseRestoreArray(K,&F);
135: if (p1 > 10*p2) {
136: err_loc = p2;
137: xm = 1.0/(PetscReal)m;
138: } else if (p1 > p2) {
139: err_loc = (p1*p2)/(p1-p2);
140: xm = 1.0/(PetscReal)m;
141: } else {
142: err_loc = p1;
143: xm = 1.0/(PetscReal)(m-1);
144: }
145: }
146: if (err_loc <= delta*t_step*tol) break;
147: else {
148: t_step = gamma*t_step*PetscPowReal(t_step*tol/err_loc,xm);
149: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_step))-1);
150: t_step = PetscCeilReal(t_step/s)*s;
151: ireject = ireject+1;
152: }
153: }
155: mx = mb + PetscMax(0,k1-1);
156: MatDenseGetArray(K,&F);
157: for (j=0;j<mx;j++) betaF[j] = beta*F[j];
158: MatDenseRestoreArray(K,&F);
159: BVSetActiveColumns(mfn->V,0,mx);
160: BVMultVec(mfn->V,1.0,0.0,x,betaF);
161: VecNorm(x,NORM_2,&beta);
163: t_now = t_now+t_step;
164: if (t_now>=t_out) mfn->reason = MFN_CONVERGED_TOL;
165: else {
166: t_new = gamma*t_step*PetscPowReal((t_step*tol)/err_loc,xm);
167: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
168: t_new = PetscCeilReal(t_new/s)*s;
169: }
170: err_loc = PetscMax(err_loc,rndoff);
171: if (mfn->its==mxstep) mfn->reason = MFN_DIVERGED_ITS;
172: MFNMonitor(mfn,mfn->its,err_loc);
173: }
174: VecScale(x,sfactor);
176: MatDestroy(&M);
177: MatDestroy(&K);
178: FNDestroy(&fn);
179: PetscFree3(betaF,H,B);
180: return(0);
181: }
183: SLEPC_EXTERN PetscErrorCode MFNCreate_Expokit(MFN mfn)
184: {
186: mfn->ops->solve = MFNSolve_Expokit;
187: mfn->ops->setup = MFNSetUp_Expokit;
188: return(0);
189: }