APBS  3.0.0
cgd.c
1 
55 #include "cgd.h"
56 
57 VPUBLIC void Vcghs(int *nx, int *ny, int *nz,
58  int *ipc, double *rpc,
59  double *ac, double *cc, double *fc,
60  double *x, double *p, double *ap, double *r,
61  int *itmax, int *iters,
62  double *errtol, double *omega,
63  int *iresid, int *iadjoint) {
64 
65  double rsnrm, pAp, denom;
66  double rhok1, rhok2, alpha, beta;
67 
68  // Setup for the looping
69  *iters = 0;
70 
71  if (*iters >= *itmax && *iresid == 0)
72  return;
73 
74  Vmresid(nx, ny, nz, ipc, rpc, ac, cc, fc, x, r);
75  denom = Vxnrm2(nx, ny, nz, r);
76 
77  if (denom == 0.0)
78  return;
79 
80  if (*iters >= *itmax)
81  return;
82 
83  while(1) {
84 
85  // Compute/check the current stopping test
86  rhok2 = Vxdot(nx, ny, nz, r, r);
87  rsnrm = VSQRT(rhok2);
88 
89  if (rsnrm / denom <= *errtol)
90  break;
91 
92  if (*iters >= *itmax)
93  break;
94 
95  // Form new direction vector from old one and residual
96  if (*iters == 0) {
97  Vxcopy(nx, ny, nz, r, p);
98  } else {
99  beta = rhok2 / rhok1;
100  alpha = 1.0 / beta;
101  Vxaxpy(nx, ny, nz, &alpha, r, p);
102  Vxscal(nx, ny, nz, &beta, p);
103  }
104 
105  // Linear case: alpha which minimizes energy norm of error
106  Vmatvec(nx, ny, nz, ipc, rpc, ac, cc, p, ap);
107  pAp = Vxdot(nx, ny, nz, p, ap);
108  alpha = rhok2 / pAp;
109 
110  // Save rhok2 for next iteration
111  rhok1 = rhok2;
112 
113  // Update solution in direction p of length alpha
114  Vxaxpy(nx, ny, nz, &alpha, p, x);
115 
116  // Update residual
117  alpha = -alpha;
118  Vxaxpy(nx, ny, nz, &alpha, ap, r);
119 
120  // some bookkeeping
121  (*iters)++;
122  }
123 }
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:57
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:173
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:318
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:57
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:57
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:426
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:112
VPUBLIC double Vxnrm2(int *nx, int *ny, int *nz, double *x)
Norm operation for a grid function with boundary values.
Definition: mikpckd.c:152