← Back to branch summary

~ubuntu-branches/ubuntu/wily/afnix/wily

« back to all changes in this revision

Viewing changes to src/mod/mth/shl/Newton.cpp

  • Committer: Package Import Robot
  • Author(s): Nobuhiro Iwamatsu
  • Date: 2015-07-11 02:00:35 UTC
  • mfrom: (10.1.1 sid)
  • Revision ID: package-import@ubuntu.com-20150711020035-2nhpztq7s15qyc0v
Tags: 2.5.1-1
http://bugs.debian.org/789968
http://bugs.debian.org/741508
http://bugs.debian.org/748146
http://bugs.debian.org/777767
* New upstream release. (Closes: #789968)
* Update debian/control.
  - Update Standards-Version to 3.9.6.
* Add support mips64(el) and ppc64el. (Closes: #741508, #748146)
* Add patches/support-gcc-5.x.patch. (Closes: #777767)
  - Fix build with gcc-5.x.
* Add patches/Disable-NET0001.als.patch.
  - Disable test of NET0001.als.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/mod/mth/shl/Newton.cpp
11
11
// - the copyright holder be liable for any  direct, indirect, incidental or -
12
12
// - special damages arising in any way out of the use of this software.     -
13
13
// ---------------------------------------------------------------------------
14
 
// - copyright (c) 1999-2012 amaury darsch                                   -
 
14
// - copyright (c) 1999-2015 amaury darsch                                   -
15
15
// ---------------------------------------------------------------------------
16
16
 
17
17
#include "Real.hpp"
19
19
#include "Vector.hpp"
20
20
#include "Newton.hpp"
21
21
#include "Linear.hpp"
 
22
#include "Boolean.hpp"
 
23
#include "Integer.hpp"
22
24
#include "Runnable.hpp"
23
25
#include "QuarkZone.hpp"
24
26
#include "Exception.hpp"
25
 
 
 
27
 
26
28
namespace afnix {
27
 
 
28
 
  // the default verification mode
29
 
  static const bool NTW_VF_DEF = false;
30
29
  // the iteration factor for the newton solver
31
 
  static const long NTW_NI_DEF = 100L;
 
30
  static const long NTW_MNI_DEF = 100L;
32
31
 
33
32
  // -------------------------------------------------------------------------
34
33
  // - class section                                                         -
35
34
  // -------------------------------------------------------------------------
36
35
 
37
 
  // create a default linear solver
 
36
  // create a default newton solver
38
37
 
39
38
  Newton::Newton (void) {
40
 
    d_avf = NTW_VF_DEF;
41
 
    d_mni = 0;
42
 
  }
43
 
  
44
 
  // create a solver by verification flag
45
 
 
46
 
  Newton::Newton (const bool avf) {
47
 
    d_avf = avf;
48
 
    d_mni = 0;
49
 
  }
50
 
 
51
 
  // create a solver by verification flag and mni value
52
 
 
53
 
  Newton::Newton (const bool avf, const long mni) {
54
 
    d_avf = avf;
55
 
    d_mni = mni;
56
 
  }
 
39
    d_mni = NTW_MNI_DEF;
 
40
    Object::iref (p_lnr = new Linear);
 
41
  }
 
42
 
 
43
  // create a newton solver with a linear solver
 
44
 
 
45
  Newton::Newton (Solver* lnr) {
 
46
    d_mni = NTW_MNI_DEF;
 
47
    Object::iref (p_lnr = lnr);
 
48
  } 
57
49
 
58
50
  // copy construct this object
59
51
 
60
52
  Newton::Newton (const Newton& that) {
61
53
    that.rdlock ();
62
54
    try {
63
 
      d_avf = that.d_avf;
64
55
      d_mni = that.d_mni;
 
56
      Object::iref (p_lnr = that.p_lnr);
65
57
      that.unlock ();
66
58
    } catch (...) {
67
59
      that.unlock ();
69
61
    }
70
62
  }
71
63
 
 
64
  // destroy this newton solver
 
65
 
 
66
  Newton::~Newton (void) {
 
67
    Object::dref (p_lnr);
 
68
  }
 
69
 
72
70
  // return the class name
73
71
  
74
72
  String Newton::repr (void) const {
90
88
    wrlock ();
91
89
    that.rdlock ();
92
90
    try {
93
 
      d_avf = that.d_avf;
94
91
      d_mni = that.d_mni;
 
92
      Object::iref (that.p_lnr);
 
93
      Object::dref (p_lnr);
 
94
      p_lnr = that.p_lnr;
95
95
      unlock ();
96
96
      that.unlock ();
97
97
      return *this;
102
102
    }
103
103
  }
104
104
 
 
105
  // set the newton linear solver
 
106
 
 
107
  void Newton::setlnr (Solver* lnr) {
 
108
    wrlock ();
 
109
    try {
 
110
      Object::iref (lnr);
 
111
      Object::dref (p_lnr);
 
112
      p_lnr = lnr;
 
113
      unlock ();
 
114
    } catch (...) {
 
115
      unlock ();
 
116
      throw;
 
117
    }
 
118
  }
 
119
 
 
120
  // get the newton linear solver
 
121
 
 
122
  Solver* Newton::getlnr (void) const {
 
123
    rdlock ();
 
124
    try {
 
125
      Solver* result = p_lnr;
 
126
      unlock ();
 
127
      return result;
 
128
    } catch (...) {
 
129
      unlock ();
 
130
      throw;
 
131
    }
 
132
  }
 
133
 
 
134
  // set the number of newton iterations
 
135
 
 
136
  void Newton::setmni (const long mni) {
 
137
    wrlock ();
 
138
    try {
 
139
      d_mni = mni;
 
140
      unlock ();
 
141
    } catch (...) {
 
142
      unlock ();
 
143
      throw;
 
144
    }
 
145
  }
 
146
 
 
147
  // get the number of newton iterations
 
148
 
 
149
  long Newton::getmni (void) const {
 
150
    rdlock ();
 
151
    try {
 
152
      long result = d_mni;
 
153
      unlock ();
 
154
      return result;
 
155
    } catch (...) {
 
156
      unlock ();
 
157
      throw;
 
158
    }
 
159
  }
 
160
 
105
161
  // find a function root with the newton method
106
162
 
107
163
  t_real Newton::solve (const Rfi& fo, const t_real xi) {
109
165
    try {
110
166
      // initialize current point
111
167
      t_real xn = xi;
112
 
      // get the number of iteration
113
 
      long ni = (d_mni <= 0) ? NTW_NI_DEF : d_mni;
114
168
      // newton loop
115
 
      for (long i = 0; i < ni; i++) {
 
169
      for (long k = 0L; k < d_mni; k++) {
116
170
        // check for zero
117
171
        t_real xz = fo.compute (xn);
118
172
        if (Math::rcmp (0, xz) == true) {
128
182
        // update delta
129
183
        xz /= xd;
130
184
        xn -= xz;
131
 
        // check for convergence
132
 
        if (Math::rcmp (xn, xz) == true) {
133
 
          unlock ();
134
 
          return xn;
135
 
        }
136
185
      }
137
186
      throw Exception ("newton-eror", "zero loop convergence failure");
138
187
    } catch (...) {
145
194
 
146
195
  Rvi* Newton::solve (const Rni& so, const Rvi& xi) {
147
196
    rdlock ();
148
 
    // get the number of iteration
149
 
    long ni = (d_mni <= 0) ? NTW_NI_DEF : d_mni;
150
197
    // initialize the current point
151
198
    Rvi* xn = dynamic_cast <Rvi*> (xi.clone ());
152
199
    Object::iref (xn);
155
202
    Rvi* rhs = nilp;
156
203
    Rmi* lhs = nilp;
157
204
    try {
158
 
      // create the linear solver
159
 
      Linear lnr (d_avf);
 
205
      // check for a linear solver
 
206
      if (p_lnr == nilp) {
 
207
        throw Exception ("newton-error", 
 
208
                         "no linear solver bound to newton solver");
 
209
      }
160
210
      // newton loop
161
 
      for (long i = 0; i < ni; i++) {
 
211
      for (long k = 0L; k < d_mni; k++) {
162
212
        // check rhs for zero
163
 
        rhs = so.getrhs (xn);
 
213
        rhs = so.getrhs (xn); Object::iref (rhs);
164
214
        if (rhs == nilp) {
165
215
          throw Exception ("newton-error", "null rhs in convergence loop");
166
216
        }
167
217
        if (Math::rcmp (0, rhs->norm ()) == true) {
168
218
          Object::tref (xn);
 
219
          Object::dref (rhs);
169
220
          unlock ();
170
221
          return xn;
171
222
        }
172
223
        // get the lhs by operand
173
 
        Rmi* lhs = so.getlhs (xn);
 
224
        Rmi* lhs = so.getlhs (xn); Object::iref (lhs);
174
225
        // try to solve the system
175
 
        xd = lnr.solve (*lhs, *rhs);
 
226
        xd = p_lnr->solve (lhs, *rhs);
176
227
        // clean elements
177
228
        Object::dref (lhs); lhs = nilp;
178
229
        Object::dref (rhs); rhs = nilp;
179
230
        // check for null delta
180
231
        if (xd == nilp) {
181
 
          throw Exception ("newton-error", "null delta in convergence loop");
 
232
          throw Exception ("newton-error", "newton linear solver failure");
182
233
        }
183
234
        // update delta
184
235
        *xn += *xd;
185
 
        // compute norm
186
 
        t_real nx = xn->norm ();
187
 
        t_real nd = xd->norm ();
188
236
        delete xd; xd = nilp;
189
 
        // check for null updated point
190
 
        if (Math::rcmp (0, nx) == true) {
191
 
          rhs = so.getrhs (xn);
192
 
          if (rhs == nilp) {
193
 
            throw Exception ("newton-error", "null rhs in convergence loop");
194
 
          }
195
 
          // check for valid point
196
 
          if (Math::rcmp (0, rhs->norm ()) == true) {
197
 
            Object::tref (xn);
198
 
            unlock ();
199
 
            return xn;
200
 
          }
201
 
          // wrong point
202
 
          throw Exception ("newton-error", "invalid point in convergence loop");
203
 
        }
204
 
        // check for convergence
205
 
        if (Math::rcmp (nd, nx) == true) {
206
 
          Object::tref (xn);
207
 
          unlock ();
208
 
          return xn;
209
 
        }
210
237
      }
211
 
      throw Exception ("newton-error", "zero loop convergence failure");
 
238
      throw Exception ("newton-error", "newton loop convergence failure");
212
239
    } catch (...) {
213
240
      delete xd;
214
241
      Object::dref (xn);
218
245
      throw;
219
246
    }
220
247
  }
221
 
  
 
248
 
222
249
  // -------------------------------------------------------------------------
223
250
  // - object section                                                        -
224
251
  // -------------------------------------------------------------------------
225
252
 
226
253
  // the quark zone
227
 
  static const long QUARK_ZONE_LENGTH = 1;
 
254
  static const long QUARK_ZONE_LENGTH = 5;
228
255
  static QuarkZone  zone (QUARK_ZONE_LENGTH);
229
256
 
230
 
  // the rvi supported quarks
231
 
  static const long QUARK_SOLVE = zone.intern ("solve");
 
257
  // the object supported quarks
 
258
  static const long QUARK_SOLVE  = zone.intern ("solve");
 
259
  static const long QUARK_SETLNR = zone.intern ("set-linear-solver");
 
260
  static const long QUARK_GETLNR = zone.intern ("get-linear-solver");
 
261
  static const long QUARK_SETMNI = zone.intern ("set-newton-iterations");
 
262
  static const long QUARK_GETMNI = zone.intern ("get-newton-iterations");
232
263
 
233
264
  // create a new object in a generic way
234
265
  
239
270
    if (argc == 0) return new Newton;
240
271
    // check for 1 argument
241
272
    if (argc == 1) {
242
 
      bool avf = argv->getbool (0);
243
 
      return new Newton (avf);
244
 
    }
245
 
    // check for 2 arguments
246
 
    if (argc == 2) {
247
 
      bool avf = argv->getbool (0);
248
 
      long mni = argv->getlong (1);
249
 
      return new Newton (avf, mni);
 
273
      Object* obj = argv->get (0);
 
274
      // check for a solver
 
275
      Solver* slv = dynamic_cast <Solver*> (obj);
 
276
      if (slv != nilp) return new Newton (slv);
 
277
      // invalid type
 
278
      throw Exception ("type-error", "invalid object for newton solver",
 
279
                       Object::repr (obj));
250
280
    }
251
281
    // invalid arguments
252
282
    throw Exception ("argument-error", "invalid arguments with newton object");
271
301
                         Vector* argv) {
272
302
    // get the number of arguments
273
303
    long argc = (argv == nilp) ? 0 : argv->length ();
274
 
    
 
304
 
 
305
    // dipatch 0 argument
 
306
    if (argc == 0) {
 
307
      if (quark == QUARK_GETMNI) return new Integer (getmni ());
 
308
      if (quark == QUARK_GETLNR) {
 
309
        rdlock ();
 
310
        try {
 
311
          Solver* result = getlnr ();
 
312
          robj->post (result);
 
313
          unlock ();
 
314
          return result;
 
315
        } catch (...) {
 
316
          unlock ();
 
317
          throw;
 
318
        }
 
319
      }
 
320
    }
 
321
    // dispatch 1 argument
 
322
    if (argc == 1) {
 
323
      if (quark == QUARK_SETMNI) {
 
324
        long mni = argv->getlong (0);
 
325
        setmni (mni);
 
326
        return nilp;
 
327
      }
 
328
      if (quark == QUARK_SETLNR) {
 
329
        Object* obj = argv->get (0);
 
330
        Solver* lnr = dynamic_cast <Solver*> (obj);
 
331
        if (lnr == nilp) {
 
332
          throw Exception ("type-error", 
 
333
                           "invalid object with set-linear-solver",
 
334
                           Object::repr (obj));
 
335
        }
 
336
        setlnr (lnr);
 
337
        return nilp;
 
338
      }
 
339
    }
275
340
    // dispatch 2 arguments
276
341
    if (argc == 2) {
277
342
      if (quark == QUARK_SOLVE) {