2: #include <petsc/private/kspimpl.h>
4: /* ---------------------------------------Method 1------------------------------------------------------------*/
5: typedef struct {
6: PetscInt method; /* 1 or 2 */
7: PetscInt curl; /* Current number of basis vectors */
8: PetscInt maxl; /* Maximum number of basis vectors */
9: PetscInt refcnt;
10: PetscBool monitor;
11: Mat mat;
12: KSP ksp;
13: PetscScalar *alpha; /* */
14: Vec *xtilde; /* Saved x vectors */
15: Vec *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: PetscNew(&itg);
30: PetscMalloc1(maxl,&itg->alpha);
31: PetscLogObjectMemory((PetscObject)ksp,sizeof(KSPFischerGuess_Method1) + maxl*sizeof(PetscScalar));
32: KSPCreateVecs(ksp,maxl,&itg->xtilde,0,NULL);
33: PetscLogObjectParents(ksp,maxl,itg->xtilde);
34: KSPCreateVecs(ksp,maxl,&itg->btilde,0,NULL);
35: PetscLogObjectParents(ksp,maxl,itg->btilde);
36: VecDuplicate(itg->xtilde[0],&itg->guess);
37: PetscLogObjectParent((PetscObject)ksp,(PetscObject)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 ",(double)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\n");
124: }
125: }
126: return(0);
127: }
129: /* ---------------------------------------Method 2------------------------------------------------------------*/
130: typedef struct {
131: PetscInt method; /* 1 or 2 */
132: PetscInt curl; /* Current number of basis vectors */
133: PetscInt maxl; /* Maximum number of basis vectors */
134: PetscInt 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: PetscNew(&itg);
153: PetscMalloc1(maxl,&itg->alpha);
154: PetscLogObjectMemory((PetscObject)ksp,sizeof(KSPFischerGuess_Method2) + maxl*sizeof(PetscScalar));
155: KSPCreateVecs(ksp,maxl,&itg->xtilde,0,NULL);
156: PetscLogObjectParents(ksp,maxl,itg->xtilde);
157: VecDuplicate(itg->xtilde[0],&itg->Ax);
158: PetscLogObjectParent((PetscObject)ksp,(PetscObject)itg->Ax);
159: VecDuplicate(itg->xtilde[0],&itg->guess);
160: PetscLogObjectParent((PetscObject)ksp,(PetscObject)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 ",(double)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\n");
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: . 1. - 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: {
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(PetscObjectComm((PetscObject)ksp),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;
294: KSPGetOperators(ksp,&(*itg)->mat,NULL);
295: return(0);
296: }
300: PetscErrorCode KSPFischerGuessSetFromOptions(KSPFischerGuess ITG)301: {
305: PetscOptionsGetBool(((PetscObject)ITG->ksp)->options,((PetscObject)ITG->ksp)->prefix,"-ksp_fischer_guess_monitor",&ITG->monitor,NULL);
306: return(0);
307: }
311: PetscErrorCode KSPFischerGuessDestroy(KSPFischerGuess *ITG)312: {
316: if (!*ITG) return(0);
317: if (--(*ITG)->refcnt) {*ITG = 0; return(0);}
319: if ((*ITG)->method == 1) {
320: KSPFischerGuessDestroy_Method1((KSPFischerGuess_Method1*)*ITG);
321: } else if ((*ITG)->method == 2) {
322: KSPFischerGuessDestroy_Method2((KSPFischerGuess_Method2*)*ITG);
323: } else SETERRQ(((PetscObject)(*ITG)->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
324: *ITG = NULL;
325: return(0);
326: }
330: PetscErrorCode KSPFischerGuessUpdate(KSPFischerGuess itg,Vec x)331: {
335: if (itg->method == 1) {
336: KSPFischerGuessUpdate_Method1((KSPFischerGuess_Method1*)itg,x);
337: } else if (itg->method == 2) {
338: KSPFischerGuessUpdate_Method2((KSPFischerGuess_Method2*)itg,x);
339: } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
340: return(0);
341: }
345: PetscErrorCode KSPFischerGuessFormGuess(KSPFischerGuess itg,Vec b,Vec x)346: {
350: if (itg->method == 1) {
351: KSPFischerGuessFormGuess_Method1((KSPFischerGuess_Method1*)itg,b,x);
352: } else if (itg->method == 2) {
353: KSPFischerGuessFormGuess_Method2((KSPFischerGuess_Method2*)itg,b,x);
354: } else SETERRQ(((PetscObject)itg->ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Method can only be 1 or 2");
355: return(0);
356: }
360: /*
361: KSPFischerGuessReset - This is called whenever KSPSetOperators() is called to tell the
362: initial guess object that the matrix is changed and so the initial guess object
363: must restart from scratch building the subspace where the guess is computed from.
364: */
365: PetscErrorCode KSPFischerGuessReset(KSPFischerGuess itg)366: {
370: itg->curl = 0;
371: KSPGetOperators(itg->ksp,&itg->mat,NULL);
372: return(0);
373: }