Actual source code: mfnexpokit.c
slepc-3.18.3 2023-03-24
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: */
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: {
30: PetscInt N;
31: PetscBool isexp;
33: MatGetSize(mfn->A,&N,NULL);
34: if (mfn->ncv==PETSC_DEFAULT) mfn->ncv = PetscMin(30,N);
35: if (mfn->max_it==PETSC_DEFAULT) mfn->max_it = 100;
36: MFNAllocateSolution(mfn,2);
38: PetscObjectTypeCompare((PetscObject)mfn->fn,FNEXP,&isexp);
40: return 0;
41: }
43: PetscErrorCode MFNSolve_Expokit(MFN mfn,Vec b,Vec x)
44: {
45: PetscInt mxstep,mxrej,m,mb,ld,i,j,ireject,mx,k1;
46: Vec v,r;
47: Mat H,M=NULL,K=NULL;
48: FN fn;
49: PetscScalar *Harray,*B,*F,*betaF,t,sgn,sfactor;
50: const PetscScalar *pK;
51: PetscReal anorm,avnorm,tol,err_loc,rndoff,t_out,t_new,t_now,t_step;
52: PetscReal xm,fact,s,p1,p2,beta,beta2,gamma,delta;
53: PetscBool breakdown;
55: m = mfn->ncv;
56: tol = mfn->tol;
57: mxstep = mfn->max_it;
58: mxrej = 10;
59: gamma = 0.9;
60: delta = 1.2;
61: mb = m;
62: FNGetScale(mfn->fn,&t,&sfactor);
63: FNDuplicate(mfn->fn,PetscObjectComm((PetscObject)mfn->fn),&fn);
64: FNSetScale(fn,1.0,1.0);
65: t_out = PetscAbsScalar(t);
66: t_now = 0.0;
67: MatNorm(mfn->A,NORM_INFINITY,&anorm);
68: rndoff = anorm*PETSC_MACHINE_EPSILON;
70: k1 = 2;
71: xm = 1.0/(PetscReal)m;
72: beta = mfn->bnorm;
73: fact = PetscPowRealInt((m+1)/2.72,m+1)*PetscSqrtReal(2.0*PETSC_PI*(m+1));
74: t_new = (1.0/anorm)*PetscPowReal((fact*tol)/(4.0*beta*anorm),xm);
75: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
76: t_new = PetscCeilReal(t_new/s)*s;
77: sgn = t/PetscAbsScalar(t);
79: VecCopy(b,x);
80: ld = m+2;
81: PetscCalloc2(m+1,&betaF,ld*ld,&B);
82: MatCreateSeqDense(PETSC_COMM_SELF,ld,ld,NULL,&H);
83: MatDenseGetArray(H,&Harray);
85: while (mfn->reason == MFN_CONVERGED_ITERATING) {
86: mfn->its++;
87: if (PetscIsInfOrNanReal(t_new)) t_new = PETSC_MAX_REAL;
88: t_step = PetscMin(t_out-t_now,t_new);
89: BVInsertVec(mfn->V,0,x);
90: BVScaleColumn(mfn->V,0,1.0/beta);
91: BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,H,0,&mb,&beta2,&breakdown);
92: if (breakdown) {
93: k1 = 0;
94: t_step = t_out-t_now;
95: }
96: if (k1!=0) {
97: Harray[m+1+ld*m] = 1.0;
98: BVGetColumn(mfn->V,m,&v);
99: BVGetColumn(mfn->V,m+1,&r);
100: MatMult(mfn->transpose_solve?mfn->AT:mfn->A,v,r);
101: BVRestoreColumn(mfn->V,m,&v);
102: BVRestoreColumn(mfn->V,m+1,&r);
103: BVNormColumn(mfn->V,m+1,NORM_2,&avnorm);
104: }
105: PetscArraycpy(B,Harray,ld*ld);
107: ireject = 0;
108: while (ireject <= mxrej) {
109: mx = mb + k1;
110: for (i=0;i<mx;i++) {
111: for (j=0;j<mx;j++) {
112: Harray[i+j*ld] = sgn*B[i+j*ld]*t_step;
113: }
114: }
115: MFN_CreateDenseMat(mx,&M);
116: MFN_CreateDenseMat(mx,&K);
117: MatDenseGetArray(M,&F);
118: for (j=0;j<mx;j++) PetscArraycpy(F+j*mx,Harray+j*ld,mx);
119: MatDenseRestoreArray(M,&F);
120: FNEvaluateFunctionMat(fn,M,K);
122: if (k1==0) {
123: err_loc = tol;
124: break;
125: } else {
126: MatDenseGetArrayRead(K,&pK);
127: p1 = PetscAbsScalar(beta*pK[m]);
128: p2 = PetscAbsScalar(beta*pK[m+1]*avnorm);
129: MatDenseRestoreArrayRead(K,&pK);
130: if (p1 > 10*p2) {
131: err_loc = p2;
132: xm = 1.0/(PetscReal)m;
133: } else if (p1 > p2) {
134: err_loc = (p1*p2)/(p1-p2);
135: xm = 1.0/(PetscReal)m;
136: } else {
137: err_loc = p1;
138: xm = 1.0/(PetscReal)(m-1);
139: }
140: }
141: if (err_loc <= delta*t_step*tol) break;
142: else {
143: t_step = gamma*t_step*PetscPowReal(t_step*tol/err_loc,xm);
144: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_step))-1);
145: t_step = PetscCeilReal(t_step/s)*s;
146: ireject = ireject+1;
147: }
148: }
150: mx = mb + PetscMax(0,k1-1);
151: MatDenseGetArrayRead(K,&pK);
152: for (j=0;j<mx;j++) betaF[j] = beta*pK[j];
153: MatDenseRestoreArrayRead(K,&pK);
154: BVSetActiveColumns(mfn->V,0,mx);
155: BVMultVec(mfn->V,1.0,0.0,x,betaF);
156: VecNorm(x,NORM_2,&beta);
158: t_now = t_now+t_step;
159: if (t_now>=t_out) mfn->reason = MFN_CONVERGED_TOL;
160: else {
161: t_new = gamma*t_step*PetscPowReal((t_step*tol)/err_loc,xm);
162: s = PetscPowReal(10.0,PetscFloorReal(PetscLog10Real(t_new))-1);
163: t_new = PetscCeilReal(t_new/s)*s;
164: }
165: err_loc = PetscMax(err_loc,rndoff);
166: if (mfn->its==mxstep) mfn->reason = MFN_DIVERGED_ITS;
167: MFNMonitor(mfn,mfn->its,err_loc);
168: }
169: VecScale(x,sfactor);
171: MatDestroy(&M);
172: MatDestroy(&K);
173: FNDestroy(&fn);
174: MatDenseRestoreArray(H,&Harray);
175: MatDestroy(&H);
176: PetscFree2(betaF,B);
177: return 0;
178: }
180: SLEPC_EXTERN PetscErrorCode MFNCreate_Expokit(MFN mfn)
181: {
182: mfn->ops->solve = MFNSolve_Expokit;
183: mfn->ops->setup = MFNSetUp_Expokit;
184: return 0;
185: }