APBS 3.0.0
Loading...
Searching...
No Matches
cgd.c
1
55#include "cgd.h"
56
57VPUBLIC 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