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