← Back to branch summary

~ubuntu-branches/ubuntu/warty/petsc/warty

« back to all changes in this revision

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

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2004-06-07 13:41:43 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20040607134143-92p586zrauvie0le
Tags: 2.2.0-2
http://bugs.debian.org/249617
* Upstream patch level 2.
* New PETSC_BOPT_EXTRA option for different BOPT and lib names, with _c++
  symlinks only for plain and single (closes: #249617).
* New DEBIAN_DIST=contrib option to link with hypre, parmetis (closes:
  #249619).
* Combined petsc-c and petsc-fortran substvars into petsc-compilers.
* Extra quote in -dev prerm eliminates "too many arguments" problem.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/sles/ksp/impls/tfqmr/tfqmr.c
1
 
/*$Id: tfqmr.c,v 1.61 2001/08/07 03:03:54 balay Exp $*/
2
 
 
3
 
/*                       
4
 
    This code implements the TFQMR (Transpose-free variant of Quasi-Minimal
5
 
    Residual) method.  Reference: Freund, 1993
6
 
 
7
 
    Note that for the complex numbers version, the VecDot() arguments
8
 
    within the code MUST remain in the order given for correct computation
9
 
    of inner products.
10
 
*/
11
 
 
12
 
#include "src/sles/ksp/kspimpl.h"
13
 
 
14
 
#undef __FUNCT__  
15
 
#define __FUNCT__ "KSPSetUp_TFQMR"
16
 
static int KSPSetUp_TFQMR(KSP ksp)
17
 
{
18
 
  int ierr;
19
 
 
20
 
  PetscFunctionBegin;
21
 
  if (ksp->pc_side == PC_SYMMETRIC){
22
 
    SETERRQ(2,"no symmetric preconditioning for KSPTFQMR");
23
 
  }
24
 
  ierr = KSPDefaultGetWork(ksp,9);CHKERRQ(ierr);
25
 
  PetscFunctionReturn(0);
26
 
}
27
 
 
28
 
#undef __FUNCT__  
29
 
#define __FUNCT__ "KSPSolve_TFQMR"
30
 
static int  KSPSolve_TFQMR(KSP ksp,int *its)
31
 
{
32
 
  int       i,maxit,m, ierr;
33
 
  PetscScalar    rho,rhoold,a,s,b,eta,etaold,psiold,cf,tmp,one = 1.0,zero = 0.0;
34
 
  PetscReal dp,dpold,w,dpest,tau,psi,cm;
35
 
  Vec       X,B,V,P,R,RP,T,T1,Q,U,D,AUQ;
36
 
 
37
 
  PetscFunctionBegin;
38
 
  maxit    = ksp->max_it;
39
 
  X        = ksp->vec_sol;
40
 
  B        = ksp->vec_rhs;
41
 
  R        = ksp->work[0];
42
 
  RP       = ksp->work[1];
43
 
  V        = ksp->work[2];
44
 
  T        = ksp->work[3];
45
 
  Q        = ksp->work[4];
46
 
  P        = ksp->work[5];
47
 
  U        = ksp->work[6];
48
 
  D        = ksp->work[7];
49
 
  T1       = ksp->work[8];
50
 
  AUQ      = V;
51
 
 
52
 
  /* Compute initial preconditioned residual */
53
 
  ierr = KSPInitialResidual(ksp,X,V,T,R,B);CHKERRQ(ierr);
54
 
 
55
 
  /* Test for nothing to do */
56
 
  ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
57
 
  ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr);
58
 
  ksp->rnorm  = dp;
59
 
  ksp->its    = 0;
60
 
  ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr);
61
 
  ierr = (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
62
 
  if (ksp->reason) {*its = 0; PetscFunctionReturn(0);}
63
 
  KSPMonitor(ksp,0,dp);
64
 
 
65
 
  /* Make the initial Rp == R */
66
 
  ierr = VecCopy(R,RP);CHKERRQ(ierr);
67
 
 
68
 
  /* Set the initial conditions */
69
 
  etaold = 0.0;
70
 
  psiold = 0.0;
71
 
  tau    = dp;
72
 
  dpold  = dp;
73
 
 
74
 
  ierr = VecDot(R,RP,&rhoold);CHKERRQ(ierr);       /* rhoold = (r,rp)     */
75
 
  ierr = VecCopy(R,U);CHKERRQ(ierr);
76
 
  ierr = VecCopy(R,P);CHKERRQ(ierr);
77
 
  ierr = KSP_PCApplyBAorAB(ksp,ksp->B,ksp->pc_side,P,V,T);CHKERRQ(ierr);
78
 
  ierr = VecSet(&zero,D);CHKERRQ(ierr);
79
 
 
80
 
  for (i=0; i<maxit; i++) {
81
 
    ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr);
82
 
    ksp->its++;
83
 
    ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr);
84
 
    ierr = VecDot(V,RP,&s);CHKERRQ(ierr);          /* s <- (v,rp)          */
85
 
    a = rhoold / s;                                 /* a <- rho / s         */
86
 
    tmp = -a; VecWAXPY(&tmp,V,U,Q);CHKERRQ(ierr);  /* q <- u - a v         */
87
 
    ierr = VecWAXPY(&one,U,Q,T);CHKERRQ(ierr);     /* t <- u + q           */
88
 
    ierr = KSP_PCApplyBAorAB(ksp,ksp->B,ksp->pc_side,T,AUQ,T1);CHKERRQ(ierr);
89
 
    ierr = VecAXPY(&tmp,AUQ,R);CHKERRQ(ierr);      /* r <- r - a K (u + q) */
90
 
    ierr = VecNorm(R,NORM_2,&dp);CHKERRQ(ierr);
91
 
    for (m=0; m<2; m++) {
92
 
      if (!m) {
93
 
        w = sqrt(dp*dpold);
94
 
      } else {
95
 
        w = dp;
96
 
      }
97
 
      psi = w / tau;
98
 
      cm  = 1.0 / sqrt(1.0 + psi * psi);
99
 
      tau = tau * psi * cm;
100
 
      eta = cm * cm * a;
101
 
      cf  = psiold * psiold * etaold / a;
102
 
      if (!m) {
103
 
        ierr = VecAYPX(&cf,U,D);CHKERRQ(ierr);
104
 
      } else {
105
 
        ierr = VecAYPX(&cf,Q,D);CHKERRQ(ierr);
106
 
      }
107
 
      ierr = VecAXPY(&eta,D,X);CHKERRQ(ierr);
108
 
 
109
 
      dpest = sqrt(m + 1.0) * tau;
110
 
      ierr = PetscObjectTakeAccess(ksp);CHKERRQ(ierr);
111
 
      ksp->rnorm                                    = dpest;
112
 
      ierr = PetscObjectGrantAccess(ksp);CHKERRQ(ierr);
113
 
      KSPLogResidualHistory(ksp,dpest);
114
 
      KSPMonitor(ksp,i+1,dpest);
115
 
      ierr = (*ksp->converged)(ksp,i+1,dpest,&ksp->reason,ksp->cnvP);CHKERRQ(ierr);
116
 
      if (ksp->reason) break;
117
 
 
118
 
      etaold = eta;
119
 
      psiold = psi;
120
 
    }
121
 
    if (ksp->reason) break;
122
 
 
123
 
    ierr = VecDot(R,RP,&rho);CHKERRQ(ierr);        /* rho <- (r,rp)       */
124
 
    b = rho / rhoold;                               /* b <- rho / rhoold   */
125
 
    ierr = VecWAXPY(&b,Q,R,U);CHKERRQ(ierr);       /* u <- r + b q        */
126
 
    ierr = VecAXPY(&b,P,Q);CHKERRQ(ierr);
127
 
    ierr = VecWAXPY(&b,Q,U,P);CHKERRQ(ierr);       /* p <- u + b(q + b p) */
128
 
    ierr = KSP_PCApplyBAorAB(ksp,ksp->B,ksp->pc_side,P,V,Q);CHKERRQ(ierr); /* v <- K p  */
129
 
 
130
 
    rhoold = rho;
131
 
    dpold  = dp;
132
 
  }
133
 
  if (i == maxit) {
134
 
    i--;
135
 
    ksp->reason = KSP_DIVERGED_ITS;
136
 
  }
137
 
 
138
 
  ierr = KSPUnwindPreconditioner(ksp,X,T);CHKERRQ(ierr);
139
 
  *its = i + 1;
140
 
  PetscFunctionReturn(0);
141
 
}
142
 
 
143
 
EXTERN_C_BEGIN
144
 
#undef __FUNCT__  
145
 
#define __FUNCT__ "KSPCreate_TFQMR"
146
 
int KSPCreate_TFQMR(KSP ksp)
147
 
{
148
 
  PetscFunctionBegin;
149
 
  ksp->data                      = (void*)0;
150
 
  ksp->pc_side                   = PC_LEFT;
151
 
  ksp->ops->setup                = KSPSetUp_TFQMR;
152
 
  ksp->ops->solve                = KSPSolve_TFQMR;
153
 
  ksp->ops->destroy              = KSPDefaultDestroy;
154
 
  ksp->ops->buildsolution        = KSPDefaultBuildSolution;
155
 
  ksp->ops->buildresidual        = KSPDefaultBuildResidual;
156
 
  ksp->ops->setfromoptions       = 0;
157
 
  ksp->ops->view                 = 0;
158
 
  PetscFunctionReturn(0);
159
 
}
160
 
EXTERN_C_END