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: "krylov"
      12             : 
      13             :    Method: Arnoldi with Eiermann-Ernst restart
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Build Arnoldi approximations using f(H) for the Hessenberg matrix H,
      18             :        restart by discarding the Krylov basis but keeping H.
      19             : 
      20             :    References:
      21             : 
      22             :        [1] M. Eiermann and O. Ernst, "A restarted Krylov subspace method
      23             :            for the evaluation of matrix functions", SIAM J. Numer. Anal.
      24             :            44(6):2481-2504, 2006.
      25             : */
      26             : 
      27             : #include <slepc/private/mfnimpl.h>
      28             : #include <slepcblaslapack.h>
      29             : 
      30          14 : static PetscErrorCode MFNSetUp_Krylov(MFN mfn)
      31             : {
      32          14 :   PetscInt       N;
      33             : 
      34          14 :   PetscFunctionBegin;
      35          14 :   PetscCall(MatGetSize(mfn->A,&N,NULL));
      36          14 :   if (mfn->ncv==PETSC_DETERMINE) mfn->ncv = PetscMin(30,N);
      37          14 :   if (mfn->max_it==PETSC_DETERMINE) mfn->max_it = 100;
      38          14 :   PetscCall(MFNAllocateSolution(mfn,1));
      39          14 :   PetscFunctionReturn(PETSC_SUCCESS);
      40             : }
      41             : 
      42          95 : static PetscErrorCode MFNSolve_Krylov(MFN mfn,Vec b,Vec x)
      43             : {
      44          95 :   PetscInt          n=0,m,ld,ldh,j;
      45          95 :   PetscBLASInt      m_,inc=1;
      46          95 :   Mat               M,G=NULL,H=NULL;
      47          95 :   Vec               F=NULL;
      48          95 :   PetscScalar       *marray,*farray,*harray;
      49          95 :   const PetscScalar *garray;
      50          95 :   PetscReal         beta,betaold=0.0,nrm=1.0;
      51          95 :   PetscBool         breakdown;
      52             : 
      53          95 :   PetscFunctionBegin;
      54          95 :   m  = mfn->ncv;
      55          95 :   ld = m+1;
      56          95 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,m,NULL,&M));
      57          95 :   PetscCall(MatDenseGetArray(M,&marray));
      58             : 
      59             :   /* set initial vector to b/||b|| */
      60          95 :   PetscCall(BVInsertVec(mfn->V,0,b));
      61          95 :   PetscCall(BVScaleColumn(mfn->V,0,1.0/mfn->bnorm));
      62          95 :   PetscCall(VecSet(x,0.0));
      63             : 
      64             :   /* Restart loop */
      65         320 :   while (mfn->reason == MFN_CONVERGED_ITERATING) {
      66         225 :     mfn->its++;
      67             : 
      68             :     /* compute Arnoldi factorization */
      69         225 :     PetscCall(BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,M,0,&m,&beta,&breakdown));
      70             : 
      71             :     /* save previous Hessenberg matrix in G; allocate new storage for H and f(H) */
      72         225 :     if (mfn->its>1) { G = H; H = NULL; }
      73         225 :     ldh = n+m;
      74         225 :     PetscCall(MFN_CreateVec(ldh,&F));
      75         225 :     PetscCall(MFN_CreateDenseMat(ldh,&H));
      76             : 
      77             :     /* glue together the previous H and the new H obtained with Arnoldi */
      78         225 :     PetscCall(MatDenseGetArray(H,&harray));
      79        5087 :     for (j=0;j<m;j++) PetscCall(PetscArraycpy(harray+n+(j+n)*ldh,marray+j*ld,m));
      80         225 :     if (mfn->its>1) {
      81         130 :       PetscCall(MatDenseGetArrayRead(G,&garray));
      82        4388 :       for (j=0;j<n;j++) PetscCall(PetscArraycpy(harray+j*ldh,garray+j*n,n));
      83         130 :       PetscCall(MatDenseRestoreArrayRead(G,&garray));
      84         130 :       PetscCall(MatDestroy(&G));
      85         130 :       harray[n+(n-1)*ldh] = betaold;
      86             :     }
      87         225 :     PetscCall(MatDenseRestoreArray(H,&harray));
      88             : 
      89         225 :     if (mfn->its==1) {
      90             :       /* set symmetry flag of H from A */
      91          95 :       PetscCall(MatPropagateSymmetryOptions(mfn->A,H));
      92             :     }
      93             : 
      94             :     /* evaluate f(H) */
      95         225 :     PetscCall(FNEvaluateFunctionMatVec(mfn->fn,H,F));
      96             : 
      97             :     /* x += ||b||*V*f(H)*e_1 */
      98         225 :     PetscCall(VecGetArray(F,&farray));
      99         225 :     PetscCall(PetscBLASIntCast(m,&m_));
     100         225 :     nrm = BLASnrm2_(&m_,farray+n,&inc);   /* relative norm of the update ||u||/||b|| */
     101         225 :     PetscCall(MFNMonitor(mfn,mfn->its,nrm));
     102        5087 :     for (j=0;j<m;j++) farray[j+n] *= mfn->bnorm;
     103         225 :     PetscCall(BVSetActiveColumns(mfn->V,0,m));
     104         225 :     PetscCall(BVMultVec(mfn->V,1.0,1.0,x,farray+n));
     105         225 :     PetscCall(VecRestoreArray(F,&farray));
     106             : 
     107             :     /* check convergence */
     108         225 :     if (mfn->its >= mfn->max_it) mfn->reason = MFN_DIVERGED_ITS;
     109         225 :     if (mfn->its>1) {
     110         130 :       if (m<mfn->ncv || breakdown || beta==0.0 || nrm<mfn->tol) mfn->reason = MFN_CONVERGED_TOL;
     111             :     }
     112             : 
     113             :     /* restart with vector v_{m+1} */
     114         225 :     if (mfn->reason == MFN_CONVERGED_ITERATING) {
     115         130 :       PetscCall(BVCopyColumn(mfn->V,m,0));
     116         130 :       n += m;
     117         130 :       betaold = beta;
     118             :     }
     119             :   }
     120             : 
     121          95 :   PetscCall(MatDestroy(&H));
     122          95 :   PetscCall(MatDestroy(&G));
     123          95 :   PetscCall(VecDestroy(&F));
     124          95 :   PetscCall(MatDenseRestoreArray(M,&marray));
     125          95 :   PetscCall(MatDestroy(&M));
     126          95 :   PetscFunctionReturn(PETSC_SUCCESS);
     127             : }
     128             : 
     129          14 : SLEPC_EXTERN PetscErrorCode MFNCreate_Krylov(MFN mfn)
     130             : {
     131          14 :   PetscFunctionBegin;
     132          14 :   mfn->ops->solve          = MFNSolve_Krylov;
     133          14 :   mfn->ops->setup          = MFNSetUp_Krylov;
     134          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     135             : }