Line data    Source code

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