2: #include <petsc-private/kspimpl.h>
4: /* ---------------------------------------Method 1------------------------------------------------------------*/
5: typedef struct {
6: PetscInt method, /* 1 or 2 */
7: curl, /* Current number of basis vectors */
8: maxl, /* Maximum number of basis vectors */
9: refcnt;
10: PetscBool monitor;
11: Mat mat;
12: KSP ksp;
13: PetscScalar *alpha; /* */
14: Vec *xtilde, /* Saved x vectors */
15: *btilde; /* Saved b vectors */
16: Vec guess;
17: } KSPFischerGuess_Method1;
22: PetscErrorCode KSPFischerGuessCreate_Method1(KSP ksp,int maxl,KSPFischerGuess_Method1 **ITG) 23: {
24: KSPFischerGuess_Method1 *itg;
25: PetscErrorCode ierr;
29: PetscMalloc(sizeof(KSPFischerGuess_Method1),&itg);
30: PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);
31: PetscLogObjectMemory(ksp,sizeof(KSPFischerGuess_Method1) + maxl*sizeof(PetscScalar));
32: KSPGetVecs(ksp,maxl,&itg->xtilde,0,PETSC_NULL);
33: PetscLogObjectParents(ksp,maxl,itg->xtilde);
34: KSPGetVecs(ksp,maxl,&itg->btilde,0,PETSC_NULL);
35: PetscLogObjectParents(ksp,maxl,itg->btilde);
36: VecDuplicate(itg->xtilde[0],&itg->guess);
37: PetscLogObjectParent(ksp,itg->guess);
38: *ITG = itg;
39: return(0);
40: }
45: PetscErrorCode KSPFischerGuessDestroy_Method1(KSPFischerGuess_Method1 *itg) 46: {
50: PetscFree(itg->alpha);
51: VecDestroyVecs(itg->maxl,&itg->btilde);
52: VecDestroyVecs(itg->maxl,&itg->xtilde);
53: VecDestroy(&itg->guess);
54: PetscFree(itg);
55: return(0);
56: }
59: /*
60: Given a basis generated already this computes a new guess x from the new right hand side b
61: Figures out the components of b in each btilde direction and adds them to x
62: */
65: PetscErrorCode KSPFischerGuessFormGuess_Method1(KSPFischerGuess_Method1 *itg,Vec b,Vec x) 66: {
68: PetscInt i;
73: VecSet(x,0.0);
74: VecMDot(b,itg->curl,itg->btilde,itg->alpha);
75: if (itg->monitor) {
76: PetscPrintf(((PetscObject)itg->ksp)->comm,"KSPFischerGuess alphas = ");
77: for (i=0; i<itg->curl; i++ ){
78: PetscPrintf(((PetscObject)itg->ksp)->comm,"%G ",PetscAbsScalar(itg->alpha[i]));
79: }
80: PetscPrintf(((PetscObject)itg->ksp)->comm,"\n");
81: }
82: VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
83: VecCopy(x,itg->guess);
84: /* Note: do not change the b right hand side as is done in the publication */
85: return(0);
86: }
90: PetscErrorCode KSPFischerGuessUpdate_Method1(KSPFischerGuess_Method1 *itg,Vec x) 91: {
92: PetscReal norm;
94: int curl = itg->curl,i;
99: if (curl == itg->maxl) {
100: KSP_MatMult(itg->ksp,itg->mat,x,itg->btilde[0]);
101: VecNormalize(itg->btilde[0],&norm);
102: VecCopy(x,itg->xtilde[0]);
103: VecScale(itg->xtilde[0],1.0/norm);
104: itg->curl = 1;
105: } else {
106: if (!curl) {
107: VecCopy(x,itg->xtilde[curl]);
108: } else {
109: VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
110: }
112: KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->btilde[curl]);
113: VecMDot(itg->btilde[curl],curl,itg->btilde,itg->alpha);
114: for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
115: VecMAXPY(itg->btilde[curl],curl,itg->alpha,itg->btilde);
116: VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
118: VecNormalize(itg->btilde[curl],&norm);
119: if (norm) {
120: VecScale(itg->xtilde[curl],1.0/norm);
121: itg->curl++;
122: } else {
123: PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous");
124: }
125: }
126: return(0);
127: }
129: /* ---------------------------------------Method 2------------------------------------------------------------*/
130: typedef struct {
131: PetscInt method, /* 1 or 2 */
132: curl, /* Current number of basis vectors */
133: maxl, /* Maximum number of basis vectors */
134: refcnt;
135: PetscBool monitor;
136: Mat mat;
137: KSP ksp;
138: PetscScalar *alpha; /* */
139: Vec *xtilde; /* Saved x vectors */
140: Vec Ax,guess;
141: } KSPFischerGuess_Method2;
145: PetscErrorCode KSPFischerGuessCreate_Method2(KSP ksp,int maxl,KSPFischerGuess_Method2 **ITG)146: {
147: KSPFischerGuess_Method2 *itg;
148: PetscErrorCode ierr;
152: PetscMalloc(sizeof(KSPFischerGuess_Method2),&itg);
153: PetscMalloc(maxl * sizeof(PetscScalar),&itg->alpha);
154: PetscLogObjectMemory(ksp,sizeof(KSPFischerGuess_Method2) + maxl*sizeof(PetscScalar));
155: KSPGetVecs(ksp,maxl,&itg->xtilde,0,PETSC_NULL);
156: PetscLogObjectParents(ksp,maxl,itg->xtilde);
157: VecDuplicate(itg->xtilde[0],&itg->Ax);
158: PetscLogObjectParent(ksp,itg->Ax);
159: VecDuplicate(itg->xtilde[0],&itg->guess);
160: PetscLogObjectParent(ksp,itg->guess);
161: *ITG = itg;
162: return(0);
163: }
168: PetscErrorCode KSPFischerGuessDestroy_Method2(KSPFischerGuess_Method2 *itg)169: {
173: PetscFree(itg->alpha);
174: VecDestroyVecs(itg->maxl,&itg->xtilde);
175: VecDestroy(&itg->Ax);
176: VecDestroy(&itg->guess);
177: PetscFree(itg);
178: return(0);
179: }
182: /*
183: Given a basis generated already this computes a new guess x from the new right hand side b
184: Figures out the components of b in each btilde direction and adds them to x
185: */
188: PetscErrorCode KSPFischerGuessFormGuess_Method2(KSPFischerGuess_Method2 *itg,Vec b,Vec x)189: {
191: PetscInt i;
196: VecSet(x,0.0);
197: VecMDot(b,itg->curl,itg->xtilde,itg->alpha);
198: if (itg->monitor) {
199: PetscPrintf(((PetscObject)itg->ksp)->comm,"KSPFischerGuess alphas = ");
200: for (i=0; i<itg->curl; i++ ){
201: PetscPrintf(((PetscObject)itg->ksp)->comm,"%G ",PetscAbsScalar(itg->alpha[i]));
202: }
203: PetscPrintf(((PetscObject)itg->ksp)->comm,"\n");
204: }
205: VecMAXPY(x,itg->curl,itg->alpha,itg->xtilde);
206: VecCopy(x,itg->guess);
207: /* Note: do not change the b right hand side as is done in the publication */
208: return(0);
209: }
213: PetscErrorCode KSPFischerGuessUpdate_Method2(KSPFischerGuess_Method2 *itg,Vec x)214: {
215: PetscScalar norm;
217: int curl = itg->curl,i;
222: if (curl == itg->maxl) {
223: KSP_MatMult(itg->ksp,itg->mat,x,itg->Ax); /* norm = sqrt(x'Ax) */
224: VecDot(x,itg->Ax,&norm);
225: VecCopy(x,itg->xtilde[0]);
226: VecScale(itg->xtilde[0],1.0/PetscSqrtScalar(norm));
227: itg->curl = 1;
228: } else {
229: if (!curl) {
230: VecCopy(x,itg->xtilde[curl]);
231: } else {
232: VecWAXPY(itg->xtilde[curl],-1.0,itg->guess,x);
233: }
234: KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax);
235: VecMDot(itg->Ax,curl,itg->xtilde,itg->alpha);
236: for (i=0; i<curl; i++) itg->alpha[i] = -itg->alpha[i];
237: VecMAXPY(itg->xtilde[curl],curl,itg->alpha,itg->xtilde);
239: KSP_MatMult(itg->ksp,itg->mat,itg->xtilde[curl],itg->Ax); /* norm = sqrt(xtilde[curl]'Axtilde[curl]) */
240: VecDot(itg->xtilde[curl],itg->Ax,&norm);
241: if (PetscAbsScalar(norm) != 0.0) {
242: VecScale(itg->xtilde[curl],1.0/PetscSqrtScalar(norm));
243: itg->curl++;
244: } else {
245: PetscInfo(itg->ksp,"Not increasing dimension of Fischer space because new direction is identical to previous");
246: }
247: }
248: return(0);
249: }
251: /* ---------------------------------------------------------------------------------------------------------*/
254: /*@C
255: KSPFischerGuessCreate - Implements Paul Fischer's initial guess algorithm Method 1 and 2 for situations where
256: a linear system is solved repeatedly
258: References:
259: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940020363_1994020363.pdf
261: Notes: the algorithm is different from the paper because we do not CHANGE the right hand side of the new
262: problem and solve the problem with an initial guess of zero, rather we solve the original new problem
263: with a nonzero initial guess (this is done so that the linear solver convergence tests are based on
264: the original RHS.) But we use the xtilde = x - xguess as the new direction so that it is not
265: mostly orthogonal to the previous solutions.
267: These are not intended to be used directly, they are called by KSP automatically when the
268: KSP option KSPSetFischerGuess(KSP,PetscInt,PetscInt) or -ksp_guess_fischer <int,int>
270: Method 2 is only for positive definite matrices, since it uses the A norm.
272: This is not currently programmed as a PETSc class because there are only two methods; if more methods
273: are introduced it should be changed. For example the Knoll guess should be included
275: Level: advanced
277: @*/
278: PetscErrorCodeKSPFischerGuessCreate(KSP ksp,PetscInt method,PetscInt maxl,KSPFischerGuess *itg)279: {
280: PetscErrorCode ierr;
283: if (method == 1) {
284: KSPFischerGuessCreate_Method1(ksp,maxl,(KSPFischerGuess_Method1 **)itg);
285: } else if (method == 2) {
286: KSPFischerGuessCreate_Method2(ksp,maxl,(KSPFischerGuess_Method2 **)itg);
287: } else SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
288: (*itg)->method = method;
289: (*itg)->curl = 0;
290: (*itg)->maxl = maxl;
291: (*itg)->ksp = ksp;
292: (*itg)->refcnt = 1;
293: KSPGetOperators(ksp,&(*itg)->mat,PETSC_NULL,PETSC_NULL);
294: return(0);
295: }
299: PetscErrorCode KSPFischerGuessSetFromOptions(KSPFischerGuess ITG)300: {
301: PetscErrorCode ierr;
304: PetscOptionsGetBool(((PetscObject)ITG->ksp)->prefix,"-ksp_fischer_guess_monitor",&ITG->monitor,PETSC_NULL);
305: return(0);
306: }
310: PetscErrorCode KSPFischerGuessDestroy(KSPFischerGuess *ITG)311: {
312: PetscErrorCode ierr;
315: if (!*ITG) return(0);
316: if (--(*ITG)->refcnt) {*ITG = 0; return(0);}
318: if ((*ITG)->method == 1) {
319: KSPFischerGuessDestroy_Method1((KSPFischerGuess_Method1 *)*ITG);
320: } else if ((*ITG)->method == 2) {
321: KSPFischerGuessDestroy_Method2((KSPFischerGuess_Method2 *)*ITG);
322: } else SETERRQ(((PetscObject)(*ITG)->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
323: *ITG = PETSC_NULL;
324: return(0);
325: }
329: PetscErrorCode KSPFischerGuessUpdate(KSPFischerGuess itg,Vec x)330: {
331: PetscErrorCode ierr;
334: if (itg->method == 1) {
335: KSPFischerGuessUpdate_Method1((KSPFischerGuess_Method1 *)itg,x);
336: } else if (itg->method == 2) {
337: KSPFischerGuessUpdate_Method2((KSPFischerGuess_Method2 *)itg,x);
338: } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
339: return(0);
340: }
344: PetscErrorCode KSPFischerGuessFormGuess(KSPFischerGuess itg,Vec b,Vec x)345: {
346: PetscErrorCode ierr;
349: if (itg->method == 1) {
350: KSPFischerGuessFormGuess_Method1((KSPFischerGuess_Method1 *)itg,b,x);
351: } else if (itg->method == 2) {
352: KSPFischerGuessFormGuess_Method2((KSPFischerGuess_Method2 *)itg,b,x);
353: } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
354: return(0);
355: }
359: /*
360: KSPFischerGuessReset - This is called whenever KSPSetOperators() is called to tell the
361: initial guess object that the matrix is changed and so the initial guess object
362: must restart from scratch building the subspace where the guess is computed from.
363: */
364: PetscErrorCode KSPFischerGuessReset(KSPFischerGuess itg)365: {
369: itg->curl = 0;
370: KSPGetOperators(itg->ksp,&itg->mat,PETSC_NULL,PETSC_NULL);
371: return(0);
372: }