← Back to branch summary

~ubuntu-branches/ubuntu/feisty/petsc/feisty

« back to all changes in this revision

Viewing changes to src/sles/ksp/impls/minres/minres.c

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2005-03-24 09:46:23 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050324094623-dfcxn8bltjms2cqq
Tags: 2.2.0-4
http://bugs.debian.org/284023
http://bugs.debian.org/293011
http://bugs.debian.org/274009
http://bugs.debian.org/267796
* Update for new mpich >> 1.2.5.3-2.
* Fixed src/inline/axpy.h for complex and UNROLL (closes: #284023, #293011).
* Added -fno-strict-aliasing to compile flags (closes: #274009).
* Switched SLES stuff to KSP in petsc.m4 (closes: 267796).

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/sles/ksp/impls/minres/minres.c
1
 
/*$Id: minres.c,v 1.18 2001/08/10 18:14:38 bsmith Exp $*/
2
 
/*                       
3
 
    This code implements the MINRES (Minimum Residual) method. 
4
 
    Reference: Paige & Saunders, 1975.
5
 
 
6
 
    Contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
7
 
 
8
 
*/
9
 
#include "src/sles/ksp/kspimpl.h"
10
 
 
11
 
typedef struct {
12
 
  PetscReal haptol;
13
 
} KSP_MINRES;
14
 
 
15
 
#undef __FUNCT__  
16
 
#define __FUNCT__ "KSPSetUp_MINRES"
17
 
int KSPSetUp_MINRES(KSP ksp)
18
 
{
19
 
  int ierr;
20
 
 
21
 
  PetscFunctionBegin;
22
 
 
23
 
  if (ksp->pc_side == PC_RIGHT) {
24
 
    SETERRQ(2,"No right preconditioning for KSPMINRES");
25
 
  } else if (ksp->pc_side == PC_SYMMETRIC) {
26
 
    SETERRQ(2,"No symmetric preconditioning for KSPMINRES");
27
 
  }
28
 
 
29
 
  ierr = KSPDefaultGetWork(ksp,9);CHKERRQ(ierr);
30
 
 
31
 
  PetscFunctionReturn(0);
32
 
}
33
 
 
34
 
 
35
 
#undef __FUNCT__  
36
 
#define __FUNCT__ "KSPSolve_MINRES"
37
 
int  KSPSolve_MINRES(KSP ksp,int *its)
38
 
{
39
 
  int          ierr,i,maxit;
40
 
  PetscScalar  alpha,malpha,beta,mbeta,ibeta,betaold,eta,c=1.0,ceta,cold=1.0,coold,s=0.0,sold=0.0,soold;
41
 
  PetscScalar  rho0,rho1,irho1,rho2,mrho2,rho3,mrho3,mone = -1.0,zero = 0.0,dp = 0.0;
42
 
  PetscReal    np;
43
 
  Vec          X,B,R,Z,U,V,W,UOLD,VOLD,WOLD,WOOLD;
44
 
  Mat          Amat,Pmat;
45
 
  MatStructure pflag;
46
 
  KSP_MINRES   *minres = (KSP_MINRES*)ksp->data;
47
 
  PetscTruth   diagonalscale;
48
 
 
49
 
  PetscFunctionBegin;
50
 
  ierr    = PCDiagonalScale(ksp->B,&diagonalscale);CHKERRQ(ierr);
51
 
  if (diagonalscale) SETERRQ1(1,"Krylov method %s does not support diagonal scaling",ksp->type_name);
52
 
 
53
 
  maxit   = ksp->max_it;
54
 
  X       = ksp->vec_sol;
55
 
  B       = ksp->vec_rhs;
56
 
  R       = ksp->work[0];
57
 
  Z       = ksp->work[1];
58
 
  U       = ksp->work[2];
59
 
  V       = ksp->work[3];
60
 
  W       = ksp->work[4];
61
 
  UOLD    = ksp->work[5];
62
 
  VOLD    = ksp->work[6];
63
 
  WOLD    = ksp->work[7];
64
 
  WOOLD   = ksp->work[8];
65
 
 
66
 
  ierr = PCGetOperators(ksp->B,&Amat,&Pmat,&pflag);CHKERRQ(ierr);
67
 
 
68
 
  ksp->its = 0;
69
 
 
70
 
  ierr = VecSet(&zero,UOLD);CHKERRQ(ierr);         /*     u_old  <-   0   */
71
 
  ierr = VecCopy(UOLD,VOLD);CHKERRQ(ierr);         /*     v_old  <-   0   */
72
 
  ierr = VecCopy(UOLD,W);CHKERRQ(ierr);            /*     w      <-   0   */
73
 
  ierr = VecCopy(UOLD,WOLD);CHKERRQ(ierr);         /*     w_old  <-   0   */
74
 
 
75
 
  if (!ksp->guess_zero) {
76
 
    ierr = KSP_MatMult(ksp,Amat,X,R);CHKERRQ(ierr); /*     r <- b - A*x    */
77
 
    ierr = VecAYPX(&mone,B,R);CHKERRQ(ierr);
78
 
  } else { 
79
 
    ierr = VecCopy(B,R);CHKERRQ(ierr);              /*     r <- b (x is 0) */
80
 
  }
81
 
 
82
 
  ierr = KSP_PCApply(ksp,ksp->B,R,Z);CHKERRQ(ierr); /*     z  <- B*r       */
83
 
 
84
 
  ierr = VecDot(R,Z,&dp);CHKERRQ(ierr);
85
 
  if (PetscAbsScalar(dp) < minres->haptol) {
86
 
    PetscLogInfo(ksp,"KSPSolve_MINRES:Detected happy breakdown %g tolerance %g\n",PetscAbsScalar(dp),minres->haptol);
87
 
    dp = PetscAbsScalar(dp); /* tiny number, can't use 0.0, cause divided by below */
88
 
    if (dp == 0.0) {
89
 
      ksp->reason = KSP_CONVERGED_ATOL;
90
 
      *its        = 0;
91
 
      PetscFunctionReturn(0);
92
 
    }
93
 
  }
94
 
 
95
 
#if !defined(PETSC_USE_COMPLEX)
96
 
  if (dp < 0.0) SETERRQ(PETSC_ERR_KSP_BRKDWN,"Indefinite preconditioner");
97
 
#endif
98
 
  dp   = PetscSqrtScalar(dp); 
99
 
  beta = dp;                                        /*  beta <- sqrt(r'*z  */
100
 
  eta  = beta;
101
 
 
102
 
  ierr = VecCopy(R,V);CHKERRQ(ierr);
103
 
  ierr = VecCopy(Z,U);CHKERRQ(ierr);
104
 
  ibeta = 1.0 / beta;
105
 
  ierr = VecScale(&ibeta,V);CHKERRQ(ierr);         /*    v <- r / beta     */
106
 
  ierr = VecScale(&ibeta,U);CHKERRQ(ierr);         /*    u <- z / beta     */
107
 
 
108
 
  ierr = VecNorm(Z,NORM_2,&np);CHKERRQ(ierr);      /*   np <- ||z||        */
109
 
 
110
 
  ierr = (*ksp->converged)(ksp,0,np,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);  /* test for convergence */
111
 
  if (ksp->reason) {*its =  0; PetscFunctionReturn(0);}
112
 
  KSPLogResidualHistory(ksp,np);
113
 
  KSPMonitor(ksp,0,np);            /* call any registered monitor routines */
114
 
  ksp->rnorm = np;  
115
 
 
116
 
  for (i=0; i<maxit; i++) {
117
 
     ksp->its = i+1;
118
 
 
119
 
/*   Lanczos  */
120
 
 
121
 
     ierr = KSP_MatMult(ksp,Amat,U,R);CHKERRQ(ierr);   /*      r <- A*u   */
122
 
     ierr = VecDot(U,R,&alpha);CHKERRQ(ierr);          /*  alpha <- r'*u  */
123
 
     ierr = KSP_PCApply(ksp,ksp->B,R,Z);CHKERRQ(ierr); /*      z <- B*r   */
124
 
 
125
 
     malpha = - alpha;
126
 
     ierr = VecAXPY(&malpha,V,R);CHKERRQ(ierr);     /*  r <- r - alpha v     */
127
 
     ierr = VecAXPY(&malpha,U,Z);CHKERRQ(ierr);     /*  z <- z - alpha u     */
128
 
     mbeta = - beta;
129
 
     ierr = VecAXPY(&mbeta,VOLD,R);CHKERRQ(ierr);   /*  r <- r - beta v_old  */
130
 
     ierr = VecAXPY(&mbeta,UOLD,Z);CHKERRQ(ierr);   /*  z <- z - beta u_old  */
131
 
 
132
 
     betaold = beta;
133
 
 
134
 
     ierr = VecDot(R,Z,&dp);CHKERRQ(ierr); 
135
 
     if (PetscAbsScalar(dp) < minres->haptol) {
136
 
       PetscLogInfo(ksp,"KSPSolve_MINRES:Detected happy breakdown %g tolerance %g\n",PetscAbsScalar(dp),minres->haptol);
137
 
       dp = PetscAbsScalar(dp); /* tiny number, can we use 0.0? */
138
 
     }
139
 
 
140
 
#if !defined(PETSC_USE_COMPLEX)
141
 
     if (dp < 0.0) SETERRQ1(PETSC_ERR_KSP_BRKDWN,"Indefinite preconditioner R'Z = %g",dp);
142
 
#endif
143
 
     beta = PetscSqrtScalar(dp);                               /*  beta <- sqrt(r'*z)   */
144
 
 
145
 
/*    QR factorisation    */
146
 
 
147
 
     coold = cold; cold = c; soold = sold; sold = s;
148
 
 
149
 
     rho0 = cold * alpha - coold * sold * betaold;
150
 
     rho1 = PetscSqrtScalar(rho0*rho0 + beta*beta);
151
 
     rho2 = sold * alpha + coold * cold * betaold;
152
 
     rho3 = soold * betaold;
153
 
 
154
 
/*     Givens rotation    */
155
 
 
156
 
     c = rho0 / rho1;
157
 
     s = beta / rho1;
158
 
 
159
 
/*    Update    */
160
 
 
161
 
     ierr = VecCopy(WOLD,WOOLD);CHKERRQ(ierr);     /*  w_oold <- w_old      */
162
 
     ierr = VecCopy(W,WOLD);CHKERRQ(ierr);         /*  w_old  <- w          */
163
 
     
164
 
     ierr = VecCopy(U,W);CHKERRQ(ierr);            /*  w      <- u          */
165
 
     mrho2 = - rho2;
166
 
     ierr = VecAXPY(&mrho2,WOLD,W);CHKERRQ(ierr);  /*  w <- w - rho2 w_old  */
167
 
     mrho3 = - rho3;
168
 
     ierr = VecAXPY(&mrho3,WOOLD,W);CHKERRQ(ierr); /*  w <- w - rho3 w_oold */
169
 
     irho1 = 1.0 / rho1;
170
 
     ierr = VecScale(&irho1,W);CHKERRQ(ierr);      /*  w <- w / rho1        */
171
 
 
172
 
     ceta = c * eta;
173
 
     ierr = VecAXPY(&ceta,W,X);CHKERRQ(ierr);      /*  x <- x + c eta w     */ 
174
 
     eta = - s * eta;
175
 
 
176
 
     ierr = VecCopy(V,VOLD);CHKERRQ(ierr);
177
 
     ierr = VecCopy(U,UOLD);CHKERRQ(ierr);
178
 
     ierr = VecCopy(R,V);CHKERRQ(ierr);
179
 
     ierr = VecCopy(Z,U);CHKERRQ(ierr);
180
 
     ibeta = 1.0 / beta;
181
 
     ierr = VecScale(&ibeta,V);CHKERRQ(ierr);      /*  v <- r / beta       */
182
 
     ierr = VecScale(&ibeta,U);CHKERRQ(ierr);      /*  u <- z / beta       */
183
 
     
184
 
     np = ksp->rnorm * PetscAbsScalar(s);
185
 
 
186
 
     ksp->rnorm = np;
187
 
     KSPLogResidualHistory(ksp,np);
188
 
     KSPMonitor(ksp,i+1,np);
189
 
     ierr = (*ksp->converged)(ksp,i+1,np,&ksp->reason,ksp->cnvP);CHKERRQ(ierr); /* test for convergence */
190
 
     if (ksp->reason) break;
191
 
  }
192
 
  if (i == maxit) {
193
 
    ksp->reason = KSP_DIVERGED_ITS;
194
 
  }
195
 
  *its = ksp->its;
196
 
  PetscFunctionReturn(0);
197
 
}
198
 
 
199
 
EXTERN_C_BEGIN
200
 
#undef __FUNCT__  
201
 
#define __FUNCT__ "KSPCreate_MINRES"
202
 
int KSPCreate_MINRES(KSP ksp)
203
 
{
204
 
  KSP_MINRES *minres;
205
 
  int ierr;
206
 
 
207
 
  PetscFunctionBegin;
208
 
 
209
 
  ksp->pc_side   = PC_LEFT;
210
 
  ierr           = PetscNew(KSP_MINRES,&minres);CHKERRQ(ierr);
211
 
  minres->haptol = 1.e-18;
212
 
  ksp->data      = (void*)minres;
213
 
 
214
 
  /*
215
 
       Sets the functions that are associated with this data structure 
216
 
       (in C++ this is the same as defining virtual functions)
217
 
  */
218
 
  ksp->ops->setup                = KSPSetUp_MINRES;
219
 
  ksp->ops->solve                = KSPSolve_MINRES;
220
 
  ksp->ops->destroy              = KSPDefaultDestroy;
221
 
  ksp->ops->setfromoptions       = 0;
222
 
  ksp->ops->buildsolution        = KSPDefaultBuildSolution;
223
 
  ksp->ops->buildresidual        = KSPDefaultBuildResidual;
224
 
 
225
 
  PetscFunctionReturn(0);
226
 
}
227
 
EXTERN_C_END
228
 
 
229
 
 
230
 
 
231
 
 
232