APBS  1.4.1
cgd.c
1 
50 #include "cgd.h"
51 
52 VPUBLIC void Vcghs(int *nx, int *ny, int *nz,
53  int *ipc, double *rpc,
54  double *ac, double *cc, double *fc,
55  double *x, double *p, double *ap, double *r,
56  int *itmax, int *iters,
57  double *errtol, double *omega,
58  int *iresid, int *iadjoint) {
59 
60  double rsnrm, pAp, denom;
61  double rhok1, rhok2, alpha, beta;
62 
63  // Setup for the looping
64  *iters = 0;
65 
66  if (*iters >= *itmax && *iresid == 0)
67  return;
68 
69  Vmresid(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
70  denom = Vxnrm2(nx, ny, nz, r);
71 
72  if (denom == 0.0)
73  return;
74 
75  if (*iters >= *itmax)
76  return;
77 
78  while(1) {
79 
80  // Compute/check the current stopping test
81  rhok2 = Vxdot(nx, ny, nz, r, r);
82  rsnrm = VSQRT(rhok2);
83 
84  if (rsnrm / denom <= *errtol)
85  break;
86 
87  if (*iters >= *itmax)
88  break;
89 
90  // Form new direction vector from old one and residual
91  if (*iters == 0) {
92  Vxcopy(nx, ny, nz, r, p);
93  } else {
94  beta = rhok2 / rhok1;
95  alpha = 1.0 / beta;
96  Vxaxpy(nx, ny, nz, &alpha, r, p);
97  Vxscal(nx, ny, nz, &beta, p);
98  }
99 
100  // Linear case: alpha which minimizes energy norm of error
101  Vmatvec(nx, ny, nz, ipc, rpc, ac, cc, p, ap);
102  pAp = Vxdot(nx, ny, nz, p, ap);
103  alpha = rhok2 / pAp;
104 
105  // Save rhok2 for next iteration
106  rhok1 = rhok2;
107 
108  // Update solution in direction p of length alpha
109  Vxaxpy(nx, ny, nz, &alpha, p, x);
110 
111  // Update residual
112  alpha = -alpha;
113  Vxaxpy(nx, ny, nz, &alpha, ap, r);
114 
115  // some bookkeeping
116  (*iters)++;
117  }
118 }
119 
VPUBLIC void Vcghs(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *p, double *ap, double *r, int *itmax, int *iters, double *errtol, double *omega, int *iresid, int *iadjoint)
A collection of useful low-level routines (timing, etc).
Definition: cgd.c:52
VPUBLIC void Vxscal(int *nx, int *ny, int *nz, double *fac, double *x)
Scale operation for a grid function with boundary values.
Definition: mikpckd.c:313
VPUBLIC void Vxcopy(int *nx, int *ny, int *nz, double *x, double *y)
A collection of useful low-level routines (timing, etc).
Definition: mikpckd.c:52
VPUBLIC double Vxdot(int *nx, int *ny, int *nz, double *x, double *y)
Inner product operation for a grid function with boundary values.
Definition: mikpckd.c:168
VPUBLIC void Vmatvec(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *x, double *y)
Matrix-vector multiplication routines.
Definition: matvecd.c:52
VPUBLIC void Vxaxpy(int *nx, int *ny, int *nz, double *alpha, double *x, double *y)
saxpy operation for a grid function with boundary values.
Definition: mikpckd.c:107
int ny
Definition: vpmgp.h:84
int nx
Definition: vpmgp.h:83
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:147
int itmax
Definition: vpmgp.h:122
VPUBLIC void Vmresid(int *nx, int *ny, int *nz, int *ipc, double *rpc, double *ac, double *cc, double *fc, double *x, double *r)
Break the matrix data-structure into diagonals and then call the residual routine.
Definition: matvecd.c:422
double errtol
Definition: vpmgp.h:121
int nz
Definition: vpmgp.h:85