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: MFN routines related to problem setup
12: */
14: #include <slepc/private/mfnimpl.h> /*I "slepcmfn.h" I*/
16: /*@
17: MFNSetUp - Sets up all the internal data structures necessary for the
18: execution of the matrix function solver.
20: Collective on MFN 22: Input Parameter:
23: . mfn - matrix function context
25: Notes:
26: This function need not be called explicitly in most cases, since MFNSolve()
27: calls it. It can be useful when one wants to measure the set-up time
28: separately from the solve time.
30: Level: developer
32: .seealso: MFNCreate(), MFNSolve(), MFNDestroy()
33: @*/
34: PetscErrorCode MFNSetUp(MFN mfn) 35: {
37: PetscInt N;
42: /* reset the convergence flag from the previous solves */
43: mfn->reason = MFN_CONVERGED_ITERATING;
45: if (mfn->setupcalled) return(0);
46: PetscLogEventBegin(MFN_SetUp,mfn,0,0,0);
48: /* Set default solver type (MFNSetFromOptions was not called) */
49: if (!((PetscObject)mfn)->type_name) {
50: MFNSetType(mfn,MFNKRYLOV);
51: }
52: if (!mfn->fn) { MFNGetFN(mfn,&mfn->fn); }
53: if (!((PetscObject)mfn->fn)->type_name) {
54: FNSetFromOptions(mfn->fn);
55: }
57: /* Check problem dimensions */
58: if (!mfn->A) SETERRQ(PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_WRONGSTATE,"MFNSetOperator must be called first");
59: MatGetSize(mfn->A,&N,NULL);
60: if (mfn->ncv > N) mfn->ncv = N;
62: /* call specific solver setup */
63: (*mfn->ops->setup)(mfn);
65: /* set tolerance if not yet set */
66: if (mfn->tol==PETSC_DEFAULT) mfn->tol = SLEPC_DEFAULT_TOL;
68: PetscLogEventEnd(MFN_SetUp,mfn,0,0,0);
69: mfn->setupcalled = 1;
70: return(0);
71: }
73: /*@
74: MFNSetOperator - Sets the matrix for which the matrix function is to be computed.
76: Collective on MFN and Mat
78: Input Parameters:
79: + mfn - the matrix function context
80: - A - the problem matrix
82: Notes:
83: It must be called before MFNSetUp(). If it is called again after MFNSetUp() then
84: the MFN object is reset.
86: Level: beginner
88: .seealso: MFNSolve(), MFNSetUp(), MFNReset()
89: @*/
90: PetscErrorCode MFNSetOperator(MFN mfn,Mat A) 91: {
93: PetscInt m,n;
100: MatGetSize(A,&m,&n);
101: if (m!=n) SETERRQ(PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
102: PetscObjectReference((PetscObject)A);
103: if (mfn->setupcalled) { MFNReset(mfn); }
104: else { MatDestroy(&mfn->A); }
105: mfn->A = A;
106: mfn->setupcalled = 0;
107: return(0);
108: }
110: /*@
111: MFNGetOperator - Gets the matrix associated with the MFN object.
113: Collective on MFN and Mat
115: Input Parameter:
116: . mfn - the MFN context
118: Output Parameters:
119: . A - the matrix for which the matrix function is to be computed
121: Level: intermediate
123: .seealso: MFNSolve(), MFNSetOperator()
124: @*/
125: PetscErrorCode MFNGetOperator(MFN mfn,Mat *A)126: {
130: *A = mfn->A;
131: return(0);
132: }
134: /*@
135: MFNAllocateSolution - Allocate memory storage for common variables such
136: as the basis vectors.
138: Collective on MFN140: Input Parameters:
141: + mfn - matrix function context
142: - extra - number of additional positions, used for methods that require a
143: working basis slightly larger than ncv
145: Developers Note:
146: This is SLEPC_EXTERN because it may be required by user plugin MFN147: implementations.
149: Level: developer
150: @*/
151: PetscErrorCode MFNAllocateSolution(MFN mfn,PetscInt extra)152: {
154: PetscInt oldsize,requested;
155: Vec t;
158: requested = mfn->ncv + extra;
160: /* oldsize is zero if this is the first time setup is called */
161: BVGetSizes(mfn->V,NULL,NULL,&oldsize);
163: /* allocate basis vectors */
164: if (!mfn->V) { MFNGetBV(mfn,&mfn->V); }
165: if (!oldsize) {
166: if (!((PetscObject)(mfn->V))->type_name) {
167: BVSetType(mfn->V,BVSVEC);
168: }
169: MatCreateVecsEmpty(mfn->A,&t,NULL);
170: BVSetSizesFromVec(mfn->V,t,requested);
171: VecDestroy(&t);
172: } else {
173: BVResize(mfn->V,requested,PETSC_FALSE);
174: }
175: return(0);
176: }