Actual source code: mfnkrylov.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: PetscErrorCode MFNSetUp_Krylov(MFN mfn)
31: {
33: PetscInt N;
36: MatGetSize(mfn->A,&N,NULL);
37: if (!mfn->ncv) mfn->ncv = PetscMin(30,N);
38: if (!mfn->max_it) mfn->max_it = 100;
39: MFNAllocateSolution(mfn,1);
40: return(0);
41: }
43: PetscErrorCode MFNBasicArnoldi(MFN mfn,PetscScalar *H,PetscInt ldh,PetscInt k,PetscInt *M,PetscReal *beta,PetscBool *breakdown)
44: {
46: PetscScalar *a;
47: PetscInt j,nc,n,m = *M;
48: Vec vj,vj1,buf;
51: BVSetActiveColumns(mfn->V,0,m);
52: for (j=k;j<m;j++) {
53: BVGetColumn(mfn->V,j,&vj);
54: BVGetColumn(mfn->V,j+1,&vj1);
55: MatMult(mfn->A,vj,vj1);
56: BVRestoreColumn(mfn->V,j,&vj);
57: BVRestoreColumn(mfn->V,j+1,&vj1);
58: BVOrthonormalizeColumn(mfn->V,j+1,PETSC_FALSE,beta,breakdown);
59: if (*breakdown) {
60: *M = j+1;
61: break;
62: }
63: }
64: /* extract Hessenberg matrix from the BV object */
65: BVGetNumConstraints(mfn->V,&nc);
66: BVGetSizes(mfn->V,NULL,NULL,&n);
67: BVGetBufferVec(mfn->V,&buf);
68: VecGetArray(buf,&a);
69: for (j=k;j<*M;j++) {
70: PetscMemcpy(H+j*ldh,a+nc+(j+1)*(nc+n),(j+2)*sizeof(PetscScalar));
71: }
72: VecRestoreArray(buf,&a);
73: return(0);
74: }
76: PetscErrorCode MFNSolve_Krylov(MFN mfn,Vec b,Vec x)
77: {
79: PetscInt n=0,m,ld,ldh,j;
80: PetscBLASInt m_,inc=1;
81: Mat G=NULL,H=NULL;
82: Vec F=NULL;
83: PetscScalar *array,*farray,*garray,*harray;
84: PetscReal beta,betaold=0.0,nrm=1.0;
85: PetscBool breakdown,set,flg,symm=PETSC_FALSE;
88: m = mfn->ncv;
89: ld = m+1;
90: PetscCalloc1(ld*ld,&array);
92: /* set initial vector to b/||b|| */
93: BVInsertVec(mfn->V,0,b);
94: BVScaleColumn(mfn->V,0,1.0/mfn->bnorm);
95: VecSet(x,0.0);
97: /* Restart loop */
98: while (mfn->reason == MFN_CONVERGED_ITERATING) {
99: mfn->its++;
101: /* compute Arnoldi factorization */
102: MFNBasicArnoldi(mfn,array,ld,0,&m,&beta,&breakdown);
104: /* save previous Hessenberg matrix in G; allocate new storage for H and f(H) */
105: if (mfn->its>1) { G = H; H = NULL; }
106: ldh = n+m;
107: MFN_CreateVec(ldh,&F);
108: MFN_CreateDenseMat(ldh,&H);
110: /* glue together the previous H and the new H obtained with Arnoldi */
111: MatDenseGetArray(H,&harray);
112: for (j=0;j<m;j++) {
113: PetscMemcpy(harray+n+(j+n)*ldh,array+j*ld,m*sizeof(PetscScalar));
114: }
115: if (mfn->its>1) {
116: MatDenseGetArray(G,&garray);
117: for (j=0;j<n;j++) {
118: PetscMemcpy(harray+j*ldh,garray+j*n,n*sizeof(PetscScalar));
119: }
120: MatDenseRestoreArray(G,&garray);
121: MatDestroy(&G);
122: harray[n+(n-1)*ldh] = betaold;
123: }
124: MatDenseRestoreArray(H,&harray);
126: if (mfn->its==1) {
127: /* set symmetry flag of H from A */
128: MatIsHermitianKnown(mfn->A,&set,&flg);
129: symm = set? flg: PETSC_FALSE;
130: if (symm) {
131: MatSetOption(H,MAT_HERMITIAN,PETSC_TRUE);
132: }
133: }
135: /* evaluate f(H) */
136: FNEvaluateFunctionMatVec(mfn->fn,H,F);
138: /* x += ||b||*V*f(H)*e_1 */
139: VecGetArray(F,&farray);
140: PetscBLASIntCast(m,&m_);
141: nrm = BLASnrm2_(&m_,farray+n,&inc); /* relative norm of the update ||u||/||b|| */
142: MFNMonitor(mfn,mfn->its,nrm);
143: for (j=0;j<m;j++) farray[j+n] *= mfn->bnorm;
144: BVSetActiveColumns(mfn->V,0,m);
145: BVMultVec(mfn->V,1.0,1.0,x,farray+n);
146: VecRestoreArray(F,&farray);
148: /* check convergence */
149: if (mfn->its >= mfn->max_it) mfn->reason = MFN_DIVERGED_ITS;
150: if (mfn->its>1) {
151: if (m<mfn->ncv || breakdown || beta==0.0 || nrm<mfn->tol) mfn->reason = MFN_CONVERGED_TOL;
152: }
154: /* restart with vector v_{m+1} */
155: if (mfn->reason == MFN_CONVERGED_ITERATING) {
156: BVCopyColumn(mfn->V,m,0);
157: n += m;
158: betaold = beta;
159: }
160: }
162: MatDestroy(&H);
163: MatDestroy(&G);
164: VecDestroy(&F);
165: PetscFree(array);
166: return(0);
167: }
169: SLEPC_EXTERN PetscErrorCode MFNCreate_Krylov(MFN mfn)
170: {
172: mfn->ops->solve = MFNSolve_Krylov;
173: mfn->ops->setup = MFNSetUp_Krylov;
174: return(0);
175: }