Peano
Newton.cpp
Go to the documentation of this file.
1 /* TwoPunctures: File "Newton.c"*/
2 
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <math.h>
7 #include <ctype.h>
8 #include <gsl/gsl_vector.h>
9 #include <gsl/gsl_linalg.h>
10 #include "libtwopunctures/TwoPunctures.h"
11 
12 namespace TP {
13 using namespace Utilities;
14 
15 const int StencilSize = 19;
16 const int N_PlaneRelax = 1;
17 const int NRELAX = 200;
18 const int Step_Relax = 1;
19 
20 /* --------------------------------------------------------------------------*/
21 double
22 TwoPunctures::norm_inf (double const * TP_RESTRICT const F,
23  int const ntotal)
24 {
25  double dmax = -1;
26  {
27  double dmax1 = -1;
28  for (int j = 0; j < ntotal; j++)
29  if (fabs (F[j]) > dmax1)
30  dmax1 = fabs (F[j]);
31  if (dmax1 > dmax)
32  dmax = dmax1;
33  }
34  return dmax;
35 }
36 /* --------------------------------------------------------------------------*/
37 void
38 TwoPunctures::resid (double * TP_RESTRICT const res,
39  int const ntotal,
40  double const * TP_RESTRICT const dv,
41  double const * TP_RESTRICT const rhs,
42  int const * TP_RESTRICT const ncols,
43  int const * TP_RESTRICT const * TP_RESTRICT const cols,
44  double const * TP_RESTRICT const * TP_RESTRICT const JFD)
45 {
46  for (int i = 0; i < ntotal; i++)
47  {
48  double JFDdv_i = 0;
49  for (int m = 0; m < ncols[i]; m++)
50  JFDdv_i += JFD[i][m] * dv[cols[i][m]];
51  res[i] = rhs[i] - JFDdv_i;
52  }
53 }
54 
55 /* -------------------------------------------------------------------------*/
56 void
57 TwoPunctures::LineRelax_al (double * TP_RESTRICT const dv,
58  int const j, int const k, int const nvar,
59  int const n1, int const n2, int const n3,
60  double const * TP_RESTRICT const rhs,
61  int const * TP_RESTRICT const ncols,
62  int const * TP_RESTRICT const * TP_RESTRICT const cols,
63  double const * TP_RESTRICT const * TP_RESTRICT const JFD)
64 {
65  int i, m, Ic, Ip, Im, col, ivar;
66 
67  gsl_vector *diag = gsl_vector_alloc(n1);
68  gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
69  gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
70  gsl_vector *b = gsl_vector_alloc(n1); /* rhs */
71  gsl_vector *x = gsl_vector_alloc(n1); /* solution vector */
72 
73  for (ivar = 0; ivar < nvar; ivar++)
74  {
75  gsl_vector_set_zero(diag);
76  gsl_vector_set_zero(e);
77  gsl_vector_set_zero(f);
78  for (i = 0; i < n1; i++)
79  {
80  Ip = Index (ivar, i + 1, j, k, nvar, n1, n2, n3);
81  Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
82  Im = Index (ivar, i - 1, j, k, nvar, n1, n2, n3);
83  gsl_vector_set(b,i,rhs[Ic]);
84  for (m = 0; m < ncols[Ic]; m++)
85  {
86  col = cols[Ic][m];
87  if (col != Ip && col != Ic && col != Im)
88  *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
89  else
90  {
91  if (col == Im && i > 0)
92  gsl_vector_set(f,i-1,JFD[Ic][m]);
93  if (col == Ic)
94  gsl_vector_set(diag,i,JFD[Ic][m]);
95  if (col == Ip && i < n1-1)
96  gsl_vector_set(e,i,JFD[Ic][m]);
97  }
98  }
99  }
100  gsl_linalg_solve_tridiag(diag, e, f, b, x);
101  for (i = 0; i < n1; i++)
102  {
103  Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
104  dv[Ic] = gsl_vector_get(x, i);
105  }
106  }
107 
108  gsl_vector_free(diag);
109  gsl_vector_free(e);
110  gsl_vector_free(f);
111  gsl_vector_free(b);
112  gsl_vector_free(x);
113 }
114 
115 /* --------------------------------------------------------------------------*/
116 void
117 TwoPunctures::LineRelax_be (double * TP_RESTRICT const dv,
118  int const i, int const k, int const nvar,
119  int const n1, int const n2, int const n3,
120  double const * TP_RESTRICT const rhs,
121  int const * TP_RESTRICT const ncols,
122  int const * TP_RESTRICT const * TP_RESTRICT const cols,
123  double const * TP_RESTRICT const * TP_RESTRICT const JFD)
124 {
125  int j, m, Ic, Ip, Im, col, ivar;
126 
127  gsl_vector *diag = gsl_vector_alloc(n2);
128  gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
129  gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
130  gsl_vector *b = gsl_vector_alloc(n2); /* rhs */
131  gsl_vector *x = gsl_vector_alloc(n2); /* solution vector */
132 
133  for (ivar = 0; ivar < nvar; ivar++)
134  {
135  gsl_vector_set_zero(diag);
136  gsl_vector_set_zero(e);
137  gsl_vector_set_zero(f);
138  for (j = 0; j < n2; j++)
139  {
140  Ip = Index (ivar, i, j + 1, k, nvar, n1, n2, n3);
141  Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
142  Im = Index (ivar, i, j - 1, k, nvar, n1, n2, n3);
143  gsl_vector_set(b,j,rhs[Ic]);
144  for (m = 0; m < ncols[Ic]; m++)
145  {
146  col = cols[Ic][m];
147  if (col != Ip && col != Ic && col != Im)
148  *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
149  else
150  {
151  if (col == Im && j > 0)
152  gsl_vector_set(f,j-1,JFD[Ic][m]);
153  if (col == Ic)
154  gsl_vector_set(diag,j,JFD[Ic][m]);
155  if (col == Ip && j < n2-1)
156  gsl_vector_set(e,j,JFD[Ic][m]);
157  }
158  }
159  }
160  gsl_linalg_solve_tridiag(diag, e, f, b, x);
161  for (j = 0; j < n2; j++)
162  {
163  Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
164  dv[Ic] = gsl_vector_get(x, j);
165  }
166  }
167  gsl_vector_free(diag);
168  gsl_vector_free(e);
169  gsl_vector_free(f);
170  gsl_vector_free(b);
171  gsl_vector_free(x);
172 }
173 
174 /* --------------------------------------------------------------------------*/
175 void
176 TwoPunctures::relax (double * TP_RESTRICT const dv,
177  int const nvar, int const n1, int const n2, int const n3,
178  double const * TP_RESTRICT const rhs,
179  int const * TP_RESTRICT const ncols,
180  int const * TP_RESTRICT const * TP_RESTRICT const cols,
181  double const * TP_RESTRICT const * TP_RESTRICT const JFD)
182 {
183  int i, j, k, n;
184 
185  for (k = 0; k < n3; k = k + 2)
186  {
187  for (n = 0; n < N_PlaneRelax; n++)
188  {
189 
190  for (i = 2; i < n1; i = i + 2)
191  LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
192 
193  for (i = 1; i < n1; i = i + 2)
194  LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
195 
196  for (j = 1; j < n2; j = j + 2)
197  LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
198 
199  for (j = 0; j < n2; j = j + 2)
200  LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
201  }
202  }
203  for (k = 1; k < n3; k = k + 2)
204  {
205  for (n = 0; n < N_PlaneRelax; n++)
206  {
207 
208  for (i = 0; i < n1; i = i + 2)
209  LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
210 
211  for (i = 1; i < n1; i = i + 2)
212  LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
213 
214  for (j = 1; j < n2; j = j + 2)
215  LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
216 
217  for (j = 0; j < n2; j = j + 2)
218  LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
219  }
220  }
221 }
222 
223 /* --------------------------------------------------------------------------*/
224 void
225 TwoPunctures::TestRelax (int nvar, int n1, int n2, int n3, derivs v,
226  double *dv)
227 {
228 
229  int ntotal = n1 * n2 * n3 * nvar, **cols, *ncols, maxcol =
230  StencilSize * nvar, j;
231  double *F, *res, **JFD;
232  derivs u;
233 
234  F = dvector (0, ntotal - 1);
235  res = dvector (0, ntotal - 1);
236  allocate_derivs (&u, ntotal);
237 
238  JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
239  cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
240  ncols = ivector (0, ntotal - 1);
241 
242  F_of_v (nvar, n1, n2, n3, v, F, u);
243 
244  SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
245 
246  for (j = 0; j < ntotal; j++)
247  dv[j] = 0;
248  resid (res, ntotal, dv, F, ncols, cols, JFD);
249  printf ("Before: |F|=%20.15e\n", (double) norm1 (res, ntotal));
250  fflush(stdout);
251  for (j = 0; j < NRELAX; j++)
252  {
253  relax (dv, nvar, n1, n2, n3, F, ncols, cols, JFD); /* solves JFD*sh = s*/
254  if (j % Step_Relax == 0)
255  {
256  resid (res, ntotal, dv, F, ncols, cols, JFD);
257  printf ("j=%d\t |F|=%20.15e\n", j, (double) norm1 (res, ntotal));
258  fflush(stdout);
259  }
260  }
261 
262  resid (res, ntotal, dv, F, ncols, cols, JFD);
263  printf ("After: |F|=%20.15e\n", (double) norm1 (res, ntotal));
264  fflush(stdout);
265 
266  free_dvector (F, 0, ntotal - 1);
267  free_dvector (res, 0, ntotal - 1);
268  free_derivs (&u, ntotal);
269 
270  free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
271  free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
272  free_ivector (ncols, 0, ntotal - 1);
273 }
274 
275 /* --------------------------------------------------------------------------*/
276 int
277 TwoPunctures::bicgstab (int const nvar, int const n1, int const n2, int const n3,
278  derivs v, derivs dv,
279  int const output, int const itmax, double const tol,
280  double * TP_RESTRICT const normres)
281 {
282 
283  int ntotal = n1 * n2 * n3 * nvar, ii;
284  double alpha = 0, beta = 0;
285  double rho = 0, rho1 = 1, rhotol = 1e-50;
286  double omega = 0, omegatol = 1e-50;
287  double *p, *rt, *s, *t, *r, *vv;
288  double **JFD;
289  int **cols, *ncols, maxcol = StencilSize * nvar;
290  double *F;
291  derivs u, ph, sh;
292 
293  F = dvector (0, ntotal - 1);
294  allocate_derivs (&u, ntotal);
295 
296  JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
297  cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
298  ncols = ivector (0, ntotal - 1);
299 
300  F_of_v (nvar, n1, n2, n3, v, F, u);
301  SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
302 
303  /* temporary storage */
304  r = dvector (0, ntotal - 1);
305  p = dvector (0, ntotal - 1);
306  allocate_derivs (&ph, ntotal);
307 /* ph = dvector(0, ntotal-1);*/
308  rt = dvector (0, ntotal - 1);
309  s = dvector (0, ntotal - 1);
310  allocate_derivs (&sh, ntotal);
311 /* sh = dvector(0, ntotal-1);*/
312  t = dvector (0, ntotal - 1);
313  vv = dvector (0, ntotal - 1);
314 
315  /* check */
316  if (output == 1) {
317  printf ("bicgstab: itmax %d, tol %e\n", itmax, (double)tol);
318  fflush(stdout);
319  }
320 
321  /* compute initial residual rt = r = F - J*dv */
322  J_times_dv (nvar, n1, n2, n3, dv, r, u);
323  for (int j = 0; j < ntotal; j++)
324  rt[j] = r[j] = F[j] - r[j];
325 
326  *normres = norm2 (r, ntotal);
327  if (output == 1) {
328  printf ("bicgstab: %5d %10.3e\n", 0, (double) *normres);
329  fflush(stdout);
330  }
331 
332  if (*normres <= tol)
333  return 0;
334 
335  /* cgs iteration */
336  for (ii = 0; ii < itmax; ii++)
337  {
338  rho = scalarproduct (rt, r, ntotal);
339  if (fabs (rho) < rhotol)
340  break;
341 
342  /* compute direction vector p */
343  if (ii == 0)
344  {
345 
346  for (int j = 0; j < ntotal; j++)
347  p[j] = r[j];
348  }
349  else
350  {
351  beta = (rho / rho1) * (alpha / omega);
352 
353  for (int j = 0; j < ntotal; j++)
354  p[j] = r[j] + beta * (p[j] - omega * vv[j]);
355  }
356 
357  /* compute direction adjusting vector ph and scalar alpha */
358 
359  for (int j = 0; j < ntotal; j++)
360  ph.d0[j] = 0;
361  for (int j = 0; j < NRELAX; j++) /* solves JFD*ph = p by relaxation*/
362  relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
363 
364  J_times_dv (nvar, n1, n2, n3, ph, vv, u); /* vv=J*ph*/
365  alpha = rho / scalarproduct (rt, vv, ntotal);
366 
367  for (int j = 0; j < ntotal; j++)
368  s[j] = r[j] - alpha * vv[j];
369 
370  /* early check of tolerance */
371  *normres = norm2 (s, ntotal);
372  if (*normres <= tol)
373  {
374 
375  for (int j = 0; j < ntotal; j++)
376  dv.d0[j] += alpha * ph.d0[j];
377  if (output == 1) {
378  printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
379  ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
380  fflush(stdout);
381  }
382  break;
383  }
384 
385  /* compute stabilizer vector sh and scalar omega */
386 
387  for (int j = 0; j < ntotal; j++)
388  sh.d0[j] = 0;
389  for (int j = 0; j < NRELAX; j++) /* solves JFD*sh = s by relaxation*/
390  relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
391 
392  J_times_dv (nvar, n1, n2, n3, sh, t, u); /* t=J*sh*/
393  omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
394 
395  /* compute new solution approximation */
396 
397  for (int j = 0; j < ntotal; j++)
398  {
399  dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
400  r[j] = s[j] - omega * t[j];
401  }
402  /* are we done? */
403  *normres = norm2 (r, ntotal);
404  if (output == 1) {
405  printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
406  ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
407  fflush(stdout);
408  }
409  if (*normres <= tol)
410  break;
411  rho1 = rho;
412  if (fabs (omega) < omegatol)
413  break;
414 
415  }
416 
417  /* free temporary storage */
418  free_dvector (r, 0, ntotal - 1);
419  free_dvector (p, 0, ntotal - 1);
420 /* free_dvector(ph, 0, ntotal-1);*/
421  free_derivs (&ph, ntotal);
422  free_dvector (rt, 0, ntotal - 1);
423  free_dvector (s, 0, ntotal - 1);
424 /* free_dvector(sh, 0, ntotal-1);*/
425  free_derivs (&sh, ntotal);
426  free_dvector (t, 0, ntotal - 1);
427  free_dvector (vv, 0, ntotal - 1);
428 
429  free_dvector (F, 0, ntotal - 1);
430  free_derivs (&u, ntotal);
431 
432  free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
433  free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
434  free_ivector (ncols, 0, ntotal - 1);
435 
436  /* iteration failed */
437  if (ii > itmax)
438  return -1;
439 
440  /* breakdown */
441  if (fabs (rho) < rhotol)
442  return -10;
443  if (fabs (omega) < omegatol)
444  return -11;
445 
446  /* success! */
447  return ii + 1;
448 }
449 
450 /* -------------------------------------------------------------------*/
451 void
452 TwoPunctures::Newton (int const nvar, int const n1, int const n2, int const n3,
453  derivs v,
454  double const tol, int const itmax)
455 {
456 
457 
458  int ntotal = n1 * n2 * n3 * nvar, ii, it;
459  double *F, dmax, normres;
460  derivs u, dv;
461 
462  F = dvector (0, ntotal - 1);
463  allocate_derivs (&dv, ntotal);
464  allocate_derivs (&u, ntotal);
465 
466 /* TestRelax(nvar, n1, n2, n3, v, dv.d0); */
467  it = 0;
468  dmax = 1;
469  while (dmax > tol && it < itmax)
470  {
471  if (it == 0)
472  {
473  F_of_v (nvar, n1, n2, n3, v, F, u);
474  dmax = norm_inf (F, ntotal);
475  }
476 
477  for (int j = 0; j < ntotal; j++)
478  dv.d0[j] = 0;
479 
480  if(verbose==1){
481  printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
482  printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
483  }
484 
485  fflush(stdout);
486  ii =
487  bicgstab (nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres);
488 
489  for (int j = 0; j < ntotal; j++)
490  v.d0[j] -= dv.d0[j];
491  F_of_v (nvar, n1, n2, n3, v, F, u);
492  dmax = norm_inf (F, ntotal);
493  it += 1;
494  }
495  if (itmax==0)
496  {
497  F_of_v (nvar, n1, n2, n3, v, F, u);
498  dmax = norm_inf (F, ntotal);
499  }
500 
501  if(verbose==1)
502  printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
503 
504  fflush(stdout);
505 
506  free_dvector (F, 0, ntotal - 1);
507  free_derivs (&dv, ntotal);
508  free_derivs (&u, ntotal);
509 }
510 
511 /* -------------------------------------------------------------------*/
512 
513 } // namespace TP
TP::TwoPunctures::bicgstab
int bicgstab(int const nvar, int const n1, int const n2, int const n3, derivs v, derivs dv, int const output, int const itmax, double const tol, double *TP_RESTRICT const normres)
Definition: Newton.cpp:277
TP::TwoPunctures::Newton
void Newton(int nvar, int n1, int n2, int n3, derivs v, double tol, int itmax)
Definition: Newton.cpp:452
TP::TwoPunctures::relax
void relax(double *TP_RESTRICT const dv, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition: Newton.cpp:176
TP::TwoPunctures::LineRelax_be
void LineRelax_be(double *TP_RESTRICT const dv, int const i, int const k, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition: Newton.cpp:117
TP::TwoPunctures::resid
void resid(double *TP_RESTRICT const res, int const ntotal, double const *TP_RESTRICT const dv, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition: Newton.cpp:38
TP::TwoPunctures::LineRelax_al
void LineRelax_al(double *TP_RESTRICT const dv, int const j, int const k, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition: Newton.cpp:57
TP::TwoPunctures::norm_inf
double norm_inf(double const *TP_RESTRICT const F, int const ntotal)
Definition: Newton.cpp:22
TP::TwoPunctures::TestRelax
void TestRelax(int nvar, int n1, int n2, int n3, derivs v, double *dv)
Definition: Newton.cpp:225
TP::Utilities::norm2
double norm2(double *v, int n)
Definition: TP_Utilities.cpp:506
TP::Utilities::dvector
double * dvector(long nl, long nh)
Definition: TP_Utilities.cpp:30
TP::Utilities::free_dvector
void free_dvector(double *v, long nl, long nh)
Definition: TP_Utilities.cpp:157
TP::Utilities::norm1
double norm1(double *v, int n)
Definition: TP_Utilities.cpp:492
TP::Utilities::scalarproduct
double scalarproduct(double *v, double *w, int n)
Definition: TP_Utilities.cpp:519
TP::Utilities::free_ivector
void free_ivector(int *v, long nl, long nh)
Definition: TP_Utilities.cpp:149
TP::Utilities::ivector
int * ivector(long nl, long nh)
Definition: TP_Utilities.cpp:16
TP::Utilities::dmatrix
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
Definition: TP_Utilities.cpp:73
TP::Utilities::free_dmatrix
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
Definition: TP_Utilities.cpp:174
TP::Utilities::free_imatrix
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
Definition: TP_Utilities.cpp:165
TP::Utilities::imatrix
int ** imatrix(long nrl, long nrh, long ncl, long nch)
Definition: TP_Utilities.cpp:44
TP
This file contains aliases for making access to the long state vector Q as used eg.
Definition: CoordTransf.cpp:13
TP::Step_Relax
const int Step_Relax
Definition: Newton.cpp:18
TP::StencilSize
const int StencilSize
Definition: Newton.cpp:15
TP::N_PlaneRelax
const int N_PlaneRelax
Definition: Newton.cpp:16
TP::NRELAX
const int NRELAX
Definition: Newton.cpp:17
euler.u
u
Definition: euler.py:113
euler.j
j
Definition: euler.py:95
swe.b
b
Definition: swe.py:82
TP::DERIVS::d0
double * d0
Definition: TwoPunctures.h:31