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: "krylov"

 13:    Method: Arnoldi with Eiermann-Ernst restart

 15:    Algorithm:

 17:        Build Arnoldi approximations using f(H) for the Hessenberg matrix H,
 18:        restart by discarding the Krylov basis but keeping H.

 20:    References:

 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: */

 27: #include <slepc/private/mfnimpl.h>
 28: #include <slepcblaslapack.h>

 30: static PetscErrorCode MFNSetUp_Krylov(MFN mfn)
 31: {
 32:   PetscInt       N;

 34:   PetscFunctionBegin;
 35:   PetscCall(MatGetSize(mfn->A,&N,NULL));
 36:   if (mfn->ncv==PETSC_DEFAULT) mfn->ncv = PetscMin(30,N);
 37:   if (mfn->max_it==PETSC_DEFAULT) mfn->max_it = 100;
 38:   PetscCall(MFNAllocateSolution(mfn,1));
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: static PetscErrorCode MFNSolve_Krylov(MFN mfn,Vec b,Vec x)
 43: {
 44:   PetscInt          n=0,m,ld,ldh,j;
 45:   PetscBLASInt      m_,inc=1;
 46:   Mat               M,G=NULL,H=NULL;
 47:   Vec               F=NULL;
 48:   PetscScalar       *marray,*farray,*harray;
 49:   const PetscScalar *garray;
 50:   PetscReal         beta,betaold=0.0,nrm=1.0;
 51:   PetscBool         breakdown;

 53:   PetscFunctionBegin;
 54:   m  = mfn->ncv;
 55:   ld = m+1;
 56:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ld,m,NULL,&M));
 57:   PetscCall(MatDenseGetArray(M,&marray));

 59:   /* set initial vector to b/||b|| */
 60:   PetscCall(BVInsertVec(mfn->V,0,b));
 61:   PetscCall(BVScaleColumn(mfn->V,0,1.0/mfn->bnorm));
 62:   PetscCall(VecSet(x,0.0));

 64:   /* Restart loop */
 65:   while (mfn->reason == MFN_CONVERGED_ITERATING) {
 66:     mfn->its++;

 68:     /* compute Arnoldi factorization */
 69:     PetscCall(BVMatArnoldi(mfn->V,mfn->transpose_solve?mfn->AT:mfn->A,M,0,&m,&beta,&breakdown));

 71:     /* save previous Hessenberg matrix in G; allocate new storage for H and f(H) */
 72:     if (mfn->its>1) { G = H; H = NULL; }
 73:     ldh = n+m;
 74:     PetscCall(MFN_CreateVec(ldh,&F));
 75:     PetscCall(MFN_CreateDenseMat(ldh,&H));

 77:     /* glue together the previous H and the new H obtained with Arnoldi */
 78:     PetscCall(MatDenseGetArray(H,&harray));
 79:     for (j=0;j<m;j++) PetscCall(PetscArraycpy(harray+n+(j+n)*ldh,marray+j*ld,m));
 80:     if (mfn->its>1) {
 81:       PetscCall(MatDenseGetArrayRead(G,&garray));
 82:       for (j=0;j<n;j++) PetscCall(PetscArraycpy(harray+j*ldh,garray+j*n,n));
 83:       PetscCall(MatDenseRestoreArrayRead(G,&garray));
 84:       PetscCall(MatDestroy(&G));
 85:       harray[n+(n-1)*ldh] = betaold;
 86:     }
 87:     PetscCall(MatDenseRestoreArray(H,&harray));

 89:     if (mfn->its==1) {
 90:       /* set symmetry flag of H from A */
 91:       PetscCall(MatPropagateSymmetryOptions(mfn->A,H));
 92:     }

 94:     /* evaluate f(H) */
 95:     PetscCall(FNEvaluateFunctionMatVec(mfn->fn,H,F));

 97:     /* x += ||b||*V*f(H)*e_1 */
 98:     PetscCall(VecGetArray(F,&farray));
 99:     PetscCall(PetscBLASIntCast(m,&m_));
100:     nrm = BLASnrm2_(&m_,farray+n,&inc);   /* relative norm of the update ||u||/||b|| */
101:     PetscCall(MFNMonitor(mfn,mfn->its,nrm));
102:     for (j=0;j<m;j++) farray[j+n] *= mfn->bnorm;
103:     PetscCall(BVSetActiveColumns(mfn->V,0,m));
104:     PetscCall(BVMultVec(mfn->V,1.0,1.0,x,farray+n));
105:     PetscCall(VecRestoreArray(F,&farray));

107:     /* check convergence */
108:     if (mfn->its >= mfn->max_it) mfn->reason = MFN_DIVERGED_ITS;
109:     if (mfn->its>1) {
110:       if (m<mfn->ncv || breakdown || beta==0.0 || nrm<mfn->tol) mfn->reason = MFN_CONVERGED_TOL;
111:     }

113:     /* restart with vector v_{m+1} */
114:     if (mfn->reason == MFN_CONVERGED_ITERATING) {
115:       PetscCall(BVCopyColumn(mfn->V,m,0));
116:       n += m;
117:       betaold = beta;
118:     }
119:   }

121:   PetscCall(MatDestroy(&H));
122:   PetscCall(MatDestroy(&G));
123:   PetscCall(VecDestroy(&F));
124:   PetscCall(MatDenseRestoreArray(M,&marray));
125:   PetscCall(MatDestroy(&M));
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: SLEPC_EXTERN PetscErrorCode MFNCreate_Krylov(MFN mfn)
130: {
131:   PetscFunctionBegin;
132:   mfn->ops->solve          = MFNSolve_Krylov;
133:   mfn->ops->setup          = MFNSetUp_Krylov;
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }