← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to src/KSPACE/fix_tune_kspace.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2013-11-20 22:41:36 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20131120224136-tzx7leh606fqnckm
Tags: 0~20131119.git7162cf0-1
* [e65b919] Imported Upstream version 0~20131119.git7162cf0
* [f7bddd4] Fix some problems, introduced by upstream recently.
* [3616dfc] Use wrap-and-sort script.
* [7e92030] Ignore quilt dir

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/KSPACE/fix_tune_kspace.cpp
 
1
/* ----------------------------------------------------------------------
 
2
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
 
3
   http://lammps.sandia.gov, Sandia National Laboratories
 
4
   Steve Plimpton, sjplimp@sandia.gov
 
5
 
 
6
   Copyright (2003) Sandia Corporation.  Under the terms of Contract
 
7
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
 
8
   certain rights in this software.  This software is distributed under
 
9
   the GNU General Public License.
 
10
 
 
11
   See the README file in the top-level LAMMPS directory.
 
12
------------------------------------------------------------------------- */
 
13
 
 
14
/* ----------------------------------------------------------------------
 
15
   Contributing author: Paul Crozier (SNL)
 
16
------------------------------------------------------------------------- */
 
17
 
 
18
#include "string.h"
 
19
#include "stdlib.h"
 
20
#include "fix_tune_kspace.h"
 
21
#include "update.h"
 
22
#include "domain.h"
 
23
#include "atom.h"
 
24
#include "comm.h"
 
25
#include "force.h"
 
26
#include "kspace.h"
 
27
#include "pair.h"
 
28
#include "error.h"
 
29
#include "memory.h"
 
30
#include "timer.h"
 
31
#include "neighbor.h"
 
32
#include "modify.h"
 
33
#include "compute.h"
 
34
#include <iostream>
 
35
#include <cmath>
 
36
#include <limits>
 
37
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
 
38
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
 
39
#define GOLD 1.618034
 
40
 
 
41
using namespace std;
 
42
using namespace LAMMPS_NS;
 
43
using namespace FixConst;
 
44
 
 
45
/* ---------------------------------------------------------------------- */
 
46
 
 
47
FixTuneKspace::FixTuneKspace(LAMMPS *lmp, int narg, char **arg) :
 
48
  Fix(lmp, narg, arg)
 
49
{
 
50
  if (narg < 3) error->all(FLERR,"Illegal fix tune/kspace command");
 
51
 
 
52
  global_freq = 1;
 
53
  firststep = 0;
 
54
  niter = 0;
 
55
  niter_adjust_rcut = 0;
 
56
  keep_bracketing = true;
 
57
  first_brent_pass = true;
 
58
  converged = false;
 
59
  need_fd2_brent = false;
 
60
 
 
61
  ewald_time = pppm_time = msm_time = 0.0;
 
62
 
 
63
  // parse arguments
 
64
 
 
65
  nevery = force->inumeric(FLERR,arg[3]);
 
66
 
 
67
  // set up reneighboring
 
68
 
 
69
  force_reneighbor = 1;
 
70
  next_reneighbor = update->ntimestep + 1;
 
71
}
 
72
 
 
73
/* ---------------------------------------------------------------------- */
 
74
 
 
75
int FixTuneKspace::setmask()
 
76
{
 
77
  int mask = 0;
 
78
  mask |= PRE_EXCHANGE;
 
79
  mask |= PRE_NEIGHBOR;
 
80
  return mask;
 
81
}
 
82
 
 
83
/* ---------------------------------------------------------------------- */
 
84
 
 
85
void FixTuneKspace::init()
 
86
{
 
87
  if (!force->kspace) 
 
88
    error->all(FLERR,"Cannot use fix tune/kspace without a kspace style");
 
89
  if (!force->pair) 
 
90
    error->all(FLERR,"Cannot use fix tune/kspace without a pair style");
 
91
 
 
92
  double old_acc = force->kspace->accuracy/force->kspace->two_charge_force;
 
93
  char old_acc_str[12];
 
94
  sprintf(old_acc_str,"%g",old_acc);
 
95
  strcpy(new_acc_str,old_acc_str);
 
96
 
 
97
  int itmp;
 
98
  double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
 
99
  pair_cut_coul = *p_cutoff;
 
100
}
 
101
 
 
102
/* ----------------------------------------------------------------------
 
103
   perform dynamic kspace parameter optimization
 
104
------------------------------------------------------------------------- */
 
105
 
 
106
void FixTuneKspace::pre_exchange()
 
107
{
 
108
  if (!nevery) return;
 
109
  if (!force->kspace) return;
 
110
  if (!force->pair) return;
 
111
  if (next_reneighbor != update->ntimestep) return;
 
112
  next_reneighbor = update->ntimestep + nevery;
 
113
 
 
114
  double time = get_timing_info();
 
115
 
 
116
  if (strcmp(force->kspace_style,"ewald") == 0) ewald_time = time;
 
117
  if (strcmp(force->kspace_style,"pppm") == 0) pppm_time = time;
 
118
  if (strcmp(force->kspace_style,"msm") == 0) msm_time = time;
 
119
 
 
120
  niter++;
 
121
  if (niter == 1) {
 
122
    // test Ewald
 
123
    store_old_kspace_settings();
 
124
    strcpy(new_kspace_style,"ewald");
 
125
    sprintf(new_pair_style,"%s/long",base_pair_style);
 
126
    update_pair_style(new_pair_style,pair_cut_coul);
 
127
    update_kspace_style(new_kspace_style,new_acc_str);
 
128
  } else if (niter == 2) {
 
129
    // test PPPM
 
130
    store_old_kspace_settings();
 
131
    strcpy(new_kspace_style,"pppm");
 
132
    sprintf(new_pair_style,"%s/long",base_pair_style);
 
133
    update_pair_style(new_pair_style,pair_cut_coul);
 
134
    update_kspace_style(new_kspace_style,new_acc_str);
 
135
  } else if (niter == 3) {
 
136
    // test MSM
 
137
    store_old_kspace_settings();
 
138
    strcpy(new_kspace_style,"msm");
 
139
    sprintf(new_pair_style,"%s/msm",base_pair_style);
 
140
    update_pair_style(new_pair_style,pair_cut_coul);
 
141
    update_kspace_style(new_kspace_style,new_acc_str);
 
142
  } else if (niter == 4) {
 
143
    store_old_kspace_settings();
 
144
    cout << "ewald_time = " << ewald_time << endl;
 
145
    cout << "pppm_time = " << pppm_time << endl;
 
146
    cout << "msm_time = " << msm_time << endl;
 
147
    // switch to fastest one
 
148
    strcpy(new_kspace_style,"ewald");
 
149
    sprintf(new_pair_style,"%s/long",base_pair_style);
 
150
    if (pppm_time < ewald_time && pppm_time < msm_time)
 
151
      strcpy(new_kspace_style,"pppm");
 
152
    else if (msm_time < pppm_time && msm_time < ewald_time) {
 
153
      strcpy(new_kspace_style,"msm");
 
154
      sprintf(new_pair_style,"%s/msm",base_pair_style);
 
155
    }
 
156
    update_pair_style(new_pair_style,pair_cut_coul);
 
157
    update_kspace_style(new_kspace_style,new_acc_str);
 
158
  } else {
 
159
    adjust_rcut(time);
 
160
  }
 
161
 
 
162
  last_spcpu = timer->elapsed(TIME_LOOP);
 
163
}
 
164
 
 
165
/* ----------------------------------------------------------------------
 
166
   figure out CPU time per timestep since last time checked
 
167
------------------------------------------------------------------------- */
 
168
 
 
169
double FixTuneKspace::get_timing_info()
 
170
{
 
171
  double dvalue;
 
172
  double new_cpu;
 
173
  int new_step = update->ntimestep;
 
174
 
 
175
  if (firststep == 0) {
 
176
    new_cpu = 0.0;
 
177
    dvalue = 0.0;
 
178
    firststep = 1;
 
179
  } else {
 
180
    new_cpu = timer->elapsed(TIME_LOOP);
 
181
    double cpu_diff = new_cpu - last_spcpu;
 
182
    int step_diff = new_step - last_step;
 
183
    if (step_diff > 0.0) dvalue = cpu_diff/step_diff;
 
184
    else dvalue = 0.0;
 
185
  }
 
186
 
 
187
  last_step = new_step;
 
188
  last_spcpu = new_cpu;
 
189
 
 
190
  return dvalue;
 
191
}
 
192
 
 
193
/* ----------------------------------------------------------------------
 
194
   store old kspace settings: style, accuracy, order, etc
 
195
------------------------------------------------------------------------- */
 
196
 
 
197
void FixTuneKspace::store_old_kspace_settings()
 
198
{
 
199
  int n = strlen(force->kspace_style) + 1;
 
200
  char *old_kspace_style = new char[n];
 
201
  strcpy(old_kspace_style,force->kspace_style);
 
202
  strcpy(new_kspace_style,old_kspace_style);
 
203
  double old_acc = force->kspace->accuracy_relative;
 
204
  char old_acc_str[12];
 
205
  sprintf(old_acc_str,"%g",old_acc);
 
206
  strcpy(new_pair_style,force->pair_style);
 
207
  strcpy(base_pair_style,force->pair_style);
 
208
  char *trunc;
 
209
  if ((trunc = strstr(base_pair_style, "/long")) != NULL) *trunc = '\0';
 
210
  if ((trunc = strstr(base_pair_style, "/msm" )) != NULL) *trunc = '\0';
 
211
 
 
212
  old_differentiation_flag = force->kspace->differentiation_flag;
 
213
  old_slabflag = force->kspace->slabflag;
 
214
  old_slab_volfactor = force->kspace->slab_volfactor;
 
215
}
 
216
 
 
217
/* ----------------------------------------------------------------------
 
218
   update the pair style if necessary, preserving the settings
 
219
------------------------------------------------------------------------- */
 
220
 
 
221
void FixTuneKspace::update_pair_style(char *new_pair_style, double pair_cut_coul)
 
222
{
 
223
  int itmp;
 
224
  double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
 
225
  *p_cutoff = pair_cut_coul;
 
226
 
 
227
  // check to see if we need to change pair styles
 
228
  if (strcmp(new_pair_style,force->pair_style) == 0) return;
 
229
 
 
230
  // create a temporary file to store current pair settings
 
231
  FILE *p_pair_settings_file;
 
232
  p_pair_settings_file = tmpfile();
 
233
  force->pair->write_restart(p_pair_settings_file);
 
234
  rewind(p_pair_settings_file);
 
235
 
 
236
  cout << "Creating new pair style: " << new_pair_style << endl;
 
237
  // delete old pair style and create new one
 
238
  force->create_pair(new_pair_style,lmp->suffix);
 
239
 
 
240
  // restore current pair settings from temporary file
 
241
  force->pair->read_restart(p_pair_settings_file);
 
242
 
 
243
  double *pcutoff = (double *) force->pair->extract("cut_coul",itmp);
 
244
  double current_cutoff = *pcutoff;
 
245
  cout << "Coulomb cutoff for real space: " << current_cutoff << endl;
 
246
 
 
247
  // close temporary file
 
248
  fclose(p_pair_settings_file);
 
249
}
 
250
 
 
251
/* ----------------------------------------------------------------------
 
252
   update the kspace style if necessary
 
253
------------------------------------------------------------------------- */
 
254
 
 
255
void FixTuneKspace::update_kspace_style(char *new_kspace_style, char *new_acc_str)
 
256
{
 
257
  // create kspace style char string
 
258
 
 
259
  int narg = 2;
 
260
  char **arg;
 
261
  arg = NULL;
 
262
  int maxarg = 100;
 
263
  arg = (char **) memory->srealloc(arg,maxarg*sizeof(char *),"tune/kspace:arg");
 
264
  int n = 12;
 
265
  arg[0] = new char[n];
 
266
  strcpy(arg[0],new_kspace_style);
 
267
  arg[1] = new char[n];
 
268
  strcpy(arg[1],new_acc_str);
 
269
 
 
270
  // delete old kspace style and create new one
 
271
 
 
272
  force->create_kspace(narg,arg,lmp->suffix);
 
273
 
 
274
  force->kspace->differentiation_flag = old_differentiation_flag;
 
275
  force->kspace->slabflag = old_slabflag;
 
276
  force->kspace->slab_volfactor = old_slab_volfactor;
 
277
 
 
278
  // initialize new kspace style, pair style, molecular styles
 
279
 
 
280
  force->init();
 
281
 
 
282
  // set up grid
 
283
  force->kspace->setup_grid();
 
284
 
 
285
  // Re-init neighbor list. Probably only needed when redefining the pair style. Should happen after pair->init() to get pair style neighbor list request registered
 
286
 
 
287
  neighbor->init();
 
288
 
 
289
  // Re-init computes to update pointers to virials, etc.
 
290
 
 
291
  for (int i = 0; i < modify->ncompute; i++) modify->compute[i]->init();
 
292
 
 
293
  memory->sfree(arg);
 
294
}
 
295
 
 
296
/* ----------------------------------------------------------------------
 
297
   find the optimal real space coulomb cutoff
 
298
------------------------------------------------------------------------- */
 
299
 
 
300
void FixTuneKspace::adjust_rcut(double time)
 
301
{
 
302
  if (strcmp(force->kspace_style,"msm") == 0) return;
 
303
  if (converged) return;
 
304
 
 
305
  double temp;
 
306
  const double TINY = 1.0e-20;
 
307
 
 
308
  // get the current cutoff
 
309
  int itmp;
 
310
  double *p_cutoff = (double *) force->pair->extract("cut_coul",itmp);
 
311
  double current_cutoff = *p_cutoff;
 
312
  cout << "Old Coulomb cutoff for real space: " << current_cutoff << endl;
 
313
 
 
314
  // use Brent's method from Numerical Recipes to find optimal real space cutoff
 
315
 
 
316
  // first time through, get ax_brent and fa_brent, and adjust cutoff
 
317
  if (keep_bracketing) {
 
318
    if (niter_adjust_rcut == 0) {
 
319
      pair_cut_coul /= 2;
 
320
    } else if (niter_adjust_rcut == 1) {
 
321
      ax_brent = current_cutoff;
 
322
      fa_brent = time;
 
323
      pair_cut_coul *= 2;
 
324
 
 
325
    // second time through, get bx_brent and fb_brent, and adjust cutoff
 
326
    } else if (niter_adjust_rcut == 2) {
 
327
      bx_brent = current_cutoff;
 
328
      fb_brent = time;
 
329
      if (fb_brent > fa_brent) {
 
330
        SWAP(ax_brent,bx_brent);
 
331
        SWAP(fb_brent,fa_brent);
 
332
        pair_cut_coul /= 4;
 
333
      } else {
 
334
        pair_cut_coul *= 2;
 
335
      }
 
336
 
 
337
    // third time through, get cx_brent and fc_brent, and adjust cutoff if needed
 
338
    } else if (niter_adjust_rcut == 3) {
 
339
      cx_brent = current_cutoff;
 
340
      fc_brent = time;
 
341
      if (fc_brent > fb_brent) keep_bracketing = false;
 
342
      else {
 
343
        double r = (bx_brent - ax_brent)*(fb_brent - fc_brent);
 
344
        double q = (bx_brent - cx_brent)*(fb_brent - fa_brent);
 
345
        dx_brent = bx_brent - ((bx_brent - cx_brent)*q - (bx_brent - ax_brent)*r)/
 
346
         (2.0*SIGN(MAX(fabs(q - r),TINY),q - r));
 
347
        pair_cut_coul = dx_brent;
 
348
      }
 
349
 
 
350
    // after third time through, bracket the minimum, and adjust cutoff
 
351
    } else if (niter_adjust_rcut > 3) {
 
352
      dx_brent = current_cutoff;
 
353
      if (need_fd2_brent) fd2_brent = time;
 
354
      else fd_brent = time;
 
355
      mnbrak();
 
356
      pair_cut_coul = dx_brent;
 
357
    }
 
358
  }
 
359
 
 
360
  if (!keep_bracketing) {
 
361
    dx_brent = current_cutoff;
 
362
    fd_brent = time;
 
363
    if (first_brent_pass) brent0();
 
364
    else brent2();
 
365
    brent1();
 
366
    pair_cut_coul = dx_brent;
 
367
  }
 
368
 
 
369
  niter_adjust_rcut++;
 
370
 
 
371
  if (pair_cut_coul <= 0.0) pair_cut_coul = fabs(MIN(ax_brent,MIN(bx_brent,(MIN(cx_brent,dx_brent))))/2.0) + TINY;
 
372
 
 
373
  if (pair_cut_coul != pair_cut_coul)
 
374
    error->all(FLERR,"Bad real space Coulomb cutoff in fix tune/kspace");
 
375
 
 
376
  // change the cutoff to pair_cut_coul
 
377
  *p_cutoff = pair_cut_coul;
 
378
 
 
379
  // report the new cutoff
 
380
  double *new_cutoff = (double *) force->pair->extract("cut_coul",itmp);
 
381
  current_cutoff = *new_cutoff;
 
382
  cout << "Adjusted Coulomb cutoff for real space: " << current_cutoff << endl;
 
383
 
 
384
  store_old_kspace_settings();
 
385
  update_pair_style(new_pair_style,pair_cut_coul);
 
386
  update_kspace_style(new_kspace_style,new_acc_str);
 
387
}
 
388
 
 
389
/* ----------------------------------------------------------------------
 
390
   bracket a minimum using parabolic extrapolation
 
391
------------------------------------------------------------------------- */
 
392
 
 
393
void FixTuneKspace::mnbrak()
 
394
{
 
395
  const double GLIMIT = 100.0, TINY = 1.0e-20;
 
396
  double temp,r,q;
 
397
  r = (bx_brent - ax_brent)*(fb_brent - fc_brent);
 
398
  q = (bx_brent - cx_brent)*(fb_brent - fa_brent);
 
399
  dx_brent = bx_brent - ((bx_brent - cx_brent)*q - (bx_brent - ax_brent)*r)/
 
400
   (2.0*SIGN(MAX(fabs(q - r),TINY),q - r));
 
401
  dxlim = bx_brent + GLIMIT*(cx_brent - bx_brent);
 
402
 
 
403
  if ((bx_brent - dx_brent)*(dx_brent - cx_brent) > 0.0) {
 
404
    if (fd_brent < fc_brent) {
 
405
      ax_brent = bx_brent;
 
406
      bx_brent = dx_brent;
 
407
      fa_brent = fb_brent;
 
408
      fb_brent = fd_brent;
 
409
      keep_bracketing = false;
 
410
      return;
 
411
    } else if (fd_brent > fb_brent) {
 
412
      cx_brent = dx_brent;
 
413
      fc_brent = fd_brent;
 
414
      keep_bracketing = false;
 
415
      return;
 
416
    }
 
417
    dx_brent = cx_brent + GOLD*(cx_brent - bx_brent);
 
418
    if (need_fd2_brent) {
 
419
      fd_brent = fd2_brent;
 
420
      need_fd2_brent = false;
 
421
    } else {
 
422
      need_fd2_brent = true;
 
423
      return;
 
424
    }
 
425
  } else if ((cx_brent - dx_brent)*(dx_brent - dxlim) > 0.0) {
 
426
    if (fd_brent < fc_brent) {
 
427
      if (need_fd2_brent) {
 
428
        need_fd2_brent = false;
 
429
      } else {
 
430
        need_fd2_brent = true;
 
431
        dx_brent += GOLD*(dx_brent - cx_brent);
 
432
        return;
 
433
      }
 
434
      shft3(bx_brent,cx_brent,dx_brent,dx_brent + GOLD*(dx_brent - cx_brent));
 
435
      shft3(fb_brent,fc_brent,fd_brent,fd2_brent);
 
436
    }
 
437
  } else if ((dx_brent - dxlim)*(dxlim - cx_brent) >= 0.0) {
 
438
    dx_brent = dxlim;
 
439
    if (need_fd2_brent) {
 
440
      fd_brent = fd2_brent;
 
441
      need_fd2_brent = false;
 
442
    } else {
 
443
      need_fd2_brent = true;
 
444
      return;
 
445
    }
 
446
  } else {
 
447
    dx_brent = cx_brent + GOLD*(cx_brent - bx_brent);
 
448
    if (need_fd2_brent) {
 
449
      fd_brent = fd2_brent;
 
450
      need_fd2_brent = false;
 
451
    } else {
 
452
      need_fd2_brent = true;
 
453
      return;
 
454
    }
 
455
  }
 
456
  shft3(ax_brent,bx_brent,cx_brent,dx_brent);
 
457
  shft3(fa_brent,fb_brent,fc_brent,fd_brent);
 
458
}
 
459
 
 
460
/* ----------------------------------------------------------------------
 
461
   Brent's method from Numerical Recipes
 
462
------------------------------------------------------------------------- */
 
463
 
 
464
void FixTuneKspace::brent0()
 
465
{
 
466
  a_brent=(ax_brent < cx_brent ? ax_brent : cx_brent);
 
467
  b_brent=(ax_brent > cx_brent ? ax_brent : cx_brent);
 
468
  x_brent=w_brent=v_brent=bx_brent;
 
469
  fw_brent=fv_brent=fx_brent=fb_brent;
 
470
}
 
471
 
 
472
/* ----------------------------------------------------------------------
 
473
   Brent's method from Numerical Recipes
 
474
------------------------------------------------------------------------- */
 
475
 
 
476
void FixTuneKspace::brent1()
 
477
{
 
478
  const int ITMAX=100;
 
479
  const double CGOLD=0.3819660;
 
480
  const double ZEPS=numeric_limits<double>::epsilon()*1.0e-3;
 
481
  double d=0.0,etemp;
 
482
  double p,q,r,tol1,tol2,xm;
 
483
  double e=0.0;
 
484
  double tol=0.001;
 
485
 
 
486
  xm=0.5*(a_brent+b_brent);
 
487
  tol2=2.0*(tol1=tol*fabs(x_brent)+ZEPS);
 
488
  if (fabs(x_brent-xm) <= (tol2-0.5*(b_brent-a_brent))) {
 
489
    converged = true;
 
490
    dx_brent = x_brent;
 
491
    return;
 
492
  }
 
493
  if (fabs(e) > tol1) {
 
494
    r=(x_brent-w_brent)*(fx_brent-fv_brent);
 
495
    q=(x_brent-v_brent)*(fx_brent-fw_brent);
 
496
    p=(x_brent-v_brent)*q-(x_brent-w_brent)*r;
 
497
    q=2.0*(q-r);
 
498
    if (q > 0.0) p = -p;
 
499
    q=fabs(q);
 
500
    etemp=e;
 
501
    e=d;
 
502
    if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a_brent-x_brent) || p >= q*(b_brent-x_brent))
 
503
      d=CGOLD*(e=(x_brent >= xm ? a_brent-x_brent : b_brent-x_brent));
 
504
    else {
 
505
      d=p/q;
 
506
      dx_brent=x_brent+d;
 
507
      if (dx_brent-a_brent < tol2 || b_brent-dx_brent < tol2)
 
508
        d=SIGN(tol1,xm-x_brent);
 
509
    }
 
510
  } else {
 
511
    d=CGOLD*(e=(x_brent >= xm ? a_brent-x_brent : b_brent-x_brent));
 
512
  }
 
513
  dx_brent=(fabs(d) >= tol1 ? x_brent+d : x_brent+SIGN(tol1,d));
 
514
 
 
515
  first_brent_pass = false;
 
516
 
 
517
  return;
 
518
}
 
519
 
 
520
/* ----------------------------------------------------------------------
 
521
   Brent's method from Numerical Recipes
 
522
------------------------------------------------------------------------- */
 
523
 
 
524
void FixTuneKspace::brent2()
 
525
{
 
526
  if (fd_brent <= fx_brent) {
 
527
    if (dx_brent >= x_brent) a_brent=x_brent; else b_brent=x_brent;
 
528
    shft3(v_brent,w_brent,x_brent,dx_brent);
 
529
    shft3(fv_brent,fw_brent,fx_brent,fd_brent);
 
530
  } else {
 
531
    if (dx_brent < x_brent) a_brent=dx_brent; else b_brent=dx_brent;
 
532
    if (fd_brent <= fw_brent || w_brent == x_brent) {
 
533
      v_brent=w_brent;
 
534
      w_brent=dx_brent;
 
535
      fv_brent=fw_brent;
 
536
      fw_brent=fd_brent;
 
537
    } else if (fd_brent <= fv_brent || v_brent == x_brent || v_brent == w_brent) {
 
538
      v_brent=dx_brent;
 
539
      fv_brent=fd_brent;
 
540
    }
 
541
  }
 
542
}
 
543