1 /*************************************************************************
2 Copyright (c) Sergey Bochkanov (ALGLIB project).
3
4 >>> SOURCE LICENSE >>>
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation (www.fsf.org); either version 2 of the
8 License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 A copy of the GNU General Public License is available at
16 http://www.fsf.org/licensing/licenses
17 >>> END OF LICENSE >>>
18 *************************************************************************/
19 #include "stdafx.h"
20 #include "solvers.h"
21
22 // disable some irrelevant warnings
23 #if (AE_COMPILER==AE_MSVC)
24 #pragma warning(disable:4100)
25 #pragma warning(disable:4127)
26 #pragma warning(disable:4702)
27 #pragma warning(disable:4996)
28 #endif
29 using namespace std;
30
31 /////////////////////////////////////////////////////////////////////////
32 //
33 // THIS SECTION CONTAINS IMPLEMENTATION OF C++ INTERFACE
34 //
35 /////////////////////////////////////////////////////////////////////////
36 namespace alglib
37 {
38
39
40 /*************************************************************************
41
42 *************************************************************************/
_densesolverreport_owner()43 _densesolverreport_owner::_densesolverreport_owner()
44 {
45 p_struct = (alglib_impl::densesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverreport), NULL);
46 if( p_struct==NULL )
47 throw ap_error("ALGLIB: malloc error");
48 if( !alglib_impl::_densesolverreport_init(p_struct, NULL, ae_false) )
49 throw ap_error("ALGLIB: malloc error");
50 }
51
_densesolverreport_owner(const _densesolverreport_owner & rhs)52 _densesolverreport_owner::_densesolverreport_owner(const _densesolverreport_owner &rhs)
53 {
54 p_struct = (alglib_impl::densesolverreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverreport), NULL);
55 if( p_struct==NULL )
56 throw ap_error("ALGLIB: malloc error");
57 if( !alglib_impl::_densesolverreport_init_copy(p_struct, const_cast<alglib_impl::densesolverreport*>(rhs.p_struct), NULL, ae_false) )
58 throw ap_error("ALGLIB: malloc error");
59 }
60
operator =(const _densesolverreport_owner & rhs)61 _densesolverreport_owner& _densesolverreport_owner::operator=(const _densesolverreport_owner &rhs)
62 {
63 if( this==&rhs )
64 return *this;
65 alglib_impl::_densesolverreport_clear(p_struct);
66 if( !alglib_impl::_densesolverreport_init_copy(p_struct, const_cast<alglib_impl::densesolverreport*>(rhs.p_struct), NULL, ae_false) )
67 throw ap_error("ALGLIB: malloc error");
68 return *this;
69 }
70
~_densesolverreport_owner()71 _densesolverreport_owner::~_densesolverreport_owner()
72 {
73 alglib_impl::_densesolverreport_clear(p_struct);
74 ae_free(p_struct);
75 }
76
c_ptr()77 alglib_impl::densesolverreport* _densesolverreport_owner::c_ptr()
78 {
79 return p_struct;
80 }
81
c_ptr() const82 alglib_impl::densesolverreport* _densesolverreport_owner::c_ptr() const
83 {
84 return const_cast<alglib_impl::densesolverreport*>(p_struct);
85 }
densesolverreport()86 densesolverreport::densesolverreport() : _densesolverreport_owner() ,r1(p_struct->r1),rinf(p_struct->rinf)
87 {
88 }
89
densesolverreport(const densesolverreport & rhs)90 densesolverreport::densesolverreport(const densesolverreport &rhs):_densesolverreport_owner(rhs) ,r1(p_struct->r1),rinf(p_struct->rinf)
91 {
92 }
93
operator =(const densesolverreport & rhs)94 densesolverreport& densesolverreport::operator=(const densesolverreport &rhs)
95 {
96 if( this==&rhs )
97 return *this;
98 _densesolverreport_owner::operator=(rhs);
99 return *this;
100 }
101
~densesolverreport()102 densesolverreport::~densesolverreport()
103 {
104 }
105
106
107 /*************************************************************************
108
109 *************************************************************************/
_densesolverlsreport_owner()110 _densesolverlsreport_owner::_densesolverlsreport_owner()
111 {
112 p_struct = (alglib_impl::densesolverlsreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverlsreport), NULL);
113 if( p_struct==NULL )
114 throw ap_error("ALGLIB: malloc error");
115 if( !alglib_impl::_densesolverlsreport_init(p_struct, NULL, ae_false) )
116 throw ap_error("ALGLIB: malloc error");
117 }
118
_densesolverlsreport_owner(const _densesolverlsreport_owner & rhs)119 _densesolverlsreport_owner::_densesolverlsreport_owner(const _densesolverlsreport_owner &rhs)
120 {
121 p_struct = (alglib_impl::densesolverlsreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::densesolverlsreport), NULL);
122 if( p_struct==NULL )
123 throw ap_error("ALGLIB: malloc error");
124 if( !alglib_impl::_densesolverlsreport_init_copy(p_struct, const_cast<alglib_impl::densesolverlsreport*>(rhs.p_struct), NULL, ae_false) )
125 throw ap_error("ALGLIB: malloc error");
126 }
127
operator =(const _densesolverlsreport_owner & rhs)128 _densesolverlsreport_owner& _densesolverlsreport_owner::operator=(const _densesolverlsreport_owner &rhs)
129 {
130 if( this==&rhs )
131 return *this;
132 alglib_impl::_densesolverlsreport_clear(p_struct);
133 if( !alglib_impl::_densesolverlsreport_init_copy(p_struct, const_cast<alglib_impl::densesolverlsreport*>(rhs.p_struct), NULL, ae_false) )
134 throw ap_error("ALGLIB: malloc error");
135 return *this;
136 }
137
~_densesolverlsreport_owner()138 _densesolverlsreport_owner::~_densesolverlsreport_owner()
139 {
140 alglib_impl::_densesolverlsreport_clear(p_struct);
141 ae_free(p_struct);
142 }
143
c_ptr()144 alglib_impl::densesolverlsreport* _densesolverlsreport_owner::c_ptr()
145 {
146 return p_struct;
147 }
148
c_ptr() const149 alglib_impl::densesolverlsreport* _densesolverlsreport_owner::c_ptr() const
150 {
151 return const_cast<alglib_impl::densesolverlsreport*>(p_struct);
152 }
densesolverlsreport()153 densesolverlsreport::densesolverlsreport() : _densesolverlsreport_owner() ,r2(p_struct->r2),cx(&p_struct->cx),n(p_struct->n),k(p_struct->k)
154 {
155 }
156
densesolverlsreport(const densesolverlsreport & rhs)157 densesolverlsreport::densesolverlsreport(const densesolverlsreport &rhs):_densesolverlsreport_owner(rhs) ,r2(p_struct->r2),cx(&p_struct->cx),n(p_struct->n),k(p_struct->k)
158 {
159 }
160
operator =(const densesolverlsreport & rhs)161 densesolverlsreport& densesolverlsreport::operator=(const densesolverlsreport &rhs)
162 {
163 if( this==&rhs )
164 return *this;
165 _densesolverlsreport_owner::operator=(rhs);
166 return *this;
167 }
168
~densesolverlsreport()169 densesolverlsreport::~densesolverlsreport()
170 {
171 }
172
173 /*************************************************************************
174 Dense solver.
175
176 This subroutine solves a system A*x=b, where A is NxN non-denegerate
177 real matrix, x and b are vectors.
178
179 Algorithm features:
180 * automatic detection of degenerate cases
181 * condition number estimation
182 * iterative refinement
183 * O(N^3) complexity
184
185 INPUT PARAMETERS
186 A - array[0..N-1,0..N-1], system matrix
187 N - size of A
188 B - array[0..N-1], right part
189
190 OUTPUT PARAMETERS
191 Info - return code:
192 * -3 A is singular, or VERY close to singular.
193 X is filled by zeros in such cases.
194 * -1 N<=0 was passed
195 * 1 task is solved (but matrix A may be ill-conditioned,
196 check R1/RInf parameters for condition numbers).
197 Rep - solver report, see below for more info
198 X - array[0..N-1], it contains:
199 * solution of A*x=b if A is non-singular (well-conditioned
200 or ill-conditioned, but not very close to singular)
201 * zeros, if A is singular or VERY close to singular
202 (in this case Info=-3).
203
204 SOLVER REPORT
205
206 Subroutine sets following fields of the Rep structure:
207 * R1 reciprocal of condition number: 1/cond(A), 1-norm.
208 * RInf reciprocal of condition number: 1/cond(A), inf-norm.
209
210 -- ALGLIB --
211 Copyright 27.01.2010 by Bochkanov Sergey
212 *************************************************************************/
rmatrixsolve(const real_2d_array & a,const ae_int_t n,const real_1d_array & b,ae_int_t & info,densesolverreport & rep,real_1d_array & x)213 void rmatrixsolve(const real_2d_array &a, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x)
214 {
215 alglib_impl::ae_state _alglib_env_state;
216 alglib_impl::ae_state_init(&_alglib_env_state);
217 try
218 {
219 alglib_impl::rmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
220 alglib_impl::ae_state_clear(&_alglib_env_state);
221 return;
222 }
223 catch(alglib_impl::ae_error_type)
224 {
225 throw ap_error(_alglib_env_state.error_msg);
226 }
227 catch(...)
228 {
229 throw;
230 }
231 }
232
233 /*************************************************************************
234 Dense solver.
235
236 Similar to RMatrixSolve() but solves task with multiple right parts (where
237 b and x are NxM matrices).
238
239 Algorithm features:
240 * automatic detection of degenerate cases
241 * condition number estimation
242 * optional iterative refinement
243 * O(N^3+M*N^2) complexity
244
245 INPUT PARAMETERS
246 A - array[0..N-1,0..N-1], system matrix
247 N - size of A
248 B - array[0..N-1,0..M-1], right part
249 M - right part size
250 RFS - iterative refinement switch:
251 * True - refinement is used.
252 Less performance, more precision.
253 * False - refinement is not used.
254 More performance, less precision.
255
256 OUTPUT PARAMETERS
257 Info - same as in RMatrixSolve
258 Rep - same as in RMatrixSolve
259 X - same as in RMatrixSolve
260
261 -- ALGLIB --
262 Copyright 27.01.2010 by Bochkanov Sergey
263 *************************************************************************/
rmatrixsolvem(const real_2d_array & a,const ae_int_t n,const real_2d_array & b,const ae_int_t m,const bool rfs,ae_int_t & info,densesolverreport & rep,real_2d_array & x)264 void rmatrixsolvem(const real_2d_array &a, const ae_int_t n, const real_2d_array &b, const ae_int_t m, const bool rfs, ae_int_t &info, densesolverreport &rep, real_2d_array &x)
265 {
266 alglib_impl::ae_state _alglib_env_state;
267 alglib_impl::ae_state_init(&_alglib_env_state);
268 try
269 {
270 alglib_impl::rmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, rfs, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
271 alglib_impl::ae_state_clear(&_alglib_env_state);
272 return;
273 }
274 catch(alglib_impl::ae_error_type)
275 {
276 throw ap_error(_alglib_env_state.error_msg);
277 }
278 catch(...)
279 {
280 throw;
281 }
282 }
283
284 /*************************************************************************
285 Dense solver.
286
287 This subroutine solves a system A*X=B, where A is NxN non-denegerate
288 real matrix given by its LU decomposition, X and B are NxM real matrices.
289
290 Algorithm features:
291 * automatic detection of degenerate cases
292 * O(N^2) complexity
293 * condition number estimation
294
295 No iterative refinement is provided because exact form of original matrix
296 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
297 need iterative refinement.
298
299 INPUT PARAMETERS
300 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
301 P - array[0..N-1], pivots array, RMatrixLU result
302 N - size of A
303 B - array[0..N-1], right part
304
305 OUTPUT PARAMETERS
306 Info - same as in RMatrixSolve
307 Rep - same as in RMatrixSolve
308 X - same as in RMatrixSolve
309
310 -- ALGLIB --
311 Copyright 27.01.2010 by Bochkanov Sergey
312 *************************************************************************/
rmatrixlusolve(const real_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const real_1d_array & b,ae_int_t & info,densesolverreport & rep,real_1d_array & x)313 void rmatrixlusolve(const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x)
314 {
315 alglib_impl::ae_state _alglib_env_state;
316 alglib_impl::ae_state_init(&_alglib_env_state);
317 try
318 {
319 alglib_impl::rmatrixlusolve(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
320 alglib_impl::ae_state_clear(&_alglib_env_state);
321 return;
322 }
323 catch(alglib_impl::ae_error_type)
324 {
325 throw ap_error(_alglib_env_state.error_msg);
326 }
327 catch(...)
328 {
329 throw;
330 }
331 }
332
333 /*************************************************************************
334 Dense solver.
335
336 Similar to RMatrixLUSolve() but solves task with multiple right parts
337 (where b and x are NxM matrices).
338
339 Algorithm features:
340 * automatic detection of degenerate cases
341 * O(M*N^2) complexity
342 * condition number estimation
343
344 No iterative refinement is provided because exact form of original matrix
345 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
346 need iterative refinement.
347
348 INPUT PARAMETERS
349 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
350 P - array[0..N-1], pivots array, RMatrixLU result
351 N - size of A
352 B - array[0..N-1,0..M-1], right part
353 M - right part size
354
355 OUTPUT PARAMETERS
356 Info - same as in RMatrixSolve
357 Rep - same as in RMatrixSolve
358 X - same as in RMatrixSolve
359
360 -- ALGLIB --
361 Copyright 27.01.2010 by Bochkanov Sergey
362 *************************************************************************/
rmatrixlusolvem(const real_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const real_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,real_2d_array & x)363 void rmatrixlusolvem(const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x)
364 {
365 alglib_impl::ae_state _alglib_env_state;
366 alglib_impl::ae_state_init(&_alglib_env_state);
367 try
368 {
369 alglib_impl::rmatrixlusolvem(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
370 alglib_impl::ae_state_clear(&_alglib_env_state);
371 return;
372 }
373 catch(alglib_impl::ae_error_type)
374 {
375 throw ap_error(_alglib_env_state.error_msg);
376 }
377 catch(...)
378 {
379 throw;
380 }
381 }
382
383 /*************************************************************************
384 Dense solver.
385
386 This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS
387 LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have
388 both A and its LU decomposition.
389
390 Algorithm features:
391 * automatic detection of degenerate cases
392 * condition number estimation
393 * iterative refinement
394 * O(N^2) complexity
395
396 INPUT PARAMETERS
397 A - array[0..N-1,0..N-1], system matrix
398 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
399 P - array[0..N-1], pivots array, RMatrixLU result
400 N - size of A
401 B - array[0..N-1], right part
402
403 OUTPUT PARAMETERS
404 Info - same as in RMatrixSolveM
405 Rep - same as in RMatrixSolveM
406 X - same as in RMatrixSolveM
407
408 -- ALGLIB --
409 Copyright 27.01.2010 by Bochkanov Sergey
410 *************************************************************************/
rmatrixmixedsolve(const real_2d_array & a,const real_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const real_1d_array & b,ae_int_t & info,densesolverreport & rep,real_1d_array & x)411 void rmatrixmixedsolve(const real_2d_array &a, const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x)
412 {
413 alglib_impl::ae_state _alglib_env_state;
414 alglib_impl::ae_state_init(&_alglib_env_state);
415 try
416 {
417 alglib_impl::rmatrixmixedsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
418 alglib_impl::ae_state_clear(&_alglib_env_state);
419 return;
420 }
421 catch(alglib_impl::ae_error_type)
422 {
423 throw ap_error(_alglib_env_state.error_msg);
424 }
425 catch(...)
426 {
427 throw;
428 }
429 }
430
431 /*************************************************************************
432 Dense solver.
433
434 Similar to RMatrixMixedSolve() but solves task with multiple right parts
435 (where b and x are NxM matrices).
436
437 Algorithm features:
438 * automatic detection of degenerate cases
439 * condition number estimation
440 * iterative refinement
441 * O(M*N^2) complexity
442
443 INPUT PARAMETERS
444 A - array[0..N-1,0..N-1], system matrix
445 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
446 P - array[0..N-1], pivots array, RMatrixLU result
447 N - size of A
448 B - array[0..N-1,0..M-1], right part
449 M - right part size
450
451 OUTPUT PARAMETERS
452 Info - same as in RMatrixSolveM
453 Rep - same as in RMatrixSolveM
454 X - same as in RMatrixSolveM
455
456 -- ALGLIB --
457 Copyright 27.01.2010 by Bochkanov Sergey
458 *************************************************************************/
rmatrixmixedsolvem(const real_2d_array & a,const real_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const real_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,real_2d_array & x)459 void rmatrixmixedsolvem(const real_2d_array &a, const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x)
460 {
461 alglib_impl::ae_state _alglib_env_state;
462 alglib_impl::ae_state_init(&_alglib_env_state);
463 try
464 {
465 alglib_impl::rmatrixmixedsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
466 alglib_impl::ae_state_clear(&_alglib_env_state);
467 return;
468 }
469 catch(alglib_impl::ae_error_type)
470 {
471 throw ap_error(_alglib_env_state.error_msg);
472 }
473 catch(...)
474 {
475 throw;
476 }
477 }
478
479 /*************************************************************************
480 Dense solver. Same as RMatrixSolveM(), but for complex matrices.
481
482 Algorithm features:
483 * automatic detection of degenerate cases
484 * condition number estimation
485 * iterative refinement
486 * O(N^3+M*N^2) complexity
487
488 INPUT PARAMETERS
489 A - array[0..N-1,0..N-1], system matrix
490 N - size of A
491 B - array[0..N-1,0..M-1], right part
492 M - right part size
493 RFS - iterative refinement switch:
494 * True - refinement is used.
495 Less performance, more precision.
496 * False - refinement is not used.
497 More performance, less precision.
498
499 OUTPUT PARAMETERS
500 Info - same as in RMatrixSolve
501 Rep - same as in RMatrixSolve
502 X - same as in RMatrixSolve
503
504 -- ALGLIB --
505 Copyright 27.01.2010 by Bochkanov Sergey
506 *************************************************************************/
cmatrixsolvem(const complex_2d_array & a,const ae_int_t n,const complex_2d_array & b,const ae_int_t m,const bool rfs,ae_int_t & info,densesolverreport & rep,complex_2d_array & x)507 void cmatrixsolvem(const complex_2d_array &a, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, const bool rfs, ae_int_t &info, densesolverreport &rep, complex_2d_array &x)
508 {
509 alglib_impl::ae_state _alglib_env_state;
510 alglib_impl::ae_state_init(&_alglib_env_state);
511 try
512 {
513 alglib_impl::cmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, rfs, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
514 alglib_impl::ae_state_clear(&_alglib_env_state);
515 return;
516 }
517 catch(alglib_impl::ae_error_type)
518 {
519 throw ap_error(_alglib_env_state.error_msg);
520 }
521 catch(...)
522 {
523 throw;
524 }
525 }
526
527 /*************************************************************************
528 Dense solver. Same as RMatrixSolve(), but for complex matrices.
529
530 Algorithm features:
531 * automatic detection of degenerate cases
532 * condition number estimation
533 * iterative refinement
534 * O(N^3) complexity
535
536 INPUT PARAMETERS
537 A - array[0..N-1,0..N-1], system matrix
538 N - size of A
539 B - array[0..N-1], right part
540
541 OUTPUT PARAMETERS
542 Info - same as in RMatrixSolve
543 Rep - same as in RMatrixSolve
544 X - same as in RMatrixSolve
545
546 -- ALGLIB --
547 Copyright 27.01.2010 by Bochkanov Sergey
548 *************************************************************************/
cmatrixsolve(const complex_2d_array & a,const ae_int_t n,const complex_1d_array & b,ae_int_t & info,densesolverreport & rep,complex_1d_array & x)549 void cmatrixsolve(const complex_2d_array &a, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x)
550 {
551 alglib_impl::ae_state _alglib_env_state;
552 alglib_impl::ae_state_init(&_alglib_env_state);
553 try
554 {
555 alglib_impl::cmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
556 alglib_impl::ae_state_clear(&_alglib_env_state);
557 return;
558 }
559 catch(alglib_impl::ae_error_type)
560 {
561 throw ap_error(_alglib_env_state.error_msg);
562 }
563 catch(...)
564 {
565 throw;
566 }
567 }
568
569 /*************************************************************************
570 Dense solver. Same as RMatrixLUSolveM(), but for complex matrices.
571
572 Algorithm features:
573 * automatic detection of degenerate cases
574 * O(M*N^2) complexity
575 * condition number estimation
576
577 No iterative refinement is provided because exact form of original matrix
578 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
579 need iterative refinement.
580
581 INPUT PARAMETERS
582 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
583 P - array[0..N-1], pivots array, RMatrixLU result
584 N - size of A
585 B - array[0..N-1,0..M-1], right part
586 M - right part size
587
588 OUTPUT PARAMETERS
589 Info - same as in RMatrixSolve
590 Rep - same as in RMatrixSolve
591 X - same as in RMatrixSolve
592
593 -- ALGLIB --
594 Copyright 27.01.2010 by Bochkanov Sergey
595 *************************************************************************/
cmatrixlusolvem(const complex_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const complex_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,complex_2d_array & x)596 void cmatrixlusolvem(const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x)
597 {
598 alglib_impl::ae_state _alglib_env_state;
599 alglib_impl::ae_state_init(&_alglib_env_state);
600 try
601 {
602 alglib_impl::cmatrixlusolvem(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
603 alglib_impl::ae_state_clear(&_alglib_env_state);
604 return;
605 }
606 catch(alglib_impl::ae_error_type)
607 {
608 throw ap_error(_alglib_env_state.error_msg);
609 }
610 catch(...)
611 {
612 throw;
613 }
614 }
615
616 /*************************************************************************
617 Dense solver. Same as RMatrixLUSolve(), but for complex matrices.
618
619 Algorithm features:
620 * automatic detection of degenerate cases
621 * O(N^2) complexity
622 * condition number estimation
623
624 No iterative refinement is provided because exact form of original matrix
625 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
626 need iterative refinement.
627
628 INPUT PARAMETERS
629 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
630 P - array[0..N-1], pivots array, CMatrixLU result
631 N - size of A
632 B - array[0..N-1], right part
633
634 OUTPUT PARAMETERS
635 Info - same as in RMatrixSolve
636 Rep - same as in RMatrixSolve
637 X - same as in RMatrixSolve
638
639 -- ALGLIB --
640 Copyright 27.01.2010 by Bochkanov Sergey
641 *************************************************************************/
cmatrixlusolve(const complex_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const complex_1d_array & b,ae_int_t & info,densesolverreport & rep,complex_1d_array & x)642 void cmatrixlusolve(const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x)
643 {
644 alglib_impl::ae_state _alglib_env_state;
645 alglib_impl::ae_state_init(&_alglib_env_state);
646 try
647 {
648 alglib_impl::cmatrixlusolve(const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
649 alglib_impl::ae_state_clear(&_alglib_env_state);
650 return;
651 }
652 catch(alglib_impl::ae_error_type)
653 {
654 throw ap_error(_alglib_env_state.error_msg);
655 }
656 catch(...)
657 {
658 throw;
659 }
660 }
661
662 /*************************************************************************
663 Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices.
664
665 Algorithm features:
666 * automatic detection of degenerate cases
667 * condition number estimation
668 * iterative refinement
669 * O(M*N^2) complexity
670
671 INPUT PARAMETERS
672 A - array[0..N-1,0..N-1], system matrix
673 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
674 P - array[0..N-1], pivots array, CMatrixLU result
675 N - size of A
676 B - array[0..N-1,0..M-1], right part
677 M - right part size
678
679 OUTPUT PARAMETERS
680 Info - same as in RMatrixSolveM
681 Rep - same as in RMatrixSolveM
682 X - same as in RMatrixSolveM
683
684 -- ALGLIB --
685 Copyright 27.01.2010 by Bochkanov Sergey
686 *************************************************************************/
cmatrixmixedsolvem(const complex_2d_array & a,const complex_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const complex_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,complex_2d_array & x)687 void cmatrixmixedsolvem(const complex_2d_array &a, const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x)
688 {
689 alglib_impl::ae_state _alglib_env_state;
690 alglib_impl::ae_state_init(&_alglib_env_state);
691 try
692 {
693 alglib_impl::cmatrixmixedsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
694 alglib_impl::ae_state_clear(&_alglib_env_state);
695 return;
696 }
697 catch(alglib_impl::ae_error_type)
698 {
699 throw ap_error(_alglib_env_state.error_msg);
700 }
701 catch(...)
702 {
703 throw;
704 }
705 }
706
707 /*************************************************************************
708 Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.
709
710 Algorithm features:
711 * automatic detection of degenerate cases
712 * condition number estimation
713 * iterative refinement
714 * O(N^2) complexity
715
716 INPUT PARAMETERS
717 A - array[0..N-1,0..N-1], system matrix
718 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
719 P - array[0..N-1], pivots array, CMatrixLU result
720 N - size of A
721 B - array[0..N-1], right part
722
723 OUTPUT PARAMETERS
724 Info - same as in RMatrixSolveM
725 Rep - same as in RMatrixSolveM
726 X - same as in RMatrixSolveM
727
728 -- ALGLIB --
729 Copyright 27.01.2010 by Bochkanov Sergey
730 *************************************************************************/
cmatrixmixedsolve(const complex_2d_array & a,const complex_2d_array & lua,const integer_1d_array & p,const ae_int_t n,const complex_1d_array & b,ae_int_t & info,densesolverreport & rep,complex_1d_array & x)731 void cmatrixmixedsolve(const complex_2d_array &a, const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x)
732 {
733 alglib_impl::ae_state _alglib_env_state;
734 alglib_impl::ae_state_init(&_alglib_env_state);
735 try
736 {
737 alglib_impl::cmatrixmixedsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), const_cast<alglib_impl::ae_matrix*>(lua.c_ptr()), const_cast<alglib_impl::ae_vector*>(p.c_ptr()), n, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
738 alglib_impl::ae_state_clear(&_alglib_env_state);
739 return;
740 }
741 catch(alglib_impl::ae_error_type)
742 {
743 throw ap_error(_alglib_env_state.error_msg);
744 }
745 catch(...)
746 {
747 throw;
748 }
749 }
750
751 /*************************************************************************
752 Dense solver. Same as RMatrixSolveM(), but for symmetric positive definite
753 matrices.
754
755 Algorithm features:
756 * automatic detection of degenerate cases
757 * condition number estimation
758 * O(N^3+M*N^2) complexity
759 * matrix is represented by its upper or lower triangle
760
761 No iterative refinement is provided because such partial representation of
762 matrix does not allow efficient calculation of extra-precise matrix-vector
763 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
764 need iterative refinement.
765
766 INPUT PARAMETERS
767 A - array[0..N-1,0..N-1], system matrix
768 N - size of A
769 IsUpper - what half of A is provided
770 B - array[0..N-1,0..M-1], right part
771 M - right part size
772
773 OUTPUT PARAMETERS
774 Info - same as in RMatrixSolve.
775 Returns -3 for non-SPD matrices.
776 Rep - same as in RMatrixSolve
777 X - same as in RMatrixSolve
778
779 -- ALGLIB --
780 Copyright 27.01.2010 by Bochkanov Sergey
781 *************************************************************************/
spdmatrixsolvem(const real_2d_array & a,const ae_int_t n,const bool isupper,const real_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,real_2d_array & x)782 void spdmatrixsolvem(const real_2d_array &a, const ae_int_t n, const bool isupper, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x)
783 {
784 alglib_impl::ae_state _alglib_env_state;
785 alglib_impl::ae_state_init(&_alglib_env_state);
786 try
787 {
788 alglib_impl::spdmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
789 alglib_impl::ae_state_clear(&_alglib_env_state);
790 return;
791 }
792 catch(alglib_impl::ae_error_type)
793 {
794 throw ap_error(_alglib_env_state.error_msg);
795 }
796 catch(...)
797 {
798 throw;
799 }
800 }
801
802 /*************************************************************************
803 Dense solver. Same as RMatrixSolve(), but for SPD matrices.
804
805 Algorithm features:
806 * automatic detection of degenerate cases
807 * condition number estimation
808 * O(N^3) complexity
809 * matrix is represented by its upper or lower triangle
810
811 No iterative refinement is provided because such partial representation of
812 matrix does not allow efficient calculation of extra-precise matrix-vector
813 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
814 need iterative refinement.
815
816 INPUT PARAMETERS
817 A - array[0..N-1,0..N-1], system matrix
818 N - size of A
819 IsUpper - what half of A is provided
820 B - array[0..N-1], right part
821
822 OUTPUT PARAMETERS
823 Info - same as in RMatrixSolve
824 Returns -3 for non-SPD matrices.
825 Rep - same as in RMatrixSolve
826 X - same as in RMatrixSolve
827
828 -- ALGLIB --
829 Copyright 27.01.2010 by Bochkanov Sergey
830 *************************************************************************/
spdmatrixsolve(const real_2d_array & a,const ae_int_t n,const bool isupper,const real_1d_array & b,ae_int_t & info,densesolverreport & rep,real_1d_array & x)831 void spdmatrixsolve(const real_2d_array &a, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x)
832 {
833 alglib_impl::ae_state _alglib_env_state;
834 alglib_impl::ae_state_init(&_alglib_env_state);
835 try
836 {
837 alglib_impl::spdmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
838 alglib_impl::ae_state_clear(&_alglib_env_state);
839 return;
840 }
841 catch(alglib_impl::ae_error_type)
842 {
843 throw ap_error(_alglib_env_state.error_msg);
844 }
845 catch(...)
846 {
847 throw;
848 }
849 }
850
851 /*************************************************************************
852 Dense solver. Same as RMatrixLUSolveM(), but for SPD matrices represented
853 by their Cholesky decomposition.
854
855 Algorithm features:
856 * automatic detection of degenerate cases
857 * O(M*N^2) complexity
858 * condition number estimation
859 * matrix is represented by its upper or lower triangle
860
861 No iterative refinement is provided because such partial representation of
862 matrix does not allow efficient calculation of extra-precise matrix-vector
863 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
864 need iterative refinement.
865
866 INPUT PARAMETERS
867 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
868 SPDMatrixCholesky result
869 N - size of CHA
870 IsUpper - what half of CHA is provided
871 B - array[0..N-1,0..M-1], right part
872 M - right part size
873
874 OUTPUT PARAMETERS
875 Info - same as in RMatrixSolve
876 Rep - same as in RMatrixSolve
877 X - same as in RMatrixSolve
878
879 -- ALGLIB --
880 Copyright 27.01.2010 by Bochkanov Sergey
881 *************************************************************************/
spdmatrixcholeskysolvem(const real_2d_array & cha,const ae_int_t n,const bool isupper,const real_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,real_2d_array & x)882 void spdmatrixcholeskysolvem(const real_2d_array &cha, const ae_int_t n, const bool isupper, const real_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, real_2d_array &x)
883 {
884 alglib_impl::ae_state _alglib_env_state;
885 alglib_impl::ae_state_init(&_alglib_env_state);
886 try
887 {
888 alglib_impl::spdmatrixcholeskysolvem(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
889 alglib_impl::ae_state_clear(&_alglib_env_state);
890 return;
891 }
892 catch(alglib_impl::ae_error_type)
893 {
894 throw ap_error(_alglib_env_state.error_msg);
895 }
896 catch(...)
897 {
898 throw;
899 }
900 }
901
902 /*************************************************************************
903 Dense solver. Same as RMatrixLUSolve(), but for SPD matrices represented
904 by their Cholesky decomposition.
905
906 Algorithm features:
907 * automatic detection of degenerate cases
908 * O(N^2) complexity
909 * condition number estimation
910 * matrix is represented by its upper or lower triangle
911
912 No iterative refinement is provided because such partial representation of
913 matrix does not allow efficient calculation of extra-precise matrix-vector
914 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
915 need iterative refinement.
916
917 INPUT PARAMETERS
918 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
919 SPDMatrixCholesky result
920 N - size of A
921 IsUpper - what half of CHA is provided
922 B - array[0..N-1], right part
923
924 OUTPUT PARAMETERS
925 Info - same as in RMatrixSolve
926 Rep - same as in RMatrixSolve
927 X - same as in RMatrixSolve
928
929 -- ALGLIB --
930 Copyright 27.01.2010 by Bochkanov Sergey
931 *************************************************************************/
spdmatrixcholeskysolve(const real_2d_array & cha,const ae_int_t n,const bool isupper,const real_1d_array & b,ae_int_t & info,densesolverreport & rep,real_1d_array & x)932 void spdmatrixcholeskysolve(const real_2d_array &cha, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, densesolverreport &rep, real_1d_array &x)
933 {
934 alglib_impl::ae_state _alglib_env_state;
935 alglib_impl::ae_state_init(&_alglib_env_state);
936 try
937 {
938 alglib_impl::spdmatrixcholeskysolve(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
939 alglib_impl::ae_state_clear(&_alglib_env_state);
940 return;
941 }
942 catch(alglib_impl::ae_error_type)
943 {
944 throw ap_error(_alglib_env_state.error_msg);
945 }
946 catch(...)
947 {
948 throw;
949 }
950 }
951
952 /*************************************************************************
953 Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
954 matrices.
955
956 Algorithm features:
957 * automatic detection of degenerate cases
958 * condition number estimation
959 * O(N^3+M*N^2) complexity
960 * matrix is represented by its upper or lower triangle
961
962 No iterative refinement is provided because such partial representation of
963 matrix does not allow efficient calculation of extra-precise matrix-vector
964 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
965 need iterative refinement.
966
967 INPUT PARAMETERS
968 A - array[0..N-1,0..N-1], system matrix
969 N - size of A
970 IsUpper - what half of A is provided
971 B - array[0..N-1,0..M-1], right part
972 M - right part size
973
974 OUTPUT PARAMETERS
975 Info - same as in RMatrixSolve.
976 Returns -3 for non-HPD matrices.
977 Rep - same as in RMatrixSolve
978 X - same as in RMatrixSolve
979
980 -- ALGLIB --
981 Copyright 27.01.2010 by Bochkanov Sergey
982 *************************************************************************/
hpdmatrixsolvem(const complex_2d_array & a,const ae_int_t n,const bool isupper,const complex_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,complex_2d_array & x)983 void hpdmatrixsolvem(const complex_2d_array &a, const ae_int_t n, const bool isupper, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x)
984 {
985 alglib_impl::ae_state _alglib_env_state;
986 alglib_impl::ae_state_init(&_alglib_env_state);
987 try
988 {
989 alglib_impl::hpdmatrixsolvem(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
990 alglib_impl::ae_state_clear(&_alglib_env_state);
991 return;
992 }
993 catch(alglib_impl::ae_error_type)
994 {
995 throw ap_error(_alglib_env_state.error_msg);
996 }
997 catch(...)
998 {
999 throw;
1000 }
1001 }
1002
1003 /*************************************************************************
1004 Dense solver. Same as RMatrixSolve(), but for Hermitian positive definite
1005 matrices.
1006
1007 Algorithm features:
1008 * automatic detection of degenerate cases
1009 * condition number estimation
1010 * O(N^3) complexity
1011 * matrix is represented by its upper or lower triangle
1012
1013 No iterative refinement is provided because such partial representation of
1014 matrix does not allow efficient calculation of extra-precise matrix-vector
1015 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1016 need iterative refinement.
1017
1018 INPUT PARAMETERS
1019 A - array[0..N-1,0..N-1], system matrix
1020 N - size of A
1021 IsUpper - what half of A is provided
1022 B - array[0..N-1], right part
1023
1024 OUTPUT PARAMETERS
1025 Info - same as in RMatrixSolve
1026 Returns -3 for non-HPD matrices.
1027 Rep - same as in RMatrixSolve
1028 X - same as in RMatrixSolve
1029
1030 -- ALGLIB --
1031 Copyright 27.01.2010 by Bochkanov Sergey
1032 *************************************************************************/
hpdmatrixsolve(const complex_2d_array & a,const ae_int_t n,const bool isupper,const complex_1d_array & b,ae_int_t & info,densesolverreport & rep,complex_1d_array & x)1033 void hpdmatrixsolve(const complex_2d_array &a, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x)
1034 {
1035 alglib_impl::ae_state _alglib_env_state;
1036 alglib_impl::ae_state_init(&_alglib_env_state);
1037 try
1038 {
1039 alglib_impl::hpdmatrixsolve(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
1040 alglib_impl::ae_state_clear(&_alglib_env_state);
1041 return;
1042 }
1043 catch(alglib_impl::ae_error_type)
1044 {
1045 throw ap_error(_alglib_env_state.error_msg);
1046 }
1047 catch(...)
1048 {
1049 throw;
1050 }
1051 }
1052
1053 /*************************************************************************
1054 Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices represented
1055 by their Cholesky decomposition.
1056
1057 Algorithm features:
1058 * automatic detection of degenerate cases
1059 * O(M*N^2) complexity
1060 * condition number estimation
1061 * matrix is represented by its upper or lower triangle
1062
1063 No iterative refinement is provided because such partial representation of
1064 matrix does not allow efficient calculation of extra-precise matrix-vector
1065 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1066 need iterative refinement.
1067
1068 INPUT PARAMETERS
1069 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1070 HPDMatrixCholesky result
1071 N - size of CHA
1072 IsUpper - what half of CHA is provided
1073 B - array[0..N-1,0..M-1], right part
1074 M - right part size
1075
1076 OUTPUT PARAMETERS
1077 Info - same as in RMatrixSolve
1078 Rep - same as in RMatrixSolve
1079 X - same as in RMatrixSolve
1080
1081 -- ALGLIB --
1082 Copyright 27.01.2010 by Bochkanov Sergey
1083 *************************************************************************/
hpdmatrixcholeskysolvem(const complex_2d_array & cha,const ae_int_t n,const bool isupper,const complex_2d_array & b,const ae_int_t m,ae_int_t & info,densesolverreport & rep,complex_2d_array & x)1084 void hpdmatrixcholeskysolvem(const complex_2d_array &cha, const ae_int_t n, const bool isupper, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, densesolverreport &rep, complex_2d_array &x)
1085 {
1086 alglib_impl::ae_state _alglib_env_state;
1087 alglib_impl::ae_state_init(&_alglib_env_state);
1088 try
1089 {
1090 alglib_impl::hpdmatrixcholeskysolvem(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_matrix*>(b.c_ptr()), m, &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_matrix*>(x.c_ptr()), &_alglib_env_state);
1091 alglib_impl::ae_state_clear(&_alglib_env_state);
1092 return;
1093 }
1094 catch(alglib_impl::ae_error_type)
1095 {
1096 throw ap_error(_alglib_env_state.error_msg);
1097 }
1098 catch(...)
1099 {
1100 throw;
1101 }
1102 }
1103
1104 /*************************************************************************
1105 Dense solver. Same as RMatrixLUSolve(), but for HPD matrices represented
1106 by their Cholesky decomposition.
1107
1108 Algorithm features:
1109 * automatic detection of degenerate cases
1110 * O(N^2) complexity
1111 * condition number estimation
1112 * matrix is represented by its upper or lower triangle
1113
1114 No iterative refinement is provided because such partial representation of
1115 matrix does not allow efficient calculation of extra-precise matrix-vector
1116 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1117 need iterative refinement.
1118
1119 INPUT PARAMETERS
1120 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1121 SPDMatrixCholesky result
1122 N - size of A
1123 IsUpper - what half of CHA is provided
1124 B - array[0..N-1], right part
1125
1126 OUTPUT PARAMETERS
1127 Info - same as in RMatrixSolve
1128 Rep - same as in RMatrixSolve
1129 X - same as in RMatrixSolve
1130
1131 -- ALGLIB --
1132 Copyright 27.01.2010 by Bochkanov Sergey
1133 *************************************************************************/
hpdmatrixcholeskysolve(const complex_2d_array & cha,const ae_int_t n,const bool isupper,const complex_1d_array & b,ae_int_t & info,densesolverreport & rep,complex_1d_array & x)1134 void hpdmatrixcholeskysolve(const complex_2d_array &cha, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, densesolverreport &rep, complex_1d_array &x)
1135 {
1136 alglib_impl::ae_state _alglib_env_state;
1137 alglib_impl::ae_state_init(&_alglib_env_state);
1138 try
1139 {
1140 alglib_impl::hpdmatrixcholeskysolve(const_cast<alglib_impl::ae_matrix*>(cha.c_ptr()), n, isupper, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), &info, const_cast<alglib_impl::densesolverreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
1141 alglib_impl::ae_state_clear(&_alglib_env_state);
1142 return;
1143 }
1144 catch(alglib_impl::ae_error_type)
1145 {
1146 throw ap_error(_alglib_env_state.error_msg);
1147 }
1148 catch(...)
1149 {
1150 throw;
1151 }
1152 }
1153
1154 /*************************************************************************
1155 Dense solver.
1156
1157 This subroutine finds solution of the linear system A*X=B with non-square,
1158 possibly degenerate A. System is solved in the least squares sense, and
1159 general least squares solution X = X0 + CX*y which minimizes |A*X-B| is
1160 returned. If A is non-degenerate, solution in the usual sense is returned
1161
1162 Algorithm features:
1163 * automatic detection of degenerate cases
1164 * iterative refinement
1165 * O(N^3) complexity
1166
1167 INPUT PARAMETERS
1168 A - array[0..NRows-1,0..NCols-1], system matrix
1169 NRows - vertical size of A
1170 NCols - horizontal size of A
1171 B - array[0..NCols-1], right part
1172 Threshold- a number in [0,1]. Singular values beyond Threshold are
1173 considered zero. Set it to 0.0, if you don't understand
1174 what it means, so the solver will choose good value on its
1175 own.
1176
1177 OUTPUT PARAMETERS
1178 Info - return code:
1179 * -4 SVD subroutine failed
1180 * -1 if NRows<=0 or NCols<=0 or Threshold<0 was passed
1181 * 1 if task is solved
1182 Rep - solver report, see below for more info
1183 X - array[0..N-1,0..M-1], it contains:
1184 * solution of A*X=B if A is non-singular (well-conditioned
1185 or ill-conditioned, but not very close to singular)
1186 * zeros, if A is singular or VERY close to singular
1187 (in this case Info=-3).
1188
1189 SOLVER REPORT
1190
1191 Subroutine sets following fields of the Rep structure:
1192 * R2 reciprocal of condition number: 1/cond(A), 2-norm.
1193 * N = NCols
1194 * K dim(Null(A))
1195 * CX array[0..N-1,0..K-1], kernel of A.
1196 Columns of CX store such vectors that A*CX[i]=0.
1197
1198 -- ALGLIB --
1199 Copyright 24.08.2009 by Bochkanov Sergey
1200 *************************************************************************/
rmatrixsolvels(const real_2d_array & a,const ae_int_t nrows,const ae_int_t ncols,const real_1d_array & b,const double threshold,ae_int_t & info,densesolverlsreport & rep,real_1d_array & x)1201 void rmatrixsolvels(const real_2d_array &a, const ae_int_t nrows, const ae_int_t ncols, const real_1d_array &b, const double threshold, ae_int_t &info, densesolverlsreport &rep, real_1d_array &x)
1202 {
1203 alglib_impl::ae_state _alglib_env_state;
1204 alglib_impl::ae_state_init(&_alglib_env_state);
1205 try
1206 {
1207 alglib_impl::rmatrixsolvels(const_cast<alglib_impl::ae_matrix*>(a.c_ptr()), nrows, ncols, const_cast<alglib_impl::ae_vector*>(b.c_ptr()), threshold, &info, const_cast<alglib_impl::densesolverlsreport*>(rep.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
1208 alglib_impl::ae_state_clear(&_alglib_env_state);
1209 return;
1210 }
1211 catch(alglib_impl::ae_error_type)
1212 {
1213 throw ap_error(_alglib_env_state.error_msg);
1214 }
1215 catch(...)
1216 {
1217 throw;
1218 }
1219 }
1220
1221 /*************************************************************************
1222
1223 *************************************************************************/
_nleqstate_owner()1224 _nleqstate_owner::_nleqstate_owner()
1225 {
1226 p_struct = (alglib_impl::nleqstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqstate), NULL);
1227 if( p_struct==NULL )
1228 throw ap_error("ALGLIB: malloc error");
1229 if( !alglib_impl::_nleqstate_init(p_struct, NULL, ae_false) )
1230 throw ap_error("ALGLIB: malloc error");
1231 }
1232
_nleqstate_owner(const _nleqstate_owner & rhs)1233 _nleqstate_owner::_nleqstate_owner(const _nleqstate_owner &rhs)
1234 {
1235 p_struct = (alglib_impl::nleqstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqstate), NULL);
1236 if( p_struct==NULL )
1237 throw ap_error("ALGLIB: malloc error");
1238 if( !alglib_impl::_nleqstate_init_copy(p_struct, const_cast<alglib_impl::nleqstate*>(rhs.p_struct), NULL, ae_false) )
1239 throw ap_error("ALGLIB: malloc error");
1240 }
1241
operator =(const _nleqstate_owner & rhs)1242 _nleqstate_owner& _nleqstate_owner::operator=(const _nleqstate_owner &rhs)
1243 {
1244 if( this==&rhs )
1245 return *this;
1246 alglib_impl::_nleqstate_clear(p_struct);
1247 if( !alglib_impl::_nleqstate_init_copy(p_struct, const_cast<alglib_impl::nleqstate*>(rhs.p_struct), NULL, ae_false) )
1248 throw ap_error("ALGLIB: malloc error");
1249 return *this;
1250 }
1251
~_nleqstate_owner()1252 _nleqstate_owner::~_nleqstate_owner()
1253 {
1254 alglib_impl::_nleqstate_clear(p_struct);
1255 ae_free(p_struct);
1256 }
1257
c_ptr()1258 alglib_impl::nleqstate* _nleqstate_owner::c_ptr()
1259 {
1260 return p_struct;
1261 }
1262
c_ptr() const1263 alglib_impl::nleqstate* _nleqstate_owner::c_ptr() const
1264 {
1265 return const_cast<alglib_impl::nleqstate*>(p_struct);
1266 }
nleqstate()1267 nleqstate::nleqstate() : _nleqstate_owner() ,needf(p_struct->needf),needfij(p_struct->needfij),xupdated(p_struct->xupdated),f(p_struct->f),fi(&p_struct->fi),j(&p_struct->j),x(&p_struct->x)
1268 {
1269 }
1270
nleqstate(const nleqstate & rhs)1271 nleqstate::nleqstate(const nleqstate &rhs):_nleqstate_owner(rhs) ,needf(p_struct->needf),needfij(p_struct->needfij),xupdated(p_struct->xupdated),f(p_struct->f),fi(&p_struct->fi),j(&p_struct->j),x(&p_struct->x)
1272 {
1273 }
1274
operator =(const nleqstate & rhs)1275 nleqstate& nleqstate::operator=(const nleqstate &rhs)
1276 {
1277 if( this==&rhs )
1278 return *this;
1279 _nleqstate_owner::operator=(rhs);
1280 return *this;
1281 }
1282
~nleqstate()1283 nleqstate::~nleqstate()
1284 {
1285 }
1286
1287
1288 /*************************************************************************
1289
1290 *************************************************************************/
_nleqreport_owner()1291 _nleqreport_owner::_nleqreport_owner()
1292 {
1293 p_struct = (alglib_impl::nleqreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqreport), NULL);
1294 if( p_struct==NULL )
1295 throw ap_error("ALGLIB: malloc error");
1296 if( !alglib_impl::_nleqreport_init(p_struct, NULL, ae_false) )
1297 throw ap_error("ALGLIB: malloc error");
1298 }
1299
_nleqreport_owner(const _nleqreport_owner & rhs)1300 _nleqreport_owner::_nleqreport_owner(const _nleqreport_owner &rhs)
1301 {
1302 p_struct = (alglib_impl::nleqreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::nleqreport), NULL);
1303 if( p_struct==NULL )
1304 throw ap_error("ALGLIB: malloc error");
1305 if( !alglib_impl::_nleqreport_init_copy(p_struct, const_cast<alglib_impl::nleqreport*>(rhs.p_struct), NULL, ae_false) )
1306 throw ap_error("ALGLIB: malloc error");
1307 }
1308
operator =(const _nleqreport_owner & rhs)1309 _nleqreport_owner& _nleqreport_owner::operator=(const _nleqreport_owner &rhs)
1310 {
1311 if( this==&rhs )
1312 return *this;
1313 alglib_impl::_nleqreport_clear(p_struct);
1314 if( !alglib_impl::_nleqreport_init_copy(p_struct, const_cast<alglib_impl::nleqreport*>(rhs.p_struct), NULL, ae_false) )
1315 throw ap_error("ALGLIB: malloc error");
1316 return *this;
1317 }
1318
~_nleqreport_owner()1319 _nleqreport_owner::~_nleqreport_owner()
1320 {
1321 alglib_impl::_nleqreport_clear(p_struct);
1322 ae_free(p_struct);
1323 }
1324
c_ptr()1325 alglib_impl::nleqreport* _nleqreport_owner::c_ptr()
1326 {
1327 return p_struct;
1328 }
1329
c_ptr() const1330 alglib_impl::nleqreport* _nleqreport_owner::c_ptr() const
1331 {
1332 return const_cast<alglib_impl::nleqreport*>(p_struct);
1333 }
nleqreport()1334 nleqreport::nleqreport() : _nleqreport_owner() ,iterationscount(p_struct->iterationscount),nfunc(p_struct->nfunc),njac(p_struct->njac),terminationtype(p_struct->terminationtype)
1335 {
1336 }
1337
nleqreport(const nleqreport & rhs)1338 nleqreport::nleqreport(const nleqreport &rhs):_nleqreport_owner(rhs) ,iterationscount(p_struct->iterationscount),nfunc(p_struct->nfunc),njac(p_struct->njac),terminationtype(p_struct->terminationtype)
1339 {
1340 }
1341
operator =(const nleqreport & rhs)1342 nleqreport& nleqreport::operator=(const nleqreport &rhs)
1343 {
1344 if( this==&rhs )
1345 return *this;
1346 _nleqreport_owner::operator=(rhs);
1347 return *this;
1348 }
1349
~nleqreport()1350 nleqreport::~nleqreport()
1351 {
1352 }
1353
1354 /*************************************************************************
1355 LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER
1356
1357 DESCRIPTION:
1358 This algorithm solves system of nonlinear equations
1359 F[0](x[0], ..., x[n-1]) = 0
1360 F[1](x[0], ..., x[n-1]) = 0
1361 ...
1362 F[M-1](x[0], ..., x[n-1]) = 0
1363 with M/N do not necessarily coincide. Algorithm converges quadratically
1364 under following conditions:
1365 * the solution set XS is nonempty
1366 * for some xs in XS there exist such neighbourhood N(xs) that:
1367 * vector function F(x) and its Jacobian J(x) are continuously
1368 differentiable on N
1369 * ||F(x)|| provides local error bound on N, i.e. there exists such
1370 c1, that ||F(x)||>c1*distance(x,XS)
1371 Note that these conditions are much more weaker than usual non-singularity
1372 conditions. For example, algorithm will converge for any affine function
1373 F (whether its Jacobian singular or not).
1374
1375
1376 REQUIREMENTS:
1377 Algorithm will request following information during its operation:
1378 * function vector F[] and Jacobian matrix at given point X
1379 * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X
1380
1381
1382 USAGE:
1383 1. User initializes algorithm state with NLEQCreateLM() call
1384 2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and
1385 other functions
1386 3. User calls NLEQSolve() function which takes algorithm state and
1387 pointers (delegates, etc.) to callback functions which calculate merit
1388 function value and Jacobian.
1389 4. User calls NLEQResults() to get solution
1390 5. Optionally, user may call NLEQRestartFrom() to solve another problem
1391 with same parameters (N/M) but another starting point and/or another
1392 function vector. NLEQRestartFrom() allows to reuse already initialized
1393 structure.
1394
1395
1396 INPUT PARAMETERS:
1397 N - space dimension, N>1:
1398 * if provided, only leading N elements of X are used
1399 * if not provided, determined automatically from size of X
1400 M - system size
1401 X - starting point
1402
1403
1404 OUTPUT PARAMETERS:
1405 State - structure which stores algorithm state
1406
1407
1408 NOTES:
1409 1. you may tune stopping conditions with NLEQSetCond() function
1410 2. if target function contains exp() or other fast growing functions, and
1411 optimization algorithm makes too large steps which leads to overflow,
1412 use NLEQSetStpMax() function to bound algorithm's steps.
1413 3. this algorithm is a slightly modified implementation of the method
1414 described in 'Levenberg-Marquardt method for constrained nonlinear
1415 equations with strong local convergence properties' by Christian Kanzow
1416 Nobuo Yamashita and Masao Fukushima and further developed in 'On the
1417 convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and
1418 Ya-Xiang Yuan.
1419
1420
1421 -- ALGLIB --
1422 Copyright 20.08.2009 by Bochkanov Sergey
1423 *************************************************************************/
nleqcreatelm(const ae_int_t n,const ae_int_t m,const real_1d_array & x,nleqstate & state)1424 void nleqcreatelm(const ae_int_t n, const ae_int_t m, const real_1d_array &x, nleqstate &state)
1425 {
1426 alglib_impl::ae_state _alglib_env_state;
1427 alglib_impl::ae_state_init(&_alglib_env_state);
1428 try
1429 {
1430 alglib_impl::nleqcreatelm(n, m, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
1431 alglib_impl::ae_state_clear(&_alglib_env_state);
1432 return;
1433 }
1434 catch(alglib_impl::ae_error_type)
1435 {
1436 throw ap_error(_alglib_env_state.error_msg);
1437 }
1438 catch(...)
1439 {
1440 throw;
1441 }
1442 }
1443
1444 /*************************************************************************
1445 LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER
1446
1447 DESCRIPTION:
1448 This algorithm solves system of nonlinear equations
1449 F[0](x[0], ..., x[n-1]) = 0
1450 F[1](x[0], ..., x[n-1]) = 0
1451 ...
1452 F[M-1](x[0], ..., x[n-1]) = 0
1453 with M/N do not necessarily coincide. Algorithm converges quadratically
1454 under following conditions:
1455 * the solution set XS is nonempty
1456 * for some xs in XS there exist such neighbourhood N(xs) that:
1457 * vector function F(x) and its Jacobian J(x) are continuously
1458 differentiable on N
1459 * ||F(x)|| provides local error bound on N, i.e. there exists such
1460 c1, that ||F(x)||>c1*distance(x,XS)
1461 Note that these conditions are much more weaker than usual non-singularity
1462 conditions. For example, algorithm will converge for any affine function
1463 F (whether its Jacobian singular or not).
1464
1465
1466 REQUIREMENTS:
1467 Algorithm will request following information during its operation:
1468 * function vector F[] and Jacobian matrix at given point X
1469 * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X
1470
1471
1472 USAGE:
1473 1. User initializes algorithm state with NLEQCreateLM() call
1474 2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and
1475 other functions
1476 3. User calls NLEQSolve() function which takes algorithm state and
1477 pointers (delegates, etc.) to callback functions which calculate merit
1478 function value and Jacobian.
1479 4. User calls NLEQResults() to get solution
1480 5. Optionally, user may call NLEQRestartFrom() to solve another problem
1481 with same parameters (N/M) but another starting point and/or another
1482 function vector. NLEQRestartFrom() allows to reuse already initialized
1483 structure.
1484
1485
1486 INPUT PARAMETERS:
1487 N - space dimension, N>1:
1488 * if provided, only leading N elements of X are used
1489 * if not provided, determined automatically from size of X
1490 M - system size
1491 X - starting point
1492
1493
1494 OUTPUT PARAMETERS:
1495 State - structure which stores algorithm state
1496
1497
1498 NOTES:
1499 1. you may tune stopping conditions with NLEQSetCond() function
1500 2. if target function contains exp() or other fast growing functions, and
1501 optimization algorithm makes too large steps which leads to overflow,
1502 use NLEQSetStpMax() function to bound algorithm's steps.
1503 3. this algorithm is a slightly modified implementation of the method
1504 described in 'Levenberg-Marquardt method for constrained nonlinear
1505 equations with strong local convergence properties' by Christian Kanzow
1506 Nobuo Yamashita and Masao Fukushima and further developed in 'On the
1507 convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and
1508 Ya-Xiang Yuan.
1509
1510
1511 -- ALGLIB --
1512 Copyright 20.08.2009 by Bochkanov Sergey
1513 *************************************************************************/
nleqcreatelm(const ae_int_t m,const real_1d_array & x,nleqstate & state)1514 void nleqcreatelm(const ae_int_t m, const real_1d_array &x, nleqstate &state)
1515 {
1516 alglib_impl::ae_state _alglib_env_state;
1517 ae_int_t n;
1518
1519 n = x.length();
1520 alglib_impl::ae_state_init(&_alglib_env_state);
1521 try
1522 {
1523 alglib_impl::nleqcreatelm(n, m, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
1524
1525 alglib_impl::ae_state_clear(&_alglib_env_state);
1526 return;
1527 }
1528 catch(alglib_impl::ae_error_type)
1529 {
1530 throw ap_error(_alglib_env_state.error_msg);
1531 }
1532 catch(...)
1533 {
1534 throw;
1535 }
1536 }
1537
1538 /*************************************************************************
1539 This function sets stopping conditions for the nonlinear solver
1540
1541 INPUT PARAMETERS:
1542 State - structure which stores algorithm state
1543 EpsF - >=0
1544 The subroutine finishes its work if on k+1-th iteration
1545 the condition ||F||<=EpsF is satisfied
1546 MaxIts - maximum number of iterations. If MaxIts=0, the number of
1547 iterations is unlimited.
1548
1549 Passing EpsF=0 and MaxIts=0 simultaneously will lead to automatic
1550 stopping criterion selection (small EpsF).
1551
1552 NOTES:
1553
1554 -- ALGLIB --
1555 Copyright 20.08.2010 by Bochkanov Sergey
1556 *************************************************************************/
nleqsetcond(const nleqstate & state,const double epsf,const ae_int_t maxits)1557 void nleqsetcond(const nleqstate &state, const double epsf, const ae_int_t maxits)
1558 {
1559 alglib_impl::ae_state _alglib_env_state;
1560 alglib_impl::ae_state_init(&_alglib_env_state);
1561 try
1562 {
1563 alglib_impl::nleqsetcond(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), epsf, maxits, &_alglib_env_state);
1564 alglib_impl::ae_state_clear(&_alglib_env_state);
1565 return;
1566 }
1567 catch(alglib_impl::ae_error_type)
1568 {
1569 throw ap_error(_alglib_env_state.error_msg);
1570 }
1571 catch(...)
1572 {
1573 throw;
1574 }
1575 }
1576
1577 /*************************************************************************
1578 This function turns on/off reporting.
1579
1580 INPUT PARAMETERS:
1581 State - structure which stores algorithm state
1582 NeedXRep- whether iteration reports are needed or not
1583
1584 If NeedXRep is True, algorithm will call rep() callback function if it is
1585 provided to NLEQSolve().
1586
1587 -- ALGLIB --
1588 Copyright 20.08.2010 by Bochkanov Sergey
1589 *************************************************************************/
nleqsetxrep(const nleqstate & state,const bool needxrep)1590 void nleqsetxrep(const nleqstate &state, const bool needxrep)
1591 {
1592 alglib_impl::ae_state _alglib_env_state;
1593 alglib_impl::ae_state_init(&_alglib_env_state);
1594 try
1595 {
1596 alglib_impl::nleqsetxrep(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), needxrep, &_alglib_env_state);
1597 alglib_impl::ae_state_clear(&_alglib_env_state);
1598 return;
1599 }
1600 catch(alglib_impl::ae_error_type)
1601 {
1602 throw ap_error(_alglib_env_state.error_msg);
1603 }
1604 catch(...)
1605 {
1606 throw;
1607 }
1608 }
1609
1610 /*************************************************************************
1611 This function sets maximum step length
1612
1613 INPUT PARAMETERS:
1614 State - structure which stores algorithm state
1615 StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
1616 want to limit step length.
1617
1618 Use this subroutine when target function contains exp() or other fast
1619 growing functions, and algorithm makes too large steps which lead to
1620 overflow. This function allows us to reject steps that are too large (and
1621 therefore expose us to the possible overflow) without actually calculating
1622 function value at the x+stp*d.
1623
1624 -- ALGLIB --
1625 Copyright 20.08.2010 by Bochkanov Sergey
1626 *************************************************************************/
nleqsetstpmax(const nleqstate & state,const double stpmax)1627 void nleqsetstpmax(const nleqstate &state, const double stpmax)
1628 {
1629 alglib_impl::ae_state _alglib_env_state;
1630 alglib_impl::ae_state_init(&_alglib_env_state);
1631 try
1632 {
1633 alglib_impl::nleqsetstpmax(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), stpmax, &_alglib_env_state);
1634 alglib_impl::ae_state_clear(&_alglib_env_state);
1635 return;
1636 }
1637 catch(alglib_impl::ae_error_type)
1638 {
1639 throw ap_error(_alglib_env_state.error_msg);
1640 }
1641 catch(...)
1642 {
1643 throw;
1644 }
1645 }
1646
1647 /*************************************************************************
1648 This function provides reverse communication interface
1649 Reverse communication interface is not documented or recommended to use.
1650 See below for functions which provide better documented API
1651 *************************************************************************/
nleqiteration(const nleqstate & state)1652 bool nleqiteration(const nleqstate &state)
1653 {
1654 alglib_impl::ae_state _alglib_env_state;
1655 alglib_impl::ae_state_init(&_alglib_env_state);
1656 try
1657 {
1658 ae_bool result = alglib_impl::nleqiteration(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), &_alglib_env_state);
1659 alglib_impl::ae_state_clear(&_alglib_env_state);
1660 return *(reinterpret_cast<bool*>(&result));
1661 }
1662 catch(alglib_impl::ae_error_type)
1663 {
1664 throw ap_error(_alglib_env_state.error_msg);
1665 }
1666 catch(...)
1667 {
1668 throw;
1669 }
1670 }
1671
1672
nleqsolve(nleqstate & state,void (* func)(const real_1d_array & x,double & func,void * ptr),void (* jac)(const real_1d_array & x,real_1d_array & fi,real_2d_array & jac,void * ptr),void (* rep)(const real_1d_array & x,double func,void * ptr),void * ptr)1673 void nleqsolve(nleqstate &state,
1674 void (*func)(const real_1d_array &x, double &func, void *ptr),
1675 void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &jac, void *ptr),
1676 void (*rep)(const real_1d_array &x, double func, void *ptr),
1677 void *ptr)
1678 {
1679 alglib_impl::ae_state _alglib_env_state;
1680 if( func==NULL )
1681 throw ap_error("ALGLIB: error in 'nleqsolve()' (func is NULL)");
1682 if( jac==NULL )
1683 throw ap_error("ALGLIB: error in 'nleqsolve()' (jac is NULL)");
1684 alglib_impl::ae_state_init(&_alglib_env_state);
1685 try
1686 {
1687 while( alglib_impl::nleqiteration(state.c_ptr(), &_alglib_env_state) )
1688 {
1689 if( state.needf )
1690 {
1691 func(state.x, state.f, ptr);
1692 continue;
1693 }
1694 if( state.needfij )
1695 {
1696 jac(state.x, state.fi, state.j, ptr);
1697 continue;
1698 }
1699 if( state.xupdated )
1700 {
1701 if( rep!=NULL )
1702 rep(state.x, state.f, ptr);
1703 continue;
1704 }
1705 throw ap_error("ALGLIB: error in 'nleqsolve' (some derivatives were not provided?)");
1706 }
1707 alglib_impl::ae_state_clear(&_alglib_env_state);
1708 }
1709 catch(alglib_impl::ae_error_type)
1710 {
1711 throw ap_error(_alglib_env_state.error_msg);
1712 }
1713 catch(...)
1714 {
1715 throw;
1716 }
1717 }
1718
1719
1720
1721 /*************************************************************************
1722 NLEQ solver results
1723
1724 INPUT PARAMETERS:
1725 State - algorithm state.
1726
1727 OUTPUT PARAMETERS:
1728 X - array[0..N-1], solution
1729 Rep - optimization report:
1730 * Rep.TerminationType completetion code:
1731 * -4 ERROR: algorithm has converged to the
1732 stationary point Xf which is local minimum of
1733 f=F[0]^2+...+F[m-1]^2, but is not solution of
1734 nonlinear system.
1735 * 1 sqrt(f)<=EpsF.
1736 * 5 MaxIts steps was taken
1737 * 7 stopping conditions are too stringent,
1738 further improvement is impossible
1739 * Rep.IterationsCount contains iterations count
1740 * NFEV countains number of function calculations
1741 * ActiveConstraints contains number of active constraints
1742
1743 -- ALGLIB --
1744 Copyright 20.08.2009 by Bochkanov Sergey
1745 *************************************************************************/
nleqresults(const nleqstate & state,real_1d_array & x,nleqreport & rep)1746 void nleqresults(const nleqstate &state, real_1d_array &x, nleqreport &rep)
1747 {
1748 alglib_impl::ae_state _alglib_env_state;
1749 alglib_impl::ae_state_init(&_alglib_env_state);
1750 try
1751 {
1752 alglib_impl::nleqresults(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqreport*>(rep.c_ptr()), &_alglib_env_state);
1753 alglib_impl::ae_state_clear(&_alglib_env_state);
1754 return;
1755 }
1756 catch(alglib_impl::ae_error_type)
1757 {
1758 throw ap_error(_alglib_env_state.error_msg);
1759 }
1760 catch(...)
1761 {
1762 throw;
1763 }
1764 }
1765
1766 /*************************************************************************
1767 NLEQ solver results
1768
1769 Buffered implementation of NLEQResults(), which uses pre-allocated buffer
1770 to store X[]. If buffer size is too small, it resizes buffer. It is
1771 intended to be used in the inner cycles of performance critical algorithms
1772 where array reallocation penalty is too large to be ignored.
1773
1774 -- ALGLIB --
1775 Copyright 20.08.2009 by Bochkanov Sergey
1776 *************************************************************************/
nleqresultsbuf(const nleqstate & state,real_1d_array & x,nleqreport & rep)1777 void nleqresultsbuf(const nleqstate &state, real_1d_array &x, nleqreport &rep)
1778 {
1779 alglib_impl::ae_state _alglib_env_state;
1780 alglib_impl::ae_state_init(&_alglib_env_state);
1781 try
1782 {
1783 alglib_impl::nleqresultsbuf(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::nleqreport*>(rep.c_ptr()), &_alglib_env_state);
1784 alglib_impl::ae_state_clear(&_alglib_env_state);
1785 return;
1786 }
1787 catch(alglib_impl::ae_error_type)
1788 {
1789 throw ap_error(_alglib_env_state.error_msg);
1790 }
1791 catch(...)
1792 {
1793 throw;
1794 }
1795 }
1796
1797 /*************************************************************************
1798 This subroutine restarts CG algorithm from new point. All optimization
1799 parameters are left unchanged.
1800
1801 This function allows to solve multiple optimization problems (which
1802 must have same number of dimensions) without object reallocation penalty.
1803
1804 INPUT PARAMETERS:
1805 State - structure used for reverse communication previously
1806 allocated with MinCGCreate call.
1807 X - new starting point.
1808 BndL - new lower bounds
1809 BndU - new upper bounds
1810
1811 -- ALGLIB --
1812 Copyright 30.07.2010 by Bochkanov Sergey
1813 *************************************************************************/
nleqrestartfrom(const nleqstate & state,const real_1d_array & x)1814 void nleqrestartfrom(const nleqstate &state, const real_1d_array &x)
1815 {
1816 alglib_impl::ae_state _alglib_env_state;
1817 alglib_impl::ae_state_init(&_alglib_env_state);
1818 try
1819 {
1820 alglib_impl::nleqrestartfrom(const_cast<alglib_impl::nleqstate*>(state.c_ptr()), const_cast<alglib_impl::ae_vector*>(x.c_ptr()), &_alglib_env_state);
1821 alglib_impl::ae_state_clear(&_alglib_env_state);
1822 return;
1823 }
1824 catch(alglib_impl::ae_error_type)
1825 {
1826 throw ap_error(_alglib_env_state.error_msg);
1827 }
1828 catch(...)
1829 {
1830 throw;
1831 }
1832 }
1833 }
1834
1835 /////////////////////////////////////////////////////////////////////////
1836 //
1837 // THIS SECTION CONTAINS IMPLEMENTATION OF COMPUTATIONAL CORE
1838 //
1839 /////////////////////////////////////////////////////////////////////////
1840 namespace alglib_impl
1841 {
1842 static void densesolver_rmatrixlusolveinternal(/* Real */ ae_matrix* lua,
1843 /* Integer */ ae_vector* p,
1844 double scalea,
1845 ae_int_t n,
1846 /* Real */ ae_matrix* a,
1847 ae_bool havea,
1848 /* Real */ ae_matrix* b,
1849 ae_int_t m,
1850 ae_int_t* info,
1851 densesolverreport* rep,
1852 /* Real */ ae_matrix* x,
1853 ae_state *_state);
1854 static void densesolver_spdmatrixcholeskysolveinternal(/* Real */ ae_matrix* cha,
1855 double sqrtscalea,
1856 ae_int_t n,
1857 ae_bool isupper,
1858 /* Real */ ae_matrix* a,
1859 ae_bool havea,
1860 /* Real */ ae_matrix* b,
1861 ae_int_t m,
1862 ae_int_t* info,
1863 densesolverreport* rep,
1864 /* Real */ ae_matrix* x,
1865 ae_state *_state);
1866 static void densesolver_cmatrixlusolveinternal(/* Complex */ ae_matrix* lua,
1867 /* Integer */ ae_vector* p,
1868 double scalea,
1869 ae_int_t n,
1870 /* Complex */ ae_matrix* a,
1871 ae_bool havea,
1872 /* Complex */ ae_matrix* b,
1873 ae_int_t m,
1874 ae_int_t* info,
1875 densesolverreport* rep,
1876 /* Complex */ ae_matrix* x,
1877 ae_state *_state);
1878 static void densesolver_hpdmatrixcholeskysolveinternal(/* Complex */ ae_matrix* cha,
1879 double sqrtscalea,
1880 ae_int_t n,
1881 ae_bool isupper,
1882 /* Complex */ ae_matrix* a,
1883 ae_bool havea,
1884 /* Complex */ ae_matrix* b,
1885 ae_int_t m,
1886 ae_int_t* info,
1887 densesolverreport* rep,
1888 /* Complex */ ae_matrix* x,
1889 ae_state *_state);
1890 static ae_int_t densesolver_densesolverrfsmax(ae_int_t n,
1891 double r1,
1892 double rinf,
1893 ae_state *_state);
1894 static ae_int_t densesolver_densesolverrfsmaxv2(ae_int_t n,
1895 double r2,
1896 ae_state *_state);
1897 static void densesolver_rbasiclusolve(/* Real */ ae_matrix* lua,
1898 /* Integer */ ae_vector* p,
1899 double scalea,
1900 ae_int_t n,
1901 /* Real */ ae_vector* xb,
1902 /* Real */ ae_vector* tmp,
1903 ae_state *_state);
1904 static void densesolver_spdbasiccholeskysolve(/* Real */ ae_matrix* cha,
1905 double sqrtscalea,
1906 ae_int_t n,
1907 ae_bool isupper,
1908 /* Real */ ae_vector* xb,
1909 /* Real */ ae_vector* tmp,
1910 ae_state *_state);
1911 static void densesolver_cbasiclusolve(/* Complex */ ae_matrix* lua,
1912 /* Integer */ ae_vector* p,
1913 double scalea,
1914 ae_int_t n,
1915 /* Complex */ ae_vector* xb,
1916 /* Complex */ ae_vector* tmp,
1917 ae_state *_state);
1918 static void densesolver_hpdbasiccholeskysolve(/* Complex */ ae_matrix* cha,
1919 double sqrtscalea,
1920 ae_int_t n,
1921 ae_bool isupper,
1922 /* Complex */ ae_vector* xb,
1923 /* Complex */ ae_vector* tmp,
1924 ae_state *_state);
1925
1926
1927 static ae_int_t nleq_armijomaxfev = 20;
1928 static void nleq_clearrequestfields(nleqstate* state, ae_state *_state);
1929 static ae_bool nleq_increaselambda(double* lambdav,
1930 double* nu,
1931 double lambdaup,
1932 ae_state *_state);
1933 static void nleq_decreaselambda(double* lambdav,
1934 double* nu,
1935 double lambdadown,
1936 ae_state *_state);
1937
1938
1939
1940
1941
1942 /*************************************************************************
1943 Dense solver.
1944
1945 This subroutine solves a system A*x=b, where A is NxN non-denegerate
1946 real matrix, x and b are vectors.
1947
1948 Algorithm features:
1949 * automatic detection of degenerate cases
1950 * condition number estimation
1951 * iterative refinement
1952 * O(N^3) complexity
1953
1954 INPUT PARAMETERS
1955 A - array[0..N-1,0..N-1], system matrix
1956 N - size of A
1957 B - array[0..N-1], right part
1958
1959 OUTPUT PARAMETERS
1960 Info - return code:
1961 * -3 A is singular, or VERY close to singular.
1962 X is filled by zeros in such cases.
1963 * -1 N<=0 was passed
1964 * 1 task is solved (but matrix A may be ill-conditioned,
1965 check R1/RInf parameters for condition numbers).
1966 Rep - solver report, see below for more info
1967 X - array[0..N-1], it contains:
1968 * solution of A*x=b if A is non-singular (well-conditioned
1969 or ill-conditioned, but not very close to singular)
1970 * zeros, if A is singular or VERY close to singular
1971 (in this case Info=-3).
1972
1973 SOLVER REPORT
1974
1975 Subroutine sets following fields of the Rep structure:
1976 * R1 reciprocal of condition number: 1/cond(A), 1-norm.
1977 * RInf reciprocal of condition number: 1/cond(A), inf-norm.
1978
1979 -- ALGLIB --
1980 Copyright 27.01.2010 by Bochkanov Sergey
1981 *************************************************************************/
rmatrixsolve(ae_matrix * a,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)1982 void rmatrixsolve(/* Real */ ae_matrix* a,
1983 ae_int_t n,
1984 /* Real */ ae_vector* b,
1985 ae_int_t* info,
1986 densesolverreport* rep,
1987 /* Real */ ae_vector* x,
1988 ae_state *_state)
1989 {
1990 ae_frame _frame_block;
1991 ae_matrix bm;
1992 ae_matrix xm;
1993
1994 ae_frame_make(_state, &_frame_block);
1995 *info = 0;
1996 _densesolverreport_clear(rep);
1997 ae_vector_clear(x);
1998 ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
1999 ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
2000
2001 if( n<=0 )
2002 {
2003 *info = -1;
2004 ae_frame_leave(_state);
2005 return;
2006 }
2007 ae_matrix_set_length(&bm, n, 1, _state);
2008 ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2009 rmatrixsolvem(a, n, &bm, 1, ae_true, info, rep, &xm, _state);
2010 ae_vector_set_length(x, n, _state);
2011 ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
2012 ae_frame_leave(_state);
2013 }
2014
2015
2016 /*************************************************************************
2017 Dense solver.
2018
2019 Similar to RMatrixSolve() but solves task with multiple right parts (where
2020 b and x are NxM matrices).
2021
2022 Algorithm features:
2023 * automatic detection of degenerate cases
2024 * condition number estimation
2025 * optional iterative refinement
2026 * O(N^3+M*N^2) complexity
2027
2028 INPUT PARAMETERS
2029 A - array[0..N-1,0..N-1], system matrix
2030 N - size of A
2031 B - array[0..N-1,0..M-1], right part
2032 M - right part size
2033 RFS - iterative refinement switch:
2034 * True - refinement is used.
2035 Less performance, more precision.
2036 * False - refinement is not used.
2037 More performance, less precision.
2038
2039 OUTPUT PARAMETERS
2040 Info - same as in RMatrixSolve
2041 Rep - same as in RMatrixSolve
2042 X - same as in RMatrixSolve
2043
2044 -- ALGLIB --
2045 Copyright 27.01.2010 by Bochkanov Sergey
2046 *************************************************************************/
rmatrixsolvem(ae_matrix * a,ae_int_t n,ae_matrix * b,ae_int_t m,ae_bool rfs,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2047 void rmatrixsolvem(/* Real */ ae_matrix* a,
2048 ae_int_t n,
2049 /* Real */ ae_matrix* b,
2050 ae_int_t m,
2051 ae_bool rfs,
2052 ae_int_t* info,
2053 densesolverreport* rep,
2054 /* Real */ ae_matrix* x,
2055 ae_state *_state)
2056 {
2057 ae_frame _frame_block;
2058 ae_matrix da;
2059 ae_matrix emptya;
2060 ae_vector p;
2061 double scalea;
2062 ae_int_t i;
2063 ae_int_t j;
2064
2065 ae_frame_make(_state, &_frame_block);
2066 *info = 0;
2067 _densesolverreport_clear(rep);
2068 ae_matrix_clear(x);
2069 ae_matrix_init(&da, 0, 0, DT_REAL, _state, ae_true);
2070 ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
2071 ae_vector_init(&p, 0, DT_INT, _state, ae_true);
2072
2073
2074 /*
2075 * prepare: check inputs, allocate space...
2076 */
2077 if( n<=0||m<=0 )
2078 {
2079 *info = -1;
2080 ae_frame_leave(_state);
2081 return;
2082 }
2083 ae_matrix_set_length(&da, n, n, _state);
2084
2085 /*
2086 * 1. scale matrix, max(|A[i,j]|)
2087 * 2. factorize scaled matrix
2088 * 3. solve
2089 */
2090 scalea = 0;
2091 for(i=0; i<=n-1; i++)
2092 {
2093 for(j=0; j<=n-1; j++)
2094 {
2095 scalea = ae_maxreal(scalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
2096 }
2097 }
2098 if( ae_fp_eq(scalea,0) )
2099 {
2100 scalea = 1;
2101 }
2102 scalea = 1/scalea;
2103 for(i=0; i<=n-1; i++)
2104 {
2105 ae_v_move(&da.ptr.pp_double[i][0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,n-1));
2106 }
2107 rmatrixlu(&da, n, n, &p, _state);
2108 if( rfs )
2109 {
2110 densesolver_rmatrixlusolveinternal(&da, &p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
2111 }
2112 else
2113 {
2114 densesolver_rmatrixlusolveinternal(&da, &p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
2115 }
2116 ae_frame_leave(_state);
2117 }
2118
2119
2120 /*************************************************************************
2121 Dense solver.
2122
2123 This subroutine solves a system A*X=B, where A is NxN non-denegerate
2124 real matrix given by its LU decomposition, X and B are NxM real matrices.
2125
2126 Algorithm features:
2127 * automatic detection of degenerate cases
2128 * O(N^2) complexity
2129 * condition number estimation
2130
2131 No iterative refinement is provided because exact form of original matrix
2132 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
2133 need iterative refinement.
2134
2135 INPUT PARAMETERS
2136 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
2137 P - array[0..N-1], pivots array, RMatrixLU result
2138 N - size of A
2139 B - array[0..N-1], right part
2140
2141 OUTPUT PARAMETERS
2142 Info - same as in RMatrixSolve
2143 Rep - same as in RMatrixSolve
2144 X - same as in RMatrixSolve
2145
2146 -- ALGLIB --
2147 Copyright 27.01.2010 by Bochkanov Sergey
2148 *************************************************************************/
rmatrixlusolve(ae_matrix * lua,ae_vector * p,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)2149 void rmatrixlusolve(/* Real */ ae_matrix* lua,
2150 /* Integer */ ae_vector* p,
2151 ae_int_t n,
2152 /* Real */ ae_vector* b,
2153 ae_int_t* info,
2154 densesolverreport* rep,
2155 /* Real */ ae_vector* x,
2156 ae_state *_state)
2157 {
2158 ae_frame _frame_block;
2159 ae_matrix bm;
2160 ae_matrix xm;
2161
2162 ae_frame_make(_state, &_frame_block);
2163 *info = 0;
2164 _densesolverreport_clear(rep);
2165 ae_vector_clear(x);
2166 ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
2167 ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
2168
2169 if( n<=0 )
2170 {
2171 *info = -1;
2172 ae_frame_leave(_state);
2173 return;
2174 }
2175 ae_matrix_set_length(&bm, n, 1, _state);
2176 ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2177 rmatrixlusolvem(lua, p, n, &bm, 1, info, rep, &xm, _state);
2178 ae_vector_set_length(x, n, _state);
2179 ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
2180 ae_frame_leave(_state);
2181 }
2182
2183
2184 /*************************************************************************
2185 Dense solver.
2186
2187 Similar to RMatrixLUSolve() but solves task with multiple right parts
2188 (where b and x are NxM matrices).
2189
2190 Algorithm features:
2191 * automatic detection of degenerate cases
2192 * O(M*N^2) complexity
2193 * condition number estimation
2194
2195 No iterative refinement is provided because exact form of original matrix
2196 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
2197 need iterative refinement.
2198
2199 INPUT PARAMETERS
2200 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
2201 P - array[0..N-1], pivots array, RMatrixLU result
2202 N - size of A
2203 B - array[0..N-1,0..M-1], right part
2204 M - right part size
2205
2206 OUTPUT PARAMETERS
2207 Info - same as in RMatrixSolve
2208 Rep - same as in RMatrixSolve
2209 X - same as in RMatrixSolve
2210
2211 -- ALGLIB --
2212 Copyright 27.01.2010 by Bochkanov Sergey
2213 *************************************************************************/
rmatrixlusolvem(ae_matrix * lua,ae_vector * p,ae_int_t n,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2214 void rmatrixlusolvem(/* Real */ ae_matrix* lua,
2215 /* Integer */ ae_vector* p,
2216 ae_int_t n,
2217 /* Real */ ae_matrix* b,
2218 ae_int_t m,
2219 ae_int_t* info,
2220 densesolverreport* rep,
2221 /* Real */ ae_matrix* x,
2222 ae_state *_state)
2223 {
2224 ae_frame _frame_block;
2225 ae_matrix emptya;
2226 ae_int_t i;
2227 ae_int_t j;
2228 double scalea;
2229
2230 ae_frame_make(_state, &_frame_block);
2231 *info = 0;
2232 _densesolverreport_clear(rep);
2233 ae_matrix_clear(x);
2234 ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
2235
2236
2237 /*
2238 * prepare: check inputs, allocate space...
2239 */
2240 if( n<=0||m<=0 )
2241 {
2242 *info = -1;
2243 ae_frame_leave(_state);
2244 return;
2245 }
2246
2247 /*
2248 * 1. scale matrix, max(|U[i,j]|)
2249 * we assume that LU is in its normal form, i.e. |L[i,j]|<=1
2250 * 2. solve
2251 */
2252 scalea = 0;
2253 for(i=0; i<=n-1; i++)
2254 {
2255 for(j=i; j<=n-1; j++)
2256 {
2257 scalea = ae_maxreal(scalea, ae_fabs(lua->ptr.pp_double[i][j], _state), _state);
2258 }
2259 }
2260 if( ae_fp_eq(scalea,0) )
2261 {
2262 scalea = 1;
2263 }
2264 scalea = 1/scalea;
2265 densesolver_rmatrixlusolveinternal(lua, p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
2266 ae_frame_leave(_state);
2267 }
2268
2269
2270 /*************************************************************************
2271 Dense solver.
2272
2273 This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS
2274 LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have
2275 both A and its LU decomposition.
2276
2277 Algorithm features:
2278 * automatic detection of degenerate cases
2279 * condition number estimation
2280 * iterative refinement
2281 * O(N^2) complexity
2282
2283 INPUT PARAMETERS
2284 A - array[0..N-1,0..N-1], system matrix
2285 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
2286 P - array[0..N-1], pivots array, RMatrixLU result
2287 N - size of A
2288 B - array[0..N-1], right part
2289
2290 OUTPUT PARAMETERS
2291 Info - same as in RMatrixSolveM
2292 Rep - same as in RMatrixSolveM
2293 X - same as in RMatrixSolveM
2294
2295 -- ALGLIB --
2296 Copyright 27.01.2010 by Bochkanov Sergey
2297 *************************************************************************/
rmatrixmixedsolve(ae_matrix * a,ae_matrix * lua,ae_vector * p,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)2298 void rmatrixmixedsolve(/* Real */ ae_matrix* a,
2299 /* Real */ ae_matrix* lua,
2300 /* Integer */ ae_vector* p,
2301 ae_int_t n,
2302 /* Real */ ae_vector* b,
2303 ae_int_t* info,
2304 densesolverreport* rep,
2305 /* Real */ ae_vector* x,
2306 ae_state *_state)
2307 {
2308 ae_frame _frame_block;
2309 ae_matrix bm;
2310 ae_matrix xm;
2311
2312 ae_frame_make(_state, &_frame_block);
2313 *info = 0;
2314 _densesolverreport_clear(rep);
2315 ae_vector_clear(x);
2316 ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
2317 ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
2318
2319 if( n<=0 )
2320 {
2321 *info = -1;
2322 ae_frame_leave(_state);
2323 return;
2324 }
2325 ae_matrix_set_length(&bm, n, 1, _state);
2326 ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
2327 rmatrixmixedsolvem(a, lua, p, n, &bm, 1, info, rep, &xm, _state);
2328 ae_vector_set_length(x, n, _state);
2329 ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
2330 ae_frame_leave(_state);
2331 }
2332
2333
2334 /*************************************************************************
2335 Dense solver.
2336
2337 Similar to RMatrixMixedSolve() but solves task with multiple right parts
2338 (where b and x are NxM matrices).
2339
2340 Algorithm features:
2341 * automatic detection of degenerate cases
2342 * condition number estimation
2343 * iterative refinement
2344 * O(M*N^2) complexity
2345
2346 INPUT PARAMETERS
2347 A - array[0..N-1,0..N-1], system matrix
2348 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
2349 P - array[0..N-1], pivots array, RMatrixLU result
2350 N - size of A
2351 B - array[0..N-1,0..M-1], right part
2352 M - right part size
2353
2354 OUTPUT PARAMETERS
2355 Info - same as in RMatrixSolveM
2356 Rep - same as in RMatrixSolveM
2357 X - same as in RMatrixSolveM
2358
2359 -- ALGLIB --
2360 Copyright 27.01.2010 by Bochkanov Sergey
2361 *************************************************************************/
rmatrixmixedsolvem(ae_matrix * a,ae_matrix * lua,ae_vector * p,ae_int_t n,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2362 void rmatrixmixedsolvem(/* Real */ ae_matrix* a,
2363 /* Real */ ae_matrix* lua,
2364 /* Integer */ ae_vector* p,
2365 ae_int_t n,
2366 /* Real */ ae_matrix* b,
2367 ae_int_t m,
2368 ae_int_t* info,
2369 densesolverreport* rep,
2370 /* Real */ ae_matrix* x,
2371 ae_state *_state)
2372 {
2373 double scalea;
2374 ae_int_t i;
2375 ae_int_t j;
2376
2377 *info = 0;
2378 _densesolverreport_clear(rep);
2379 ae_matrix_clear(x);
2380
2381
2382 /*
2383 * prepare: check inputs, allocate space...
2384 */
2385 if( n<=0||m<=0 )
2386 {
2387 *info = -1;
2388 return;
2389 }
2390
2391 /*
2392 * 1. scale matrix, max(|A[i,j]|)
2393 * 2. factorize scaled matrix
2394 * 3. solve
2395 */
2396 scalea = 0;
2397 for(i=0; i<=n-1; i++)
2398 {
2399 for(j=0; j<=n-1; j++)
2400 {
2401 scalea = ae_maxreal(scalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
2402 }
2403 }
2404 if( ae_fp_eq(scalea,0) )
2405 {
2406 scalea = 1;
2407 }
2408 scalea = 1/scalea;
2409 densesolver_rmatrixlusolveinternal(lua, p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
2410 }
2411
2412
2413 /*************************************************************************
2414 Dense solver. Same as RMatrixSolveM(), but for complex matrices.
2415
2416 Algorithm features:
2417 * automatic detection of degenerate cases
2418 * condition number estimation
2419 * iterative refinement
2420 * O(N^3+M*N^2) complexity
2421
2422 INPUT PARAMETERS
2423 A - array[0..N-1,0..N-1], system matrix
2424 N - size of A
2425 B - array[0..N-1,0..M-1], right part
2426 M - right part size
2427 RFS - iterative refinement switch:
2428 * True - refinement is used.
2429 Less performance, more precision.
2430 * False - refinement is not used.
2431 More performance, less precision.
2432
2433 OUTPUT PARAMETERS
2434 Info - same as in RMatrixSolve
2435 Rep - same as in RMatrixSolve
2436 X - same as in RMatrixSolve
2437
2438 -- ALGLIB --
2439 Copyright 27.01.2010 by Bochkanov Sergey
2440 *************************************************************************/
cmatrixsolvem(ae_matrix * a,ae_int_t n,ae_matrix * b,ae_int_t m,ae_bool rfs,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2441 void cmatrixsolvem(/* Complex */ ae_matrix* a,
2442 ae_int_t n,
2443 /* Complex */ ae_matrix* b,
2444 ae_int_t m,
2445 ae_bool rfs,
2446 ae_int_t* info,
2447 densesolverreport* rep,
2448 /* Complex */ ae_matrix* x,
2449 ae_state *_state)
2450 {
2451 ae_frame _frame_block;
2452 ae_matrix da;
2453 ae_matrix emptya;
2454 ae_vector p;
2455 double scalea;
2456 ae_int_t i;
2457 ae_int_t j;
2458
2459 ae_frame_make(_state, &_frame_block);
2460 *info = 0;
2461 _densesolverreport_clear(rep);
2462 ae_matrix_clear(x);
2463 ae_matrix_init(&da, 0, 0, DT_COMPLEX, _state, ae_true);
2464 ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
2465 ae_vector_init(&p, 0, DT_INT, _state, ae_true);
2466
2467
2468 /*
2469 * prepare: check inputs, allocate space...
2470 */
2471 if( n<=0||m<=0 )
2472 {
2473 *info = -1;
2474 ae_frame_leave(_state);
2475 return;
2476 }
2477 ae_matrix_set_length(&da, n, n, _state);
2478
2479 /*
2480 * 1. scale matrix, max(|A[i,j]|)
2481 * 2. factorize scaled matrix
2482 * 3. solve
2483 */
2484 scalea = 0;
2485 for(i=0; i<=n-1; i++)
2486 {
2487 for(j=0; j<=n-1; j++)
2488 {
2489 scalea = ae_maxreal(scalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
2490 }
2491 }
2492 if( ae_fp_eq(scalea,0) )
2493 {
2494 scalea = 1;
2495 }
2496 scalea = 1/scalea;
2497 for(i=0; i<=n-1; i++)
2498 {
2499 ae_v_cmove(&da.ptr.pp_complex[i][0], 1, &a->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,n-1));
2500 }
2501 cmatrixlu(&da, n, n, &p, _state);
2502 if( rfs )
2503 {
2504 densesolver_cmatrixlusolveinternal(&da, &p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
2505 }
2506 else
2507 {
2508 densesolver_cmatrixlusolveinternal(&da, &p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
2509 }
2510 ae_frame_leave(_state);
2511 }
2512
2513
2514 /*************************************************************************
2515 Dense solver. Same as RMatrixSolve(), but for complex matrices.
2516
2517 Algorithm features:
2518 * automatic detection of degenerate cases
2519 * condition number estimation
2520 * iterative refinement
2521 * O(N^3) complexity
2522
2523 INPUT PARAMETERS
2524 A - array[0..N-1,0..N-1], system matrix
2525 N - size of A
2526 B - array[0..N-1], right part
2527
2528 OUTPUT PARAMETERS
2529 Info - same as in RMatrixSolve
2530 Rep - same as in RMatrixSolve
2531 X - same as in RMatrixSolve
2532
2533 -- ALGLIB --
2534 Copyright 27.01.2010 by Bochkanov Sergey
2535 *************************************************************************/
cmatrixsolve(ae_matrix * a,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)2536 void cmatrixsolve(/* Complex */ ae_matrix* a,
2537 ae_int_t n,
2538 /* Complex */ ae_vector* b,
2539 ae_int_t* info,
2540 densesolverreport* rep,
2541 /* Complex */ ae_vector* x,
2542 ae_state *_state)
2543 {
2544 ae_frame _frame_block;
2545 ae_matrix bm;
2546 ae_matrix xm;
2547
2548 ae_frame_make(_state, &_frame_block);
2549 *info = 0;
2550 _densesolverreport_clear(rep);
2551 ae_vector_clear(x);
2552 ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
2553 ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
2554
2555 if( n<=0 )
2556 {
2557 *info = -1;
2558 ae_frame_leave(_state);
2559 return;
2560 }
2561 ae_matrix_set_length(&bm, n, 1, _state);
2562 ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
2563 cmatrixsolvem(a, n, &bm, 1, ae_true, info, rep, &xm, _state);
2564 ae_vector_set_length(x, n, _state);
2565 ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
2566 ae_frame_leave(_state);
2567 }
2568
2569
2570 /*************************************************************************
2571 Dense solver. Same as RMatrixLUSolveM(), but for complex matrices.
2572
2573 Algorithm features:
2574 * automatic detection of degenerate cases
2575 * O(M*N^2) complexity
2576 * condition number estimation
2577
2578 No iterative refinement is provided because exact form of original matrix
2579 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
2580 need iterative refinement.
2581
2582 INPUT PARAMETERS
2583 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
2584 P - array[0..N-1], pivots array, RMatrixLU result
2585 N - size of A
2586 B - array[0..N-1,0..M-1], right part
2587 M - right part size
2588
2589 OUTPUT PARAMETERS
2590 Info - same as in RMatrixSolve
2591 Rep - same as in RMatrixSolve
2592 X - same as in RMatrixSolve
2593
2594 -- ALGLIB --
2595 Copyright 27.01.2010 by Bochkanov Sergey
2596 *************************************************************************/
cmatrixlusolvem(ae_matrix * lua,ae_vector * p,ae_int_t n,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2597 void cmatrixlusolvem(/* Complex */ ae_matrix* lua,
2598 /* Integer */ ae_vector* p,
2599 ae_int_t n,
2600 /* Complex */ ae_matrix* b,
2601 ae_int_t m,
2602 ae_int_t* info,
2603 densesolverreport* rep,
2604 /* Complex */ ae_matrix* x,
2605 ae_state *_state)
2606 {
2607 ae_frame _frame_block;
2608 ae_matrix emptya;
2609 ae_int_t i;
2610 ae_int_t j;
2611 double scalea;
2612
2613 ae_frame_make(_state, &_frame_block);
2614 *info = 0;
2615 _densesolverreport_clear(rep);
2616 ae_matrix_clear(x);
2617 ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
2618
2619
2620 /*
2621 * prepare: check inputs, allocate space...
2622 */
2623 if( n<=0||m<=0 )
2624 {
2625 *info = -1;
2626 ae_frame_leave(_state);
2627 return;
2628 }
2629
2630 /*
2631 * 1. scale matrix, max(|U[i,j]|)
2632 * we assume that LU is in its normal form, i.e. |L[i,j]|<=1
2633 * 2. solve
2634 */
2635 scalea = 0;
2636 for(i=0; i<=n-1; i++)
2637 {
2638 for(j=i; j<=n-1; j++)
2639 {
2640 scalea = ae_maxreal(scalea, ae_c_abs(lua->ptr.pp_complex[i][j], _state), _state);
2641 }
2642 }
2643 if( ae_fp_eq(scalea,0) )
2644 {
2645 scalea = 1;
2646 }
2647 scalea = 1/scalea;
2648 densesolver_cmatrixlusolveinternal(lua, p, scalea, n, &emptya, ae_false, b, m, info, rep, x, _state);
2649 ae_frame_leave(_state);
2650 }
2651
2652
2653 /*************************************************************************
2654 Dense solver. Same as RMatrixLUSolve(), but for complex matrices.
2655
2656 Algorithm features:
2657 * automatic detection of degenerate cases
2658 * O(N^2) complexity
2659 * condition number estimation
2660
2661 No iterative refinement is provided because exact form of original matrix
2662 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
2663 need iterative refinement.
2664
2665 INPUT PARAMETERS
2666 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
2667 P - array[0..N-1], pivots array, CMatrixLU result
2668 N - size of A
2669 B - array[0..N-1], right part
2670
2671 OUTPUT PARAMETERS
2672 Info - same as in RMatrixSolve
2673 Rep - same as in RMatrixSolve
2674 X - same as in RMatrixSolve
2675
2676 -- ALGLIB --
2677 Copyright 27.01.2010 by Bochkanov Sergey
2678 *************************************************************************/
cmatrixlusolve(ae_matrix * lua,ae_vector * p,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)2679 void cmatrixlusolve(/* Complex */ ae_matrix* lua,
2680 /* Integer */ ae_vector* p,
2681 ae_int_t n,
2682 /* Complex */ ae_vector* b,
2683 ae_int_t* info,
2684 densesolverreport* rep,
2685 /* Complex */ ae_vector* x,
2686 ae_state *_state)
2687 {
2688 ae_frame _frame_block;
2689 ae_matrix bm;
2690 ae_matrix xm;
2691
2692 ae_frame_make(_state, &_frame_block);
2693 *info = 0;
2694 _densesolverreport_clear(rep);
2695 ae_vector_clear(x);
2696 ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
2697 ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
2698
2699 if( n<=0 )
2700 {
2701 *info = -1;
2702 ae_frame_leave(_state);
2703 return;
2704 }
2705 ae_matrix_set_length(&bm, n, 1, _state);
2706 ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
2707 cmatrixlusolvem(lua, p, n, &bm, 1, info, rep, &xm, _state);
2708 ae_vector_set_length(x, n, _state);
2709 ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
2710 ae_frame_leave(_state);
2711 }
2712
2713
2714 /*************************************************************************
2715 Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices.
2716
2717 Algorithm features:
2718 * automatic detection of degenerate cases
2719 * condition number estimation
2720 * iterative refinement
2721 * O(M*N^2) complexity
2722
2723 INPUT PARAMETERS
2724 A - array[0..N-1,0..N-1], system matrix
2725 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
2726 P - array[0..N-1], pivots array, CMatrixLU result
2727 N - size of A
2728 B - array[0..N-1,0..M-1], right part
2729 M - right part size
2730
2731 OUTPUT PARAMETERS
2732 Info - same as in RMatrixSolveM
2733 Rep - same as in RMatrixSolveM
2734 X - same as in RMatrixSolveM
2735
2736 -- ALGLIB --
2737 Copyright 27.01.2010 by Bochkanov Sergey
2738 *************************************************************************/
cmatrixmixedsolvem(ae_matrix * a,ae_matrix * lua,ae_vector * p,ae_int_t n,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2739 void cmatrixmixedsolvem(/* Complex */ ae_matrix* a,
2740 /* Complex */ ae_matrix* lua,
2741 /* Integer */ ae_vector* p,
2742 ae_int_t n,
2743 /* Complex */ ae_matrix* b,
2744 ae_int_t m,
2745 ae_int_t* info,
2746 densesolverreport* rep,
2747 /* Complex */ ae_matrix* x,
2748 ae_state *_state)
2749 {
2750 double scalea;
2751 ae_int_t i;
2752 ae_int_t j;
2753
2754 *info = 0;
2755 _densesolverreport_clear(rep);
2756 ae_matrix_clear(x);
2757
2758
2759 /*
2760 * prepare: check inputs, allocate space...
2761 */
2762 if( n<=0||m<=0 )
2763 {
2764 *info = -1;
2765 return;
2766 }
2767
2768 /*
2769 * 1. scale matrix, max(|A[i,j]|)
2770 * 2. factorize scaled matrix
2771 * 3. solve
2772 */
2773 scalea = 0;
2774 for(i=0; i<=n-1; i++)
2775 {
2776 for(j=0; j<=n-1; j++)
2777 {
2778 scalea = ae_maxreal(scalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
2779 }
2780 }
2781 if( ae_fp_eq(scalea,0) )
2782 {
2783 scalea = 1;
2784 }
2785 scalea = 1/scalea;
2786 densesolver_cmatrixlusolveinternal(lua, p, scalea, n, a, ae_true, b, m, info, rep, x, _state);
2787 }
2788
2789
2790 /*************************************************************************
2791 Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.
2792
2793 Algorithm features:
2794 * automatic detection of degenerate cases
2795 * condition number estimation
2796 * iterative refinement
2797 * O(N^2) complexity
2798
2799 INPUT PARAMETERS
2800 A - array[0..N-1,0..N-1], system matrix
2801 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
2802 P - array[0..N-1], pivots array, CMatrixLU result
2803 N - size of A
2804 B - array[0..N-1], right part
2805
2806 OUTPUT PARAMETERS
2807 Info - same as in RMatrixSolveM
2808 Rep - same as in RMatrixSolveM
2809 X - same as in RMatrixSolveM
2810
2811 -- ALGLIB --
2812 Copyright 27.01.2010 by Bochkanov Sergey
2813 *************************************************************************/
cmatrixmixedsolve(ae_matrix * a,ae_matrix * lua,ae_vector * p,ae_int_t n,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)2814 void cmatrixmixedsolve(/* Complex */ ae_matrix* a,
2815 /* Complex */ ae_matrix* lua,
2816 /* Integer */ ae_vector* p,
2817 ae_int_t n,
2818 /* Complex */ ae_vector* b,
2819 ae_int_t* info,
2820 densesolverreport* rep,
2821 /* Complex */ ae_vector* x,
2822 ae_state *_state)
2823 {
2824 ae_frame _frame_block;
2825 ae_matrix bm;
2826 ae_matrix xm;
2827
2828 ae_frame_make(_state, &_frame_block);
2829 *info = 0;
2830 _densesolverreport_clear(rep);
2831 ae_vector_clear(x);
2832 ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
2833 ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
2834
2835 if( n<=0 )
2836 {
2837 *info = -1;
2838 ae_frame_leave(_state);
2839 return;
2840 }
2841 ae_matrix_set_length(&bm, n, 1, _state);
2842 ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
2843 cmatrixmixedsolvem(a, lua, p, n, &bm, 1, info, rep, &xm, _state);
2844 ae_vector_set_length(x, n, _state);
2845 ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
2846 ae_frame_leave(_state);
2847 }
2848
2849
2850 /*************************************************************************
2851 Dense solver. Same as RMatrixSolveM(), but for symmetric positive definite
2852 matrices.
2853
2854 Algorithm features:
2855 * automatic detection of degenerate cases
2856 * condition number estimation
2857 * O(N^3+M*N^2) complexity
2858 * matrix is represented by its upper or lower triangle
2859
2860 No iterative refinement is provided because such partial representation of
2861 matrix does not allow efficient calculation of extra-precise matrix-vector
2862 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
2863 need iterative refinement.
2864
2865 INPUT PARAMETERS
2866 A - array[0..N-1,0..N-1], system matrix
2867 N - size of A
2868 IsUpper - what half of A is provided
2869 B - array[0..N-1,0..M-1], right part
2870 M - right part size
2871
2872 OUTPUT PARAMETERS
2873 Info - same as in RMatrixSolve.
2874 Returns -3 for non-SPD matrices.
2875 Rep - same as in RMatrixSolve
2876 X - same as in RMatrixSolve
2877
2878 -- ALGLIB --
2879 Copyright 27.01.2010 by Bochkanov Sergey
2880 *************************************************************************/
spdmatrixsolvem(ae_matrix * a,ae_int_t n,ae_bool isupper,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)2881 void spdmatrixsolvem(/* Real */ ae_matrix* a,
2882 ae_int_t n,
2883 ae_bool isupper,
2884 /* Real */ ae_matrix* b,
2885 ae_int_t m,
2886 ae_int_t* info,
2887 densesolverreport* rep,
2888 /* Real */ ae_matrix* x,
2889 ae_state *_state)
2890 {
2891 ae_frame _frame_block;
2892 ae_matrix da;
2893 double sqrtscalea;
2894 ae_int_t i;
2895 ae_int_t j;
2896 ae_int_t j1;
2897 ae_int_t j2;
2898
2899 ae_frame_make(_state, &_frame_block);
2900 *info = 0;
2901 _densesolverreport_clear(rep);
2902 ae_matrix_clear(x);
2903 ae_matrix_init(&da, 0, 0, DT_REAL, _state, ae_true);
2904
2905
2906 /*
2907 * prepare: check inputs, allocate space...
2908 */
2909 if( n<=0||m<=0 )
2910 {
2911 *info = -1;
2912 ae_frame_leave(_state);
2913 return;
2914 }
2915 ae_matrix_set_length(&da, n, n, _state);
2916
2917 /*
2918 * 1. scale matrix, max(|A[i,j]|)
2919 * 2. factorize scaled matrix
2920 * 3. solve
2921 */
2922 sqrtscalea = 0;
2923 for(i=0; i<=n-1; i++)
2924 {
2925 if( isupper )
2926 {
2927 j1 = i;
2928 j2 = n-1;
2929 }
2930 else
2931 {
2932 j1 = 0;
2933 j2 = i;
2934 }
2935 for(j=j1; j<=j2; j++)
2936 {
2937 sqrtscalea = ae_maxreal(sqrtscalea, ae_fabs(a->ptr.pp_double[i][j], _state), _state);
2938 }
2939 }
2940 if( ae_fp_eq(sqrtscalea,0) )
2941 {
2942 sqrtscalea = 1;
2943 }
2944 sqrtscalea = 1/sqrtscalea;
2945 sqrtscalea = ae_sqrt(sqrtscalea, _state);
2946 for(i=0; i<=n-1; i++)
2947 {
2948 if( isupper )
2949 {
2950 j1 = i;
2951 j2 = n-1;
2952 }
2953 else
2954 {
2955 j1 = 0;
2956 j2 = i;
2957 }
2958 ae_v_move(&da.ptr.pp_double[i][j1], 1, &a->ptr.pp_double[i][j1], 1, ae_v_len(j1,j2));
2959 }
2960 if( !spdmatrixcholesky(&da, n, isupper, _state) )
2961 {
2962 ae_matrix_set_length(x, n, m, _state);
2963 for(i=0; i<=n-1; i++)
2964 {
2965 for(j=0; j<=m-1; j++)
2966 {
2967 x->ptr.pp_double[i][j] = 0;
2968 }
2969 }
2970 rep->r1 = 0;
2971 rep->rinf = 0;
2972 *info = -3;
2973 ae_frame_leave(_state);
2974 return;
2975 }
2976 *info = 1;
2977 densesolver_spdmatrixcholeskysolveinternal(&da, sqrtscalea, n, isupper, a, ae_true, b, m, info, rep, x, _state);
2978 ae_frame_leave(_state);
2979 }
2980
2981
2982 /*************************************************************************
2983 Dense solver. Same as RMatrixSolve(), but for SPD matrices.
2984
2985 Algorithm features:
2986 * automatic detection of degenerate cases
2987 * condition number estimation
2988 * O(N^3) complexity
2989 * matrix is represented by its upper or lower triangle
2990
2991 No iterative refinement is provided because such partial representation of
2992 matrix does not allow efficient calculation of extra-precise matrix-vector
2993 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
2994 need iterative refinement.
2995
2996 INPUT PARAMETERS
2997 A - array[0..N-1,0..N-1], system matrix
2998 N - size of A
2999 IsUpper - what half of A is provided
3000 B - array[0..N-1], right part
3001
3002 OUTPUT PARAMETERS
3003 Info - same as in RMatrixSolve
3004 Returns -3 for non-SPD matrices.
3005 Rep - same as in RMatrixSolve
3006 X - same as in RMatrixSolve
3007
3008 -- ALGLIB --
3009 Copyright 27.01.2010 by Bochkanov Sergey
3010 *************************************************************************/
spdmatrixsolve(ae_matrix * a,ae_int_t n,ae_bool isupper,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)3011 void spdmatrixsolve(/* Real */ ae_matrix* a,
3012 ae_int_t n,
3013 ae_bool isupper,
3014 /* Real */ ae_vector* b,
3015 ae_int_t* info,
3016 densesolverreport* rep,
3017 /* Real */ ae_vector* x,
3018 ae_state *_state)
3019 {
3020 ae_frame _frame_block;
3021 ae_matrix bm;
3022 ae_matrix xm;
3023
3024 ae_frame_make(_state, &_frame_block);
3025 *info = 0;
3026 _densesolverreport_clear(rep);
3027 ae_vector_clear(x);
3028 ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
3029 ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
3030
3031 if( n<=0 )
3032 {
3033 *info = -1;
3034 ae_frame_leave(_state);
3035 return;
3036 }
3037 ae_matrix_set_length(&bm, n, 1, _state);
3038 ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3039 spdmatrixsolvem(a, n, isupper, &bm, 1, info, rep, &xm, _state);
3040 ae_vector_set_length(x, n, _state);
3041 ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
3042 ae_frame_leave(_state);
3043 }
3044
3045
3046 /*************************************************************************
3047 Dense solver. Same as RMatrixLUSolveM(), but for SPD matrices represented
3048 by their Cholesky decomposition.
3049
3050 Algorithm features:
3051 * automatic detection of degenerate cases
3052 * O(M*N^2) complexity
3053 * condition number estimation
3054 * matrix is represented by its upper or lower triangle
3055
3056 No iterative refinement is provided because such partial representation of
3057 matrix does not allow efficient calculation of extra-precise matrix-vector
3058 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3059 need iterative refinement.
3060
3061 INPUT PARAMETERS
3062 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
3063 SPDMatrixCholesky result
3064 N - size of CHA
3065 IsUpper - what half of CHA is provided
3066 B - array[0..N-1,0..M-1], right part
3067 M - right part size
3068
3069 OUTPUT PARAMETERS
3070 Info - same as in RMatrixSolve
3071 Rep - same as in RMatrixSolve
3072 X - same as in RMatrixSolve
3073
3074 -- ALGLIB --
3075 Copyright 27.01.2010 by Bochkanov Sergey
3076 *************************************************************************/
spdmatrixcholeskysolvem(ae_matrix * cha,ae_int_t n,ae_bool isupper,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)3077 void spdmatrixcholeskysolvem(/* Real */ ae_matrix* cha,
3078 ae_int_t n,
3079 ae_bool isupper,
3080 /* Real */ ae_matrix* b,
3081 ae_int_t m,
3082 ae_int_t* info,
3083 densesolverreport* rep,
3084 /* Real */ ae_matrix* x,
3085 ae_state *_state)
3086 {
3087 ae_frame _frame_block;
3088 ae_matrix emptya;
3089 double sqrtscalea;
3090 ae_int_t i;
3091 ae_int_t j;
3092 ae_int_t j1;
3093 ae_int_t j2;
3094
3095 ae_frame_make(_state, &_frame_block);
3096 *info = 0;
3097 _densesolverreport_clear(rep);
3098 ae_matrix_clear(x);
3099 ae_matrix_init(&emptya, 0, 0, DT_REAL, _state, ae_true);
3100
3101
3102 /*
3103 * prepare: check inputs, allocate space...
3104 */
3105 if( n<=0||m<=0 )
3106 {
3107 *info = -1;
3108 ae_frame_leave(_state);
3109 return;
3110 }
3111
3112 /*
3113 * 1. scale matrix, max(|U[i,j]|)
3114 * 2. factorize scaled matrix
3115 * 3. solve
3116 */
3117 sqrtscalea = 0;
3118 for(i=0; i<=n-1; i++)
3119 {
3120 if( isupper )
3121 {
3122 j1 = i;
3123 j2 = n-1;
3124 }
3125 else
3126 {
3127 j1 = 0;
3128 j2 = i;
3129 }
3130 for(j=j1; j<=j2; j++)
3131 {
3132 sqrtscalea = ae_maxreal(sqrtscalea, ae_fabs(cha->ptr.pp_double[i][j], _state), _state);
3133 }
3134 }
3135 if( ae_fp_eq(sqrtscalea,0) )
3136 {
3137 sqrtscalea = 1;
3138 }
3139 sqrtscalea = 1/sqrtscalea;
3140 densesolver_spdmatrixcholeskysolveinternal(cha, sqrtscalea, n, isupper, &emptya, ae_false, b, m, info, rep, x, _state);
3141 ae_frame_leave(_state);
3142 }
3143
3144
3145 /*************************************************************************
3146 Dense solver. Same as RMatrixLUSolve(), but for SPD matrices represented
3147 by their Cholesky decomposition.
3148
3149 Algorithm features:
3150 * automatic detection of degenerate cases
3151 * O(N^2) complexity
3152 * condition number estimation
3153 * matrix is represented by its upper or lower triangle
3154
3155 No iterative refinement is provided because such partial representation of
3156 matrix does not allow efficient calculation of extra-precise matrix-vector
3157 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3158 need iterative refinement.
3159
3160 INPUT PARAMETERS
3161 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
3162 SPDMatrixCholesky result
3163 N - size of A
3164 IsUpper - what half of CHA is provided
3165 B - array[0..N-1], right part
3166
3167 OUTPUT PARAMETERS
3168 Info - same as in RMatrixSolve
3169 Rep - same as in RMatrixSolve
3170 X - same as in RMatrixSolve
3171
3172 -- ALGLIB --
3173 Copyright 27.01.2010 by Bochkanov Sergey
3174 *************************************************************************/
spdmatrixcholeskysolve(ae_matrix * cha,ae_int_t n,ae_bool isupper,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)3175 void spdmatrixcholeskysolve(/* Real */ ae_matrix* cha,
3176 ae_int_t n,
3177 ae_bool isupper,
3178 /* Real */ ae_vector* b,
3179 ae_int_t* info,
3180 densesolverreport* rep,
3181 /* Real */ ae_vector* x,
3182 ae_state *_state)
3183 {
3184 ae_frame _frame_block;
3185 ae_matrix bm;
3186 ae_matrix xm;
3187
3188 ae_frame_make(_state, &_frame_block);
3189 *info = 0;
3190 _densesolverreport_clear(rep);
3191 ae_vector_clear(x);
3192 ae_matrix_init(&bm, 0, 0, DT_REAL, _state, ae_true);
3193 ae_matrix_init(&xm, 0, 0, DT_REAL, _state, ae_true);
3194
3195 if( n<=0 )
3196 {
3197 *info = -1;
3198 ae_frame_leave(_state);
3199 return;
3200 }
3201 ae_matrix_set_length(&bm, n, 1, _state);
3202 ae_v_move(&bm.ptr.pp_double[0][0], bm.stride, &b->ptr.p_double[0], 1, ae_v_len(0,n-1));
3203 spdmatrixcholeskysolvem(cha, n, isupper, &bm, 1, info, rep, &xm, _state);
3204 ae_vector_set_length(x, n, _state);
3205 ae_v_move(&x->ptr.p_double[0], 1, &xm.ptr.pp_double[0][0], xm.stride, ae_v_len(0,n-1));
3206 ae_frame_leave(_state);
3207 }
3208
3209
3210 /*************************************************************************
3211 Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
3212 matrices.
3213
3214 Algorithm features:
3215 * automatic detection of degenerate cases
3216 * condition number estimation
3217 * O(N^3+M*N^2) complexity
3218 * matrix is represented by its upper or lower triangle
3219
3220 No iterative refinement is provided because such partial representation of
3221 matrix does not allow efficient calculation of extra-precise matrix-vector
3222 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3223 need iterative refinement.
3224
3225 INPUT PARAMETERS
3226 A - array[0..N-1,0..N-1], system matrix
3227 N - size of A
3228 IsUpper - what half of A is provided
3229 B - array[0..N-1,0..M-1], right part
3230 M - right part size
3231
3232 OUTPUT PARAMETERS
3233 Info - same as in RMatrixSolve.
3234 Returns -3 for non-HPD matrices.
3235 Rep - same as in RMatrixSolve
3236 X - same as in RMatrixSolve
3237
3238 -- ALGLIB --
3239 Copyright 27.01.2010 by Bochkanov Sergey
3240 *************************************************************************/
hpdmatrixsolvem(ae_matrix * a,ae_int_t n,ae_bool isupper,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)3241 void hpdmatrixsolvem(/* Complex */ ae_matrix* a,
3242 ae_int_t n,
3243 ae_bool isupper,
3244 /* Complex */ ae_matrix* b,
3245 ae_int_t m,
3246 ae_int_t* info,
3247 densesolverreport* rep,
3248 /* Complex */ ae_matrix* x,
3249 ae_state *_state)
3250 {
3251 ae_frame _frame_block;
3252 ae_matrix da;
3253 double sqrtscalea;
3254 ae_int_t i;
3255 ae_int_t j;
3256 ae_int_t j1;
3257 ae_int_t j2;
3258
3259 ae_frame_make(_state, &_frame_block);
3260 *info = 0;
3261 _densesolverreport_clear(rep);
3262 ae_matrix_clear(x);
3263 ae_matrix_init(&da, 0, 0, DT_COMPLEX, _state, ae_true);
3264
3265
3266 /*
3267 * prepare: check inputs, allocate space...
3268 */
3269 if( n<=0||m<=0 )
3270 {
3271 *info = -1;
3272 ae_frame_leave(_state);
3273 return;
3274 }
3275 ae_matrix_set_length(&da, n, n, _state);
3276
3277 /*
3278 * 1. scale matrix, max(|A[i,j]|)
3279 * 2. factorize scaled matrix
3280 * 3. solve
3281 */
3282 sqrtscalea = 0;
3283 for(i=0; i<=n-1; i++)
3284 {
3285 if( isupper )
3286 {
3287 j1 = i;
3288 j2 = n-1;
3289 }
3290 else
3291 {
3292 j1 = 0;
3293 j2 = i;
3294 }
3295 for(j=j1; j<=j2; j++)
3296 {
3297 sqrtscalea = ae_maxreal(sqrtscalea, ae_c_abs(a->ptr.pp_complex[i][j], _state), _state);
3298 }
3299 }
3300 if( ae_fp_eq(sqrtscalea,0) )
3301 {
3302 sqrtscalea = 1;
3303 }
3304 sqrtscalea = 1/sqrtscalea;
3305 sqrtscalea = ae_sqrt(sqrtscalea, _state);
3306 for(i=0; i<=n-1; i++)
3307 {
3308 if( isupper )
3309 {
3310 j1 = i;
3311 j2 = n-1;
3312 }
3313 else
3314 {
3315 j1 = 0;
3316 j2 = i;
3317 }
3318 ae_v_cmove(&da.ptr.pp_complex[i][j1], 1, &a->ptr.pp_complex[i][j1], 1, "N", ae_v_len(j1,j2));
3319 }
3320 if( !hpdmatrixcholesky(&da, n, isupper, _state) )
3321 {
3322 ae_matrix_set_length(x, n, m, _state);
3323 for(i=0; i<=n-1; i++)
3324 {
3325 for(j=0; j<=m-1; j++)
3326 {
3327 x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
3328 }
3329 }
3330 rep->r1 = 0;
3331 rep->rinf = 0;
3332 *info = -3;
3333 ae_frame_leave(_state);
3334 return;
3335 }
3336 *info = 1;
3337 densesolver_hpdmatrixcholeskysolveinternal(&da, sqrtscalea, n, isupper, a, ae_true, b, m, info, rep, x, _state);
3338 ae_frame_leave(_state);
3339 }
3340
3341
3342 /*************************************************************************
3343 Dense solver. Same as RMatrixSolve(), but for Hermitian positive definite
3344 matrices.
3345
3346 Algorithm features:
3347 * automatic detection of degenerate cases
3348 * condition number estimation
3349 * O(N^3) complexity
3350 * matrix is represented by its upper or lower triangle
3351
3352 No iterative refinement is provided because such partial representation of
3353 matrix does not allow efficient calculation of extra-precise matrix-vector
3354 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3355 need iterative refinement.
3356
3357 INPUT PARAMETERS
3358 A - array[0..N-1,0..N-1], system matrix
3359 N - size of A
3360 IsUpper - what half of A is provided
3361 B - array[0..N-1], right part
3362
3363 OUTPUT PARAMETERS
3364 Info - same as in RMatrixSolve
3365 Returns -3 for non-HPD matrices.
3366 Rep - same as in RMatrixSolve
3367 X - same as in RMatrixSolve
3368
3369 -- ALGLIB --
3370 Copyright 27.01.2010 by Bochkanov Sergey
3371 *************************************************************************/
hpdmatrixsolve(ae_matrix * a,ae_int_t n,ae_bool isupper,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)3372 void hpdmatrixsolve(/* Complex */ ae_matrix* a,
3373 ae_int_t n,
3374 ae_bool isupper,
3375 /* Complex */ ae_vector* b,
3376 ae_int_t* info,
3377 densesolverreport* rep,
3378 /* Complex */ ae_vector* x,
3379 ae_state *_state)
3380 {
3381 ae_frame _frame_block;
3382 ae_matrix bm;
3383 ae_matrix xm;
3384
3385 ae_frame_make(_state, &_frame_block);
3386 *info = 0;
3387 _densesolverreport_clear(rep);
3388 ae_vector_clear(x);
3389 ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
3390 ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
3391
3392 if( n<=0 )
3393 {
3394 *info = -1;
3395 ae_frame_leave(_state);
3396 return;
3397 }
3398 ae_matrix_set_length(&bm, n, 1, _state);
3399 ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
3400 hpdmatrixsolvem(a, n, isupper, &bm, 1, info, rep, &xm, _state);
3401 ae_vector_set_length(x, n, _state);
3402 ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
3403 ae_frame_leave(_state);
3404 }
3405
3406
3407 /*************************************************************************
3408 Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices represented
3409 by their Cholesky decomposition.
3410
3411 Algorithm features:
3412 * automatic detection of degenerate cases
3413 * O(M*N^2) complexity
3414 * condition number estimation
3415 * matrix is represented by its upper or lower triangle
3416
3417 No iterative refinement is provided because such partial representation of
3418 matrix does not allow efficient calculation of extra-precise matrix-vector
3419 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3420 need iterative refinement.
3421
3422 INPUT PARAMETERS
3423 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
3424 HPDMatrixCholesky result
3425 N - size of CHA
3426 IsUpper - what half of CHA is provided
3427 B - array[0..N-1,0..M-1], right part
3428 M - right part size
3429
3430 OUTPUT PARAMETERS
3431 Info - same as in RMatrixSolve
3432 Rep - same as in RMatrixSolve
3433 X - same as in RMatrixSolve
3434
3435 -- ALGLIB --
3436 Copyright 27.01.2010 by Bochkanov Sergey
3437 *************************************************************************/
hpdmatrixcholeskysolvem(ae_matrix * cha,ae_int_t n,ae_bool isupper,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)3438 void hpdmatrixcholeskysolvem(/* Complex */ ae_matrix* cha,
3439 ae_int_t n,
3440 ae_bool isupper,
3441 /* Complex */ ae_matrix* b,
3442 ae_int_t m,
3443 ae_int_t* info,
3444 densesolverreport* rep,
3445 /* Complex */ ae_matrix* x,
3446 ae_state *_state)
3447 {
3448 ae_frame _frame_block;
3449 ae_matrix emptya;
3450 double sqrtscalea;
3451 ae_int_t i;
3452 ae_int_t j;
3453 ae_int_t j1;
3454 ae_int_t j2;
3455
3456 ae_frame_make(_state, &_frame_block);
3457 *info = 0;
3458 _densesolverreport_clear(rep);
3459 ae_matrix_clear(x);
3460 ae_matrix_init(&emptya, 0, 0, DT_COMPLEX, _state, ae_true);
3461
3462
3463 /*
3464 * prepare: check inputs, allocate space...
3465 */
3466 if( n<=0||m<=0 )
3467 {
3468 *info = -1;
3469 ae_frame_leave(_state);
3470 return;
3471 }
3472
3473 /*
3474 * 1. scale matrix, max(|U[i,j]|)
3475 * 2. factorize scaled matrix
3476 * 3. solve
3477 */
3478 sqrtscalea = 0;
3479 for(i=0; i<=n-1; i++)
3480 {
3481 if( isupper )
3482 {
3483 j1 = i;
3484 j2 = n-1;
3485 }
3486 else
3487 {
3488 j1 = 0;
3489 j2 = i;
3490 }
3491 for(j=j1; j<=j2; j++)
3492 {
3493 sqrtscalea = ae_maxreal(sqrtscalea, ae_c_abs(cha->ptr.pp_complex[i][j], _state), _state);
3494 }
3495 }
3496 if( ae_fp_eq(sqrtscalea,0) )
3497 {
3498 sqrtscalea = 1;
3499 }
3500 sqrtscalea = 1/sqrtscalea;
3501 densesolver_hpdmatrixcholeskysolveinternal(cha, sqrtscalea, n, isupper, &emptya, ae_false, b, m, info, rep, x, _state);
3502 ae_frame_leave(_state);
3503 }
3504
3505
3506 /*************************************************************************
3507 Dense solver. Same as RMatrixLUSolve(), but for HPD matrices represented
3508 by their Cholesky decomposition.
3509
3510 Algorithm features:
3511 * automatic detection of degenerate cases
3512 * O(N^2) complexity
3513 * condition number estimation
3514 * matrix is represented by its upper or lower triangle
3515
3516 No iterative refinement is provided because such partial representation of
3517 matrix does not allow efficient calculation of extra-precise matrix-vector
3518 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
3519 need iterative refinement.
3520
3521 INPUT PARAMETERS
3522 CHA - array[0..N-1,0..N-1], Cholesky decomposition,
3523 SPDMatrixCholesky result
3524 N - size of A
3525 IsUpper - what half of CHA is provided
3526 B - array[0..N-1], right part
3527
3528 OUTPUT PARAMETERS
3529 Info - same as in RMatrixSolve
3530 Rep - same as in RMatrixSolve
3531 X - same as in RMatrixSolve
3532
3533 -- ALGLIB --
3534 Copyright 27.01.2010 by Bochkanov Sergey
3535 *************************************************************************/
hpdmatrixcholeskysolve(ae_matrix * cha,ae_int_t n,ae_bool isupper,ae_vector * b,ae_int_t * info,densesolverreport * rep,ae_vector * x,ae_state * _state)3536 void hpdmatrixcholeskysolve(/* Complex */ ae_matrix* cha,
3537 ae_int_t n,
3538 ae_bool isupper,
3539 /* Complex */ ae_vector* b,
3540 ae_int_t* info,
3541 densesolverreport* rep,
3542 /* Complex */ ae_vector* x,
3543 ae_state *_state)
3544 {
3545 ae_frame _frame_block;
3546 ae_matrix bm;
3547 ae_matrix xm;
3548
3549 ae_frame_make(_state, &_frame_block);
3550 *info = 0;
3551 _densesolverreport_clear(rep);
3552 ae_vector_clear(x);
3553 ae_matrix_init(&bm, 0, 0, DT_COMPLEX, _state, ae_true);
3554 ae_matrix_init(&xm, 0, 0, DT_COMPLEX, _state, ae_true);
3555
3556 if( n<=0 )
3557 {
3558 *info = -1;
3559 ae_frame_leave(_state);
3560 return;
3561 }
3562 ae_matrix_set_length(&bm, n, 1, _state);
3563 ae_v_cmove(&bm.ptr.pp_complex[0][0], bm.stride, &b->ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
3564 hpdmatrixcholeskysolvem(cha, n, isupper, &bm, 1, info, rep, &xm, _state);
3565 ae_vector_set_length(x, n, _state);
3566 ae_v_cmove(&x->ptr.p_complex[0], 1, &xm.ptr.pp_complex[0][0], xm.stride, "N", ae_v_len(0,n-1));
3567 ae_frame_leave(_state);
3568 }
3569
3570
3571 /*************************************************************************
3572 Dense solver.
3573
3574 This subroutine finds solution of the linear system A*X=B with non-square,
3575 possibly degenerate A. System is solved in the least squares sense, and
3576 general least squares solution X = X0 + CX*y which minimizes |A*X-B| is
3577 returned. If A is non-degenerate, solution in the usual sense is returned
3578
3579 Algorithm features:
3580 * automatic detection of degenerate cases
3581 * iterative refinement
3582 * O(N^3) complexity
3583
3584 INPUT PARAMETERS
3585 A - array[0..NRows-1,0..NCols-1], system matrix
3586 NRows - vertical size of A
3587 NCols - horizontal size of A
3588 B - array[0..NCols-1], right part
3589 Threshold- a number in [0,1]. Singular values beyond Threshold are
3590 considered zero. Set it to 0.0, if you don't understand
3591 what it means, so the solver will choose good value on its
3592 own.
3593
3594 OUTPUT PARAMETERS
3595 Info - return code:
3596 * -4 SVD subroutine failed
3597 * -1 if NRows<=0 or NCols<=0 or Threshold<0 was passed
3598 * 1 if task is solved
3599 Rep - solver report, see below for more info
3600 X - array[0..N-1,0..M-1], it contains:
3601 * solution of A*X=B if A is non-singular (well-conditioned
3602 or ill-conditioned, but not very close to singular)
3603 * zeros, if A is singular or VERY close to singular
3604 (in this case Info=-3).
3605
3606 SOLVER REPORT
3607
3608 Subroutine sets following fields of the Rep structure:
3609 * R2 reciprocal of condition number: 1/cond(A), 2-norm.
3610 * N = NCols
3611 * K dim(Null(A))
3612 * CX array[0..N-1,0..K-1], kernel of A.
3613 Columns of CX store such vectors that A*CX[i]=0.
3614
3615 -- ALGLIB --
3616 Copyright 24.08.2009 by Bochkanov Sergey
3617 *************************************************************************/
rmatrixsolvels(ae_matrix * a,ae_int_t nrows,ae_int_t ncols,ae_vector * b,double threshold,ae_int_t * info,densesolverlsreport * rep,ae_vector * x,ae_state * _state)3618 void rmatrixsolvels(/* Real */ ae_matrix* a,
3619 ae_int_t nrows,
3620 ae_int_t ncols,
3621 /* Real */ ae_vector* b,
3622 double threshold,
3623 ae_int_t* info,
3624 densesolverlsreport* rep,
3625 /* Real */ ae_vector* x,
3626 ae_state *_state)
3627 {
3628 ae_frame _frame_block;
3629 ae_vector sv;
3630 ae_matrix u;
3631 ae_matrix vt;
3632 ae_vector rp;
3633 ae_vector utb;
3634 ae_vector sutb;
3635 ae_vector tmp;
3636 ae_vector ta;
3637 ae_vector tx;
3638 ae_vector buf;
3639 ae_vector w;
3640 ae_int_t i;
3641 ae_int_t j;
3642 ae_int_t nsv;
3643 ae_int_t kernelidx;
3644 double v;
3645 double verr;
3646 ae_bool svdfailed;
3647 ae_bool zeroa;
3648 ae_int_t rfs;
3649 ae_int_t nrfs;
3650 ae_bool terminatenexttime;
3651 ae_bool smallerr;
3652
3653 ae_frame_make(_state, &_frame_block);
3654 *info = 0;
3655 _densesolverlsreport_clear(rep);
3656 ae_vector_clear(x);
3657 ae_vector_init(&sv, 0, DT_REAL, _state, ae_true);
3658 ae_matrix_init(&u, 0, 0, DT_REAL, _state, ae_true);
3659 ae_matrix_init(&vt, 0, 0, DT_REAL, _state, ae_true);
3660 ae_vector_init(&rp, 0, DT_REAL, _state, ae_true);
3661 ae_vector_init(&utb, 0, DT_REAL, _state, ae_true);
3662 ae_vector_init(&sutb, 0, DT_REAL, _state, ae_true);
3663 ae_vector_init(&tmp, 0, DT_REAL, _state, ae_true);
3664 ae_vector_init(&ta, 0, DT_REAL, _state, ae_true);
3665 ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
3666 ae_vector_init(&buf, 0, DT_REAL, _state, ae_true);
3667 ae_vector_init(&w, 0, DT_REAL, _state, ae_true);
3668
3669 if( (nrows<=0||ncols<=0)||ae_fp_less(threshold,0) )
3670 {
3671 *info = -1;
3672 ae_frame_leave(_state);
3673 return;
3674 }
3675 if( ae_fp_eq(threshold,0) )
3676 {
3677 threshold = 1000*ae_machineepsilon;
3678 }
3679
3680 /*
3681 * Factorize A first
3682 */
3683 svdfailed = !rmatrixsvd(a, nrows, ncols, 1, 2, 2, &sv, &u, &vt, _state);
3684 zeroa = ae_fp_eq(sv.ptr.p_double[0],0);
3685 if( svdfailed||zeroa )
3686 {
3687 if( svdfailed )
3688 {
3689 *info = -4;
3690 }
3691 else
3692 {
3693 *info = 1;
3694 }
3695 ae_vector_set_length(x, ncols, _state);
3696 for(i=0; i<=ncols-1; i++)
3697 {
3698 x->ptr.p_double[i] = 0;
3699 }
3700 rep->n = ncols;
3701 rep->k = ncols;
3702 ae_matrix_set_length(&rep->cx, ncols, ncols, _state);
3703 for(i=0; i<=ncols-1; i++)
3704 {
3705 for(j=0; j<=ncols-1; j++)
3706 {
3707 if( i==j )
3708 {
3709 rep->cx.ptr.pp_double[i][j] = 1;
3710 }
3711 else
3712 {
3713 rep->cx.ptr.pp_double[i][j] = 0;
3714 }
3715 }
3716 }
3717 rep->r2 = 0;
3718 ae_frame_leave(_state);
3719 return;
3720 }
3721 nsv = ae_minint(ncols, nrows, _state);
3722 if( nsv==ncols )
3723 {
3724 rep->r2 = sv.ptr.p_double[nsv-1]/sv.ptr.p_double[0];
3725 }
3726 else
3727 {
3728 rep->r2 = 0;
3729 }
3730 rep->n = ncols;
3731 *info = 1;
3732
3733 /*
3734 * Iterative refinement of xc combined with solution:
3735 * 1. xc = 0
3736 * 2. calculate r = bc-A*xc using extra-precise dot product
3737 * 3. solve A*y = r
3738 * 4. update x:=x+r
3739 * 5. goto 2
3740 *
3741 * This cycle is executed until one of two things happens:
3742 * 1. maximum number of iterations reached
3743 * 2. last iteration decreased error to the lower limit
3744 */
3745 ae_vector_set_length(&utb, nsv, _state);
3746 ae_vector_set_length(&sutb, nsv, _state);
3747 ae_vector_set_length(x, ncols, _state);
3748 ae_vector_set_length(&tmp, ncols, _state);
3749 ae_vector_set_length(&ta, ncols+1, _state);
3750 ae_vector_set_length(&tx, ncols+1, _state);
3751 ae_vector_set_length(&buf, ncols+1, _state);
3752 for(i=0; i<=ncols-1; i++)
3753 {
3754 x->ptr.p_double[i] = 0;
3755 }
3756 kernelidx = nsv;
3757 for(i=0; i<=nsv-1; i++)
3758 {
3759 if( ae_fp_less_eq(sv.ptr.p_double[i],threshold*sv.ptr.p_double[0]) )
3760 {
3761 kernelidx = i;
3762 break;
3763 }
3764 }
3765 rep->k = ncols-kernelidx;
3766 nrfs = densesolver_densesolverrfsmaxv2(ncols, rep->r2, _state);
3767 terminatenexttime = ae_false;
3768 ae_vector_set_length(&rp, nrows, _state);
3769 for(rfs=0; rfs<=nrfs; rfs++)
3770 {
3771 if( terminatenexttime )
3772 {
3773 break;
3774 }
3775
3776 /*
3777 * calculate right part
3778 */
3779 if( rfs==0 )
3780 {
3781 ae_v_move(&rp.ptr.p_double[0], 1, &b->ptr.p_double[0], 1, ae_v_len(0,nrows-1));
3782 }
3783 else
3784 {
3785 smallerr = ae_true;
3786 for(i=0; i<=nrows-1; i++)
3787 {
3788 ae_v_move(&ta.ptr.p_double[0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,ncols-1));
3789 ta.ptr.p_double[ncols] = -1;
3790 ae_v_move(&tx.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,ncols-1));
3791 tx.ptr.p_double[ncols] = b->ptr.p_double[i];
3792 xdot(&ta, &tx, ncols+1, &buf, &v, &verr, _state);
3793 rp.ptr.p_double[i] = -v;
3794 smallerr = smallerr&&ae_fp_less(ae_fabs(v, _state),4*verr);
3795 }
3796 if( smallerr )
3797 {
3798 terminatenexttime = ae_true;
3799 }
3800 }
3801
3802 /*
3803 * solve A*dx = rp
3804 */
3805 for(i=0; i<=ncols-1; i++)
3806 {
3807 tmp.ptr.p_double[i] = 0;
3808 }
3809 for(i=0; i<=nsv-1; i++)
3810 {
3811 utb.ptr.p_double[i] = 0;
3812 }
3813 for(i=0; i<=nrows-1; i++)
3814 {
3815 v = rp.ptr.p_double[i];
3816 ae_v_addd(&utb.ptr.p_double[0], 1, &u.ptr.pp_double[i][0], 1, ae_v_len(0,nsv-1), v);
3817 }
3818 for(i=0; i<=nsv-1; i++)
3819 {
3820 if( i<kernelidx )
3821 {
3822 sutb.ptr.p_double[i] = utb.ptr.p_double[i]/sv.ptr.p_double[i];
3823 }
3824 else
3825 {
3826 sutb.ptr.p_double[i] = 0;
3827 }
3828 }
3829 for(i=0; i<=nsv-1; i++)
3830 {
3831 v = sutb.ptr.p_double[i];
3832 ae_v_addd(&tmp.ptr.p_double[0], 1, &vt.ptr.pp_double[i][0], 1, ae_v_len(0,ncols-1), v);
3833 }
3834
3835 /*
3836 * update x: x:=x+dx
3837 */
3838 ae_v_add(&x->ptr.p_double[0], 1, &tmp.ptr.p_double[0], 1, ae_v_len(0,ncols-1));
3839 }
3840
3841 /*
3842 * fill CX
3843 */
3844 if( rep->k>0 )
3845 {
3846 ae_matrix_set_length(&rep->cx, ncols, rep->k, _state);
3847 for(i=0; i<=rep->k-1; i++)
3848 {
3849 ae_v_move(&rep->cx.ptr.pp_double[0][i], rep->cx.stride, &vt.ptr.pp_double[kernelidx+i][0], 1, ae_v_len(0,ncols-1));
3850 }
3851 }
3852 ae_frame_leave(_state);
3853 }
3854
3855
3856 /*************************************************************************
3857 Internal LU solver
3858
3859 -- ALGLIB --
3860 Copyright 27.01.2010 by Bochkanov Sergey
3861 *************************************************************************/
densesolver_rmatrixlusolveinternal(ae_matrix * lua,ae_vector * p,double scalea,ae_int_t n,ae_matrix * a,ae_bool havea,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)3862 static void densesolver_rmatrixlusolveinternal(/* Real */ ae_matrix* lua,
3863 /* Integer */ ae_vector* p,
3864 double scalea,
3865 ae_int_t n,
3866 /* Real */ ae_matrix* a,
3867 ae_bool havea,
3868 /* Real */ ae_matrix* b,
3869 ae_int_t m,
3870 ae_int_t* info,
3871 densesolverreport* rep,
3872 /* Real */ ae_matrix* x,
3873 ae_state *_state)
3874 {
3875 ae_frame _frame_block;
3876 ae_int_t i;
3877 ae_int_t j;
3878 ae_int_t k;
3879 ae_int_t rfs;
3880 ae_int_t nrfs;
3881 ae_vector xc;
3882 ae_vector y;
3883 ae_vector bc;
3884 ae_vector xa;
3885 ae_vector xb;
3886 ae_vector tx;
3887 double v;
3888 double verr;
3889 double mxb;
3890 double scaleright;
3891 ae_bool smallerr;
3892 ae_bool terminatenexttime;
3893
3894 ae_frame_make(_state, &_frame_block);
3895 *info = 0;
3896 _densesolverreport_clear(rep);
3897 ae_matrix_clear(x);
3898 ae_vector_init(&xc, 0, DT_REAL, _state, ae_true);
3899 ae_vector_init(&y, 0, DT_REAL, _state, ae_true);
3900 ae_vector_init(&bc, 0, DT_REAL, _state, ae_true);
3901 ae_vector_init(&xa, 0, DT_REAL, _state, ae_true);
3902 ae_vector_init(&xb, 0, DT_REAL, _state, ae_true);
3903 ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
3904
3905 ae_assert(ae_fp_greater(scalea,0), "Assertion failed", _state);
3906
3907 /*
3908 * prepare: check inputs, allocate space...
3909 */
3910 if( n<=0||m<=0 )
3911 {
3912 *info = -1;
3913 ae_frame_leave(_state);
3914 return;
3915 }
3916 for(i=0; i<=n-1; i++)
3917 {
3918 if( p->ptr.p_int[i]>n-1||p->ptr.p_int[i]<i )
3919 {
3920 *info = -1;
3921 ae_frame_leave(_state);
3922 return;
3923 }
3924 }
3925 ae_matrix_set_length(x, n, m, _state);
3926 ae_vector_set_length(&y, n, _state);
3927 ae_vector_set_length(&xc, n, _state);
3928 ae_vector_set_length(&bc, n, _state);
3929 ae_vector_set_length(&tx, n+1, _state);
3930 ae_vector_set_length(&xa, n+1, _state);
3931 ae_vector_set_length(&xb, n+1, _state);
3932
3933 /*
3934 * estimate condition number, test for near singularity
3935 */
3936 rep->r1 = rmatrixlurcond1(lua, n, _state);
3937 rep->rinf = rmatrixlurcondinf(lua, n, _state);
3938 if( ae_fp_less(rep->r1,rcondthreshold(_state))||ae_fp_less(rep->rinf,rcondthreshold(_state)) )
3939 {
3940 for(i=0; i<=n-1; i++)
3941 {
3942 for(j=0; j<=m-1; j++)
3943 {
3944 x->ptr.pp_double[i][j] = 0;
3945 }
3946 }
3947 rep->r1 = 0;
3948 rep->rinf = 0;
3949 *info = -3;
3950 ae_frame_leave(_state);
3951 return;
3952 }
3953 *info = 1;
3954
3955 /*
3956 * solve
3957 */
3958 for(k=0; k<=m-1; k++)
3959 {
3960
3961 /*
3962 * copy B to contiguous storage
3963 */
3964 ae_v_move(&bc.ptr.p_double[0], 1, &b->ptr.pp_double[0][k], b->stride, ae_v_len(0,n-1));
3965
3966 /*
3967 * Scale right part:
3968 * * MX stores max(|Bi|)
3969 * * ScaleRight stores actual scaling applied to B when solving systems
3970 * it is chosen to make |scaleRight*b| close to 1.
3971 */
3972 mxb = 0;
3973 for(i=0; i<=n-1; i++)
3974 {
3975 mxb = ae_maxreal(mxb, ae_fabs(bc.ptr.p_double[i], _state), _state);
3976 }
3977 if( ae_fp_eq(mxb,0) )
3978 {
3979 mxb = 1;
3980 }
3981 scaleright = 1/mxb;
3982
3983 /*
3984 * First, non-iterative part of solution process.
3985 * We use separate code for this task because
3986 * XDot is quite slow and we want to save time.
3987 */
3988 ae_v_moved(&xc.ptr.p_double[0], 1, &bc.ptr.p_double[0], 1, ae_v_len(0,n-1), scaleright);
3989 densesolver_rbasiclusolve(lua, p, scalea, n, &xc, &tx, _state);
3990
3991 /*
3992 * Iterative refinement of xc:
3993 * * calculate r = bc-A*xc using extra-precise dot product
3994 * * solve A*y = r
3995 * * update x:=x+r
3996 *
3997 * This cycle is executed until one of two things happens:
3998 * 1. maximum number of iterations reached
3999 * 2. last iteration decreased error to the lower limit
4000 */
4001 if( havea )
4002 {
4003 nrfs = densesolver_densesolverrfsmax(n, rep->r1, rep->rinf, _state);
4004 terminatenexttime = ae_false;
4005 for(rfs=0; rfs<=nrfs-1; rfs++)
4006 {
4007 if( terminatenexttime )
4008 {
4009 break;
4010 }
4011
4012 /*
4013 * generate right part
4014 */
4015 smallerr = ae_true;
4016 ae_v_move(&xb.ptr.p_double[0], 1, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1));
4017 for(i=0; i<=n-1; i++)
4018 {
4019 ae_v_moved(&xa.ptr.p_double[0], 1, &a->ptr.pp_double[i][0], 1, ae_v_len(0,n-1), scalea);
4020 xa.ptr.p_double[n] = -1;
4021 xb.ptr.p_double[n] = scaleright*bc.ptr.p_double[i];
4022 xdot(&xa, &xb, n+1, &tx, &v, &verr, _state);
4023 y.ptr.p_double[i] = -v;
4024 smallerr = smallerr&&ae_fp_less(ae_fabs(v, _state),4*verr);
4025 }
4026 if( smallerr )
4027 {
4028 terminatenexttime = ae_true;
4029 }
4030
4031 /*
4032 * solve and update
4033 */
4034 densesolver_rbasiclusolve(lua, p, scalea, n, &y, &tx, _state);
4035 ae_v_add(&xc.ptr.p_double[0], 1, &y.ptr.p_double[0], 1, ae_v_len(0,n-1));
4036 }
4037 }
4038
4039 /*
4040 * Store xc.
4041 * Post-scale result.
4042 */
4043 v = scalea*mxb;
4044 ae_v_moved(&x->ptr.pp_double[0][k], x->stride, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1), v);
4045 }
4046 ae_frame_leave(_state);
4047 }
4048
4049
4050 /*************************************************************************
4051 Internal Cholesky solver
4052
4053 -- ALGLIB --
4054 Copyright 27.01.2010 by Bochkanov Sergey
4055 *************************************************************************/
densesolver_spdmatrixcholeskysolveinternal(ae_matrix * cha,double sqrtscalea,ae_int_t n,ae_bool isupper,ae_matrix * a,ae_bool havea,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)4056 static void densesolver_spdmatrixcholeskysolveinternal(/* Real */ ae_matrix* cha,
4057 double sqrtscalea,
4058 ae_int_t n,
4059 ae_bool isupper,
4060 /* Real */ ae_matrix* a,
4061 ae_bool havea,
4062 /* Real */ ae_matrix* b,
4063 ae_int_t m,
4064 ae_int_t* info,
4065 densesolverreport* rep,
4066 /* Real */ ae_matrix* x,
4067 ae_state *_state)
4068 {
4069 ae_frame _frame_block;
4070 ae_int_t i;
4071 ae_int_t j;
4072 ae_int_t k;
4073 ae_vector xc;
4074 ae_vector y;
4075 ae_vector bc;
4076 ae_vector xa;
4077 ae_vector xb;
4078 ae_vector tx;
4079 double v;
4080 double mxb;
4081 double scaleright;
4082
4083 ae_frame_make(_state, &_frame_block);
4084 *info = 0;
4085 _densesolverreport_clear(rep);
4086 ae_matrix_clear(x);
4087 ae_vector_init(&xc, 0, DT_REAL, _state, ae_true);
4088 ae_vector_init(&y, 0, DT_REAL, _state, ae_true);
4089 ae_vector_init(&bc, 0, DT_REAL, _state, ae_true);
4090 ae_vector_init(&xa, 0, DT_REAL, _state, ae_true);
4091 ae_vector_init(&xb, 0, DT_REAL, _state, ae_true);
4092 ae_vector_init(&tx, 0, DT_REAL, _state, ae_true);
4093
4094 ae_assert(ae_fp_greater(sqrtscalea,0), "Assertion failed", _state);
4095
4096 /*
4097 * prepare: check inputs, allocate space...
4098 */
4099 if( n<=0||m<=0 )
4100 {
4101 *info = -1;
4102 ae_frame_leave(_state);
4103 return;
4104 }
4105 ae_matrix_set_length(x, n, m, _state);
4106 ae_vector_set_length(&y, n, _state);
4107 ae_vector_set_length(&xc, n, _state);
4108 ae_vector_set_length(&bc, n, _state);
4109 ae_vector_set_length(&tx, n+1, _state);
4110 ae_vector_set_length(&xa, n+1, _state);
4111 ae_vector_set_length(&xb, n+1, _state);
4112
4113 /*
4114 * estimate condition number, test for near singularity
4115 */
4116 rep->r1 = spdmatrixcholeskyrcond(cha, n, isupper, _state);
4117 rep->rinf = rep->r1;
4118 if( ae_fp_less(rep->r1,rcondthreshold(_state)) )
4119 {
4120 for(i=0; i<=n-1; i++)
4121 {
4122 for(j=0; j<=m-1; j++)
4123 {
4124 x->ptr.pp_double[i][j] = 0;
4125 }
4126 }
4127 rep->r1 = 0;
4128 rep->rinf = 0;
4129 *info = -3;
4130 ae_frame_leave(_state);
4131 return;
4132 }
4133 *info = 1;
4134
4135 /*
4136 * solve
4137 */
4138 for(k=0; k<=m-1; k++)
4139 {
4140
4141 /*
4142 * copy B to contiguous storage
4143 */
4144 ae_v_move(&bc.ptr.p_double[0], 1, &b->ptr.pp_double[0][k], b->stride, ae_v_len(0,n-1));
4145
4146 /*
4147 * Scale right part:
4148 * * MX stores max(|Bi|)
4149 * * ScaleRight stores actual scaling applied to B when solving systems
4150 * it is chosen to make |scaleRight*b| close to 1.
4151 */
4152 mxb = 0;
4153 for(i=0; i<=n-1; i++)
4154 {
4155 mxb = ae_maxreal(mxb, ae_fabs(bc.ptr.p_double[i], _state), _state);
4156 }
4157 if( ae_fp_eq(mxb,0) )
4158 {
4159 mxb = 1;
4160 }
4161 scaleright = 1/mxb;
4162
4163 /*
4164 * First, non-iterative part of solution process.
4165 * We use separate code for this task because
4166 * XDot is quite slow and we want to save time.
4167 */
4168 ae_v_moved(&xc.ptr.p_double[0], 1, &bc.ptr.p_double[0], 1, ae_v_len(0,n-1), scaleright);
4169 densesolver_spdbasiccholeskysolve(cha, sqrtscalea, n, isupper, &xc, &tx, _state);
4170
4171 /*
4172 * Store xc.
4173 * Post-scale result.
4174 */
4175 v = ae_sqr(sqrtscalea, _state)*mxb;
4176 ae_v_moved(&x->ptr.pp_double[0][k], x->stride, &xc.ptr.p_double[0], 1, ae_v_len(0,n-1), v);
4177 }
4178 ae_frame_leave(_state);
4179 }
4180
4181
4182 /*************************************************************************
4183 Internal LU solver
4184
4185 -- ALGLIB --
4186 Copyright 27.01.2010 by Bochkanov Sergey
4187 *************************************************************************/
densesolver_cmatrixlusolveinternal(ae_matrix * lua,ae_vector * p,double scalea,ae_int_t n,ae_matrix * a,ae_bool havea,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)4188 static void densesolver_cmatrixlusolveinternal(/* Complex */ ae_matrix* lua,
4189 /* Integer */ ae_vector* p,
4190 double scalea,
4191 ae_int_t n,
4192 /* Complex */ ae_matrix* a,
4193 ae_bool havea,
4194 /* Complex */ ae_matrix* b,
4195 ae_int_t m,
4196 ae_int_t* info,
4197 densesolverreport* rep,
4198 /* Complex */ ae_matrix* x,
4199 ae_state *_state)
4200 {
4201 ae_frame _frame_block;
4202 ae_int_t i;
4203 ae_int_t j;
4204 ae_int_t k;
4205 ae_int_t rfs;
4206 ae_int_t nrfs;
4207 ae_vector xc;
4208 ae_vector y;
4209 ae_vector bc;
4210 ae_vector xa;
4211 ae_vector xb;
4212 ae_vector tx;
4213 ae_vector tmpbuf;
4214 ae_complex v;
4215 double verr;
4216 double mxb;
4217 double scaleright;
4218 ae_bool smallerr;
4219 ae_bool terminatenexttime;
4220
4221 ae_frame_make(_state, &_frame_block);
4222 *info = 0;
4223 _densesolverreport_clear(rep);
4224 ae_matrix_clear(x);
4225 ae_vector_init(&xc, 0, DT_COMPLEX, _state, ae_true);
4226 ae_vector_init(&y, 0, DT_COMPLEX, _state, ae_true);
4227 ae_vector_init(&bc, 0, DT_COMPLEX, _state, ae_true);
4228 ae_vector_init(&xa, 0, DT_COMPLEX, _state, ae_true);
4229 ae_vector_init(&xb, 0, DT_COMPLEX, _state, ae_true);
4230 ae_vector_init(&tx, 0, DT_COMPLEX, _state, ae_true);
4231 ae_vector_init(&tmpbuf, 0, DT_REAL, _state, ae_true);
4232
4233 ae_assert(ae_fp_greater(scalea,0), "Assertion failed", _state);
4234
4235 /*
4236 * prepare: check inputs, allocate space...
4237 */
4238 if( n<=0||m<=0 )
4239 {
4240 *info = -1;
4241 ae_frame_leave(_state);
4242 return;
4243 }
4244 for(i=0; i<=n-1; i++)
4245 {
4246 if( p->ptr.p_int[i]>n-1||p->ptr.p_int[i]<i )
4247 {
4248 *info = -1;
4249 ae_frame_leave(_state);
4250 return;
4251 }
4252 }
4253 ae_matrix_set_length(x, n, m, _state);
4254 ae_vector_set_length(&y, n, _state);
4255 ae_vector_set_length(&xc, n, _state);
4256 ae_vector_set_length(&bc, n, _state);
4257 ae_vector_set_length(&tx, n, _state);
4258 ae_vector_set_length(&xa, n+1, _state);
4259 ae_vector_set_length(&xb, n+1, _state);
4260 ae_vector_set_length(&tmpbuf, 2*n+2, _state);
4261
4262 /*
4263 * estimate condition number, test for near singularity
4264 */
4265 rep->r1 = cmatrixlurcond1(lua, n, _state);
4266 rep->rinf = cmatrixlurcondinf(lua, n, _state);
4267 if( ae_fp_less(rep->r1,rcondthreshold(_state))||ae_fp_less(rep->rinf,rcondthreshold(_state)) )
4268 {
4269 for(i=0; i<=n-1; i++)
4270 {
4271 for(j=0; j<=m-1; j++)
4272 {
4273 x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
4274 }
4275 }
4276 rep->r1 = 0;
4277 rep->rinf = 0;
4278 *info = -3;
4279 ae_frame_leave(_state);
4280 return;
4281 }
4282 *info = 1;
4283
4284 /*
4285 * solve
4286 */
4287 for(k=0; k<=m-1; k++)
4288 {
4289
4290 /*
4291 * copy B to contiguous storage
4292 */
4293 ae_v_cmove(&bc.ptr.p_complex[0], 1, &b->ptr.pp_complex[0][k], b->stride, "N", ae_v_len(0,n-1));
4294
4295 /*
4296 * Scale right part:
4297 * * MX stores max(|Bi|)
4298 * * ScaleRight stores actual scaling applied to B when solving systems
4299 * it is chosen to make |scaleRight*b| close to 1.
4300 */
4301 mxb = 0;
4302 for(i=0; i<=n-1; i++)
4303 {
4304 mxb = ae_maxreal(mxb, ae_c_abs(bc.ptr.p_complex[i], _state), _state);
4305 }
4306 if( ae_fp_eq(mxb,0) )
4307 {
4308 mxb = 1;
4309 }
4310 scaleright = 1/mxb;
4311
4312 /*
4313 * First, non-iterative part of solution process.
4314 * We use separate code for this task because
4315 * XDot is quite slow and we want to save time.
4316 */
4317 ae_v_cmoved(&xc.ptr.p_complex[0], 1, &bc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), scaleright);
4318 densesolver_cbasiclusolve(lua, p, scalea, n, &xc, &tx, _state);
4319
4320 /*
4321 * Iterative refinement of xc:
4322 * * calculate r = bc-A*xc using extra-precise dot product
4323 * * solve A*y = r
4324 * * update x:=x+r
4325 *
4326 * This cycle is executed until one of two things happens:
4327 * 1. maximum number of iterations reached
4328 * 2. last iteration decreased error to the lower limit
4329 */
4330 if( havea )
4331 {
4332 nrfs = densesolver_densesolverrfsmax(n, rep->r1, rep->rinf, _state);
4333 terminatenexttime = ae_false;
4334 for(rfs=0; rfs<=nrfs-1; rfs++)
4335 {
4336 if( terminatenexttime )
4337 {
4338 break;
4339 }
4340
4341 /*
4342 * generate right part
4343 */
4344 smallerr = ae_true;
4345 ae_v_cmove(&xb.ptr.p_complex[0], 1, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
4346 for(i=0; i<=n-1; i++)
4347 {
4348 ae_v_cmoved(&xa.ptr.p_complex[0], 1, &a->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,n-1), scalea);
4349 xa.ptr.p_complex[n] = ae_complex_from_d(-1);
4350 xb.ptr.p_complex[n] = ae_c_mul_d(bc.ptr.p_complex[i],scaleright);
4351 xcdot(&xa, &xb, n+1, &tmpbuf, &v, &verr, _state);
4352 y.ptr.p_complex[i] = ae_c_neg(v);
4353 smallerr = smallerr&&ae_fp_less(ae_c_abs(v, _state),4*verr);
4354 }
4355 if( smallerr )
4356 {
4357 terminatenexttime = ae_true;
4358 }
4359
4360 /*
4361 * solve and update
4362 */
4363 densesolver_cbasiclusolve(lua, p, scalea, n, &y, &tx, _state);
4364 ae_v_cadd(&xc.ptr.p_complex[0], 1, &y.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1));
4365 }
4366 }
4367
4368 /*
4369 * Store xc.
4370 * Post-scale result.
4371 */
4372 v = ae_complex_from_d(scalea*mxb);
4373 ae_v_cmovec(&x->ptr.pp_complex[0][k], x->stride, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), v);
4374 }
4375 ae_frame_leave(_state);
4376 }
4377
4378
4379 /*************************************************************************
4380 Internal Cholesky solver
4381
4382 -- ALGLIB --
4383 Copyright 27.01.2010 by Bochkanov Sergey
4384 *************************************************************************/
densesolver_hpdmatrixcholeskysolveinternal(ae_matrix * cha,double sqrtscalea,ae_int_t n,ae_bool isupper,ae_matrix * a,ae_bool havea,ae_matrix * b,ae_int_t m,ae_int_t * info,densesolverreport * rep,ae_matrix * x,ae_state * _state)4385 static void densesolver_hpdmatrixcholeskysolveinternal(/* Complex */ ae_matrix* cha,
4386 double sqrtscalea,
4387 ae_int_t n,
4388 ae_bool isupper,
4389 /* Complex */ ae_matrix* a,
4390 ae_bool havea,
4391 /* Complex */ ae_matrix* b,
4392 ae_int_t m,
4393 ae_int_t* info,
4394 densesolverreport* rep,
4395 /* Complex */ ae_matrix* x,
4396 ae_state *_state)
4397 {
4398 ae_frame _frame_block;
4399 ae_int_t i;
4400 ae_int_t j;
4401 ae_int_t k;
4402 ae_vector xc;
4403 ae_vector y;
4404 ae_vector bc;
4405 ae_vector xa;
4406 ae_vector xb;
4407 ae_vector tx;
4408 double v;
4409 double mxb;
4410 double scaleright;
4411
4412 ae_frame_make(_state, &_frame_block);
4413 *info = 0;
4414 _densesolverreport_clear(rep);
4415 ae_matrix_clear(x);
4416 ae_vector_init(&xc, 0, DT_COMPLEX, _state, ae_true);
4417 ae_vector_init(&y, 0, DT_COMPLEX, _state, ae_true);
4418 ae_vector_init(&bc, 0, DT_COMPLEX, _state, ae_true);
4419 ae_vector_init(&xa, 0, DT_COMPLEX, _state, ae_true);
4420 ae_vector_init(&xb, 0, DT_COMPLEX, _state, ae_true);
4421 ae_vector_init(&tx, 0, DT_COMPLEX, _state, ae_true);
4422
4423 ae_assert(ae_fp_greater(sqrtscalea,0), "Assertion failed", _state);
4424
4425 /*
4426 * prepare: check inputs, allocate space...
4427 */
4428 if( n<=0||m<=0 )
4429 {
4430 *info = -1;
4431 ae_frame_leave(_state);
4432 return;
4433 }
4434 ae_matrix_set_length(x, n, m, _state);
4435 ae_vector_set_length(&y, n, _state);
4436 ae_vector_set_length(&xc, n, _state);
4437 ae_vector_set_length(&bc, n, _state);
4438 ae_vector_set_length(&tx, n+1, _state);
4439 ae_vector_set_length(&xa, n+1, _state);
4440 ae_vector_set_length(&xb, n+1, _state);
4441
4442 /*
4443 * estimate condition number, test for near singularity
4444 */
4445 rep->r1 = hpdmatrixcholeskyrcond(cha, n, isupper, _state);
4446 rep->rinf = rep->r1;
4447 if( ae_fp_less(rep->r1,rcondthreshold(_state)) )
4448 {
4449 for(i=0; i<=n-1; i++)
4450 {
4451 for(j=0; j<=m-1; j++)
4452 {
4453 x->ptr.pp_complex[i][j] = ae_complex_from_d(0);
4454 }
4455 }
4456 rep->r1 = 0;
4457 rep->rinf = 0;
4458 *info = -3;
4459 ae_frame_leave(_state);
4460 return;
4461 }
4462 *info = 1;
4463
4464 /*
4465 * solve
4466 */
4467 for(k=0; k<=m-1; k++)
4468 {
4469
4470 /*
4471 * copy B to contiguous storage
4472 */
4473 ae_v_cmove(&bc.ptr.p_complex[0], 1, &b->ptr.pp_complex[0][k], b->stride, "N", ae_v_len(0,n-1));
4474
4475 /*
4476 * Scale right part:
4477 * * MX stores max(|Bi|)
4478 * * ScaleRight stores actual scaling applied to B when solving systems
4479 * it is chosen to make |scaleRight*b| close to 1.
4480 */
4481 mxb = 0;
4482 for(i=0; i<=n-1; i++)
4483 {
4484 mxb = ae_maxreal(mxb, ae_c_abs(bc.ptr.p_complex[i], _state), _state);
4485 }
4486 if( ae_fp_eq(mxb,0) )
4487 {
4488 mxb = 1;
4489 }
4490 scaleright = 1/mxb;
4491
4492 /*
4493 * First, non-iterative part of solution process.
4494 * We use separate code for this task because
4495 * XDot is quite slow and we want to save time.
4496 */
4497 ae_v_cmoved(&xc.ptr.p_complex[0], 1, &bc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), scaleright);
4498 densesolver_hpdbasiccholeskysolve(cha, sqrtscalea, n, isupper, &xc, &tx, _state);
4499
4500 /*
4501 * Store xc.
4502 * Post-scale result.
4503 */
4504 v = ae_sqr(sqrtscalea, _state)*mxb;
4505 ae_v_cmoved(&x->ptr.pp_complex[0][k], x->stride, &xc.ptr.p_complex[0], 1, "N", ae_v_len(0,n-1), v);
4506 }
4507 ae_frame_leave(_state);
4508 }
4509
4510
4511 /*************************************************************************
4512 Internal subroutine.
4513 Returns maximum count of RFS iterations as function of:
4514 1. machine epsilon
4515 2. task size.
4516 3. condition number
4517
4518 -- ALGLIB --
4519 Copyright 27.01.2010 by Bochkanov Sergey
4520 *************************************************************************/
densesolver_densesolverrfsmax(ae_int_t n,double r1,double rinf,ae_state * _state)4521 static ae_int_t densesolver_densesolverrfsmax(ae_int_t n,
4522 double r1,
4523 double rinf,
4524 ae_state *_state)
4525 {
4526 ae_int_t result;
4527
4528
4529 result = 5;
4530 return result;
4531 }
4532
4533
4534 /*************************************************************************
4535 Internal subroutine.
4536 Returns maximum count of RFS iterations as function of:
4537 1. machine epsilon
4538 2. task size.
4539 3. norm-2 condition number
4540
4541 -- ALGLIB --
4542 Copyright 27.01.2010 by Bochkanov Sergey
4543 *************************************************************************/
densesolver_densesolverrfsmaxv2(ae_int_t n,double r2,ae_state * _state)4544 static ae_int_t densesolver_densesolverrfsmaxv2(ae_int_t n,
4545 double r2,
4546 ae_state *_state)
4547 {
4548 ae_int_t result;
4549
4550
4551 result = densesolver_densesolverrfsmax(n, 0, 0, _state);
4552 return result;
4553 }
4554
4555
4556 /*************************************************************************
4557 Basic LU solver for ScaleA*PLU*x = y.
4558
4559 This subroutine assumes that:
4560 * L is well-scaled, and it is U which needs scaling by ScaleA.
4561 * A=PLU is well-conditioned, so no zero divisions or overflow may occur
4562
4563 -- ALGLIB --
4564 Copyright 27.01.2010 by Bochkanov Sergey
4565 *************************************************************************/
densesolver_rbasiclusolve(ae_matrix * lua,ae_vector * p,double scalea,ae_int_t n,ae_vector * xb,ae_vector * tmp,ae_state * _state)4566 static void densesolver_rbasiclusolve(/* Real */ ae_matrix* lua,
4567 /* Integer */ ae_vector* p,
4568 double scalea,
4569 ae_int_t n,
4570 /* Real */ ae_vector* xb,
4571 /* Real */ ae_vector* tmp,
4572 ae_state *_state)
4573 {
4574 ae_int_t i;
4575 double v;
4576
4577
4578 for(i=0; i<=n-1; i++)
4579 {
4580 if( p->ptr.p_int[i]!=i )
4581 {
4582 v = xb->ptr.p_double[i];
4583 xb->ptr.p_double[i] = xb->ptr.p_double[p->ptr.p_int[i]];
4584 xb->ptr.p_double[p->ptr.p_int[i]] = v;
4585 }
4586 }
4587 for(i=1; i<=n-1; i++)
4588 {
4589 v = ae_v_dotproduct(&lua->ptr.pp_double[i][0], 1, &xb->ptr.p_double[0], 1, ae_v_len(0,i-1));
4590 xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
4591 }
4592 xb->ptr.p_double[n-1] = xb->ptr.p_double[n-1]/(scalea*lua->ptr.pp_double[n-1][n-1]);
4593 for(i=n-2; i>=0; i--)
4594 {
4595 ae_v_moved(&tmp->ptr.p_double[i+1], 1, &lua->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), scalea);
4596 v = ae_v_dotproduct(&tmp->ptr.p_double[i+1], 1, &xb->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1));
4597 xb->ptr.p_double[i] = (xb->ptr.p_double[i]-v)/(scalea*lua->ptr.pp_double[i][i]);
4598 }
4599 }
4600
4601
4602 /*************************************************************************
4603 Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y.
4604
4605 This subroutine assumes that:
4606 * A*ScaleA is well scaled
4607 * A is well-conditioned, so no zero divisions or overflow may occur
4608
4609 -- ALGLIB --
4610 Copyright 27.01.2010 by Bochkanov Sergey
4611 *************************************************************************/
densesolver_spdbasiccholeskysolve(ae_matrix * cha,double sqrtscalea,ae_int_t n,ae_bool isupper,ae_vector * xb,ae_vector * tmp,ae_state * _state)4612 static void densesolver_spdbasiccholeskysolve(/* Real */ ae_matrix* cha,
4613 double sqrtscalea,
4614 ae_int_t n,
4615 ae_bool isupper,
4616 /* Real */ ae_vector* xb,
4617 /* Real */ ae_vector* tmp,
4618 ae_state *_state)
4619 {
4620 ae_int_t i;
4621 double v;
4622
4623
4624
4625 /*
4626 * A = L*L' or A=U'*U
4627 */
4628 if( isupper )
4629 {
4630
4631 /*
4632 * Solve U'*y=b first.
4633 */
4634 for(i=0; i<=n-1; i++)
4635 {
4636 xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
4637 if( i<n-1 )
4638 {
4639 v = xb->ptr.p_double[i];
4640 ae_v_moved(&tmp->ptr.p_double[i+1], 1, &cha->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), sqrtscalea);
4641 ae_v_subd(&xb->ptr.p_double[i+1], 1, &tmp->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1), v);
4642 }
4643 }
4644
4645 /*
4646 * Solve U*x=y then.
4647 */
4648 for(i=n-1; i>=0; i--)
4649 {
4650 if( i<n-1 )
4651 {
4652 ae_v_moved(&tmp->ptr.p_double[i+1], 1, &cha->ptr.pp_double[i][i+1], 1, ae_v_len(i+1,n-1), sqrtscalea);
4653 v = ae_v_dotproduct(&tmp->ptr.p_double[i+1], 1, &xb->ptr.p_double[i+1], 1, ae_v_len(i+1,n-1));
4654 xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
4655 }
4656 xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
4657 }
4658 }
4659 else
4660 {
4661
4662 /*
4663 * Solve L*y=b first
4664 */
4665 for(i=0; i<=n-1; i++)
4666 {
4667 if( i>0 )
4668 {
4669 ae_v_moved(&tmp->ptr.p_double[0], 1, &cha->ptr.pp_double[i][0], 1, ae_v_len(0,i-1), sqrtscalea);
4670 v = ae_v_dotproduct(&tmp->ptr.p_double[0], 1, &xb->ptr.p_double[0], 1, ae_v_len(0,i-1));
4671 xb->ptr.p_double[i] = xb->ptr.p_double[i]-v;
4672 }
4673 xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
4674 }
4675
4676 /*
4677 * Solve L'*x=y then.
4678 */
4679 for(i=n-1; i>=0; i--)
4680 {
4681 xb->ptr.p_double[i] = xb->ptr.p_double[i]/(sqrtscalea*cha->ptr.pp_double[i][i]);
4682 if( i>0 )
4683 {
4684 v = xb->ptr.p_double[i];
4685 ae_v_moved(&tmp->ptr.p_double[0], 1, &cha->ptr.pp_double[i][0], 1, ae_v_len(0,i-1), sqrtscalea);
4686 ae_v_subd(&xb->ptr.p_double[0], 1, &tmp->ptr.p_double[0], 1, ae_v_len(0,i-1), v);
4687 }
4688 }
4689 }
4690 }
4691
4692
4693 /*************************************************************************
4694 Basic LU solver for ScaleA*PLU*x = y.
4695
4696 This subroutine assumes that:
4697 * L is well-scaled, and it is U which needs scaling by ScaleA.
4698 * A=PLU is well-conditioned, so no zero divisions or overflow may occur
4699
4700 -- ALGLIB --
4701 Copyright 27.01.2010 by Bochkanov Sergey
4702 *************************************************************************/
densesolver_cbasiclusolve(ae_matrix * lua,ae_vector * p,double scalea,ae_int_t n,ae_vector * xb,ae_vector * tmp,ae_state * _state)4703 static void densesolver_cbasiclusolve(/* Complex */ ae_matrix* lua,
4704 /* Integer */ ae_vector* p,
4705 double scalea,
4706 ae_int_t n,
4707 /* Complex */ ae_vector* xb,
4708 /* Complex */ ae_vector* tmp,
4709 ae_state *_state)
4710 {
4711 ae_int_t i;
4712 ae_complex v;
4713
4714
4715 for(i=0; i<=n-1; i++)
4716 {
4717 if( p->ptr.p_int[i]!=i )
4718 {
4719 v = xb->ptr.p_complex[i];
4720 xb->ptr.p_complex[i] = xb->ptr.p_complex[p->ptr.p_int[i]];
4721 xb->ptr.p_complex[p->ptr.p_int[i]] = v;
4722 }
4723 }
4724 for(i=1; i<=n-1; i++)
4725 {
4726 v = ae_v_cdotproduct(&lua->ptr.pp_complex[i][0], 1, "N", &xb->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1));
4727 xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
4728 }
4729 xb->ptr.p_complex[n-1] = ae_c_div(xb->ptr.p_complex[n-1],ae_c_mul_d(lua->ptr.pp_complex[n-1][n-1],scalea));
4730 for(i=n-2; i>=0; i--)
4731 {
4732 ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &lua->ptr.pp_complex[i][i+1], 1, "N", ae_v_len(i+1,n-1), scalea);
4733 v = ae_v_cdotproduct(&tmp->ptr.p_complex[i+1], 1, "N", &xb->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1));
4734 xb->ptr.p_complex[i] = ae_c_div(ae_c_sub(xb->ptr.p_complex[i],v),ae_c_mul_d(lua->ptr.pp_complex[i][i],scalea));
4735 }
4736 }
4737
4738
4739 /*************************************************************************
4740 Basic Cholesky solver for ScaleA*Cholesky(A)'*x = y.
4741
4742 This subroutine assumes that:
4743 * A*ScaleA is well scaled
4744 * A is well-conditioned, so no zero divisions or overflow may occur
4745
4746 -- ALGLIB --
4747 Copyright 27.01.2010 by Bochkanov Sergey
4748 *************************************************************************/
densesolver_hpdbasiccholeskysolve(ae_matrix * cha,double sqrtscalea,ae_int_t n,ae_bool isupper,ae_vector * xb,ae_vector * tmp,ae_state * _state)4749 static void densesolver_hpdbasiccholeskysolve(/* Complex */ ae_matrix* cha,
4750 double sqrtscalea,
4751 ae_int_t n,
4752 ae_bool isupper,
4753 /* Complex */ ae_vector* xb,
4754 /* Complex */ ae_vector* tmp,
4755 ae_state *_state)
4756 {
4757 ae_int_t i;
4758 ae_complex v;
4759
4760
4761
4762 /*
4763 * A = L*L' or A=U'*U
4764 */
4765 if( isupper )
4766 {
4767
4768 /*
4769 * Solve U'*y=b first.
4770 */
4771 for(i=0; i<=n-1; i++)
4772 {
4773 xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(ae_c_conj(cha->ptr.pp_complex[i][i], _state),sqrtscalea));
4774 if( i<n-1 )
4775 {
4776 v = xb->ptr.p_complex[i];
4777 ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &cha->ptr.pp_complex[i][i+1], 1, "Conj", ae_v_len(i+1,n-1), sqrtscalea);
4778 ae_v_csubc(&xb->ptr.p_complex[i+1], 1, &tmp->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1), v);
4779 }
4780 }
4781
4782 /*
4783 * Solve U*x=y then.
4784 */
4785 for(i=n-1; i>=0; i--)
4786 {
4787 if( i<n-1 )
4788 {
4789 ae_v_cmoved(&tmp->ptr.p_complex[i+1], 1, &cha->ptr.pp_complex[i][i+1], 1, "N", ae_v_len(i+1,n-1), sqrtscalea);
4790 v = ae_v_cdotproduct(&tmp->ptr.p_complex[i+1], 1, "N", &xb->ptr.p_complex[i+1], 1, "N", ae_v_len(i+1,n-1));
4791 xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
4792 }
4793 xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(cha->ptr.pp_complex[i][i],sqrtscalea));
4794 }
4795 }
4796 else
4797 {
4798
4799 /*
4800 * Solve L*y=b first
4801 */
4802 for(i=0; i<=n-1; i++)
4803 {
4804 if( i>0 )
4805 {
4806 ae_v_cmoved(&tmp->ptr.p_complex[0], 1, &cha->ptr.pp_complex[i][0], 1, "N", ae_v_len(0,i-1), sqrtscalea);
4807 v = ae_v_cdotproduct(&tmp->ptr.p_complex[0], 1, "N", &xb->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1));
4808 xb->ptr.p_complex[i] = ae_c_sub(xb->ptr.p_complex[i],v);
4809 }
4810 xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(cha->ptr.pp_complex[i][i],sqrtscalea));
4811 }
4812
4813 /*
4814 * Solve L'*x=y then.
4815 */
4816 for(i=n-1; i>=0; i--)
4817 {
4818 xb->ptr.p_complex[i] = ae_c_div(xb->ptr.p_complex[i],ae_c_mul_d(ae_c_conj(cha->ptr.pp_complex[i][i], _state),sqrtscalea));
4819 if( i>0 )
4820 {
4821 v = xb->ptr.p_complex[i];
4822 ae_v_cmoved(&tmp->ptr.p_complex[0], 1, &cha->ptr.pp_complex[i][0], 1, "Conj", ae_v_len(0,i-1), sqrtscalea);
4823 ae_v_csubc(&xb->ptr.p_complex[0], 1, &tmp->ptr.p_complex[0], 1, "N", ae_v_len(0,i-1), v);
4824 }
4825 }
4826 }
4827 }
4828
4829
_densesolverreport_init(densesolverreport * p,ae_state * _state,ae_bool make_automatic)4830 ae_bool _densesolverreport_init(densesolverreport* p, ae_state *_state, ae_bool make_automatic)
4831 {
4832 return ae_true;
4833 }
4834
4835
_densesolverreport_init_copy(densesolverreport * dst,densesolverreport * src,ae_state * _state,ae_bool make_automatic)4836 ae_bool _densesolverreport_init_copy(densesolverreport* dst, densesolverreport* src, ae_state *_state, ae_bool make_automatic)
4837 {
4838 dst->r1 = src->r1;
4839 dst->rinf = src->rinf;
4840 return ae_true;
4841 }
4842
4843
_densesolverreport_clear(densesolverreport * p)4844 void _densesolverreport_clear(densesolverreport* p)
4845 {
4846 }
4847
4848
_densesolverlsreport_init(densesolverlsreport * p,ae_state * _state,ae_bool make_automatic)4849 ae_bool _densesolverlsreport_init(densesolverlsreport* p, ae_state *_state, ae_bool make_automatic)
4850 {
4851 if( !ae_matrix_init(&p->cx, 0, 0, DT_REAL, _state, make_automatic) )
4852 return ae_false;
4853 return ae_true;
4854 }
4855
4856
_densesolverlsreport_init_copy(densesolverlsreport * dst,densesolverlsreport * src,ae_state * _state,ae_bool make_automatic)4857 ae_bool _densesolverlsreport_init_copy(densesolverlsreport* dst, densesolverlsreport* src, ae_state *_state, ae_bool make_automatic)
4858 {
4859 dst->r2 = src->r2;
4860 if( !ae_matrix_init_copy(&dst->cx, &src->cx, _state, make_automatic) )
4861 return ae_false;
4862 dst->n = src->n;
4863 dst->k = src->k;
4864 return ae_true;
4865 }
4866
4867
_densesolverlsreport_clear(densesolverlsreport * p)4868 void _densesolverlsreport_clear(densesolverlsreport* p)
4869 {
4870 ae_matrix_clear(&p->cx);
4871 }
4872
4873
4874
4875
4876 /*************************************************************************
4877 LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER
4878
4879 DESCRIPTION:
4880 This algorithm solves system of nonlinear equations
4881 F[0](x[0], ..., x[n-1]) = 0
4882 F[1](x[0], ..., x[n-1]) = 0
4883 ...
4884 F[M-1](x[0], ..., x[n-1]) = 0
4885 with M/N do not necessarily coincide. Algorithm converges quadratically
4886 under following conditions:
4887 * the solution set XS is nonempty
4888 * for some xs in XS there exist such neighbourhood N(xs) that:
4889 * vector function F(x) and its Jacobian J(x) are continuously
4890 differentiable on N
4891 * ||F(x)|| provides local error bound on N, i.e. there exists such
4892 c1, that ||F(x)||>c1*distance(x,XS)
4893 Note that these conditions are much more weaker than usual non-singularity
4894 conditions. For example, algorithm will converge for any affine function
4895 F (whether its Jacobian singular or not).
4896
4897
4898 REQUIREMENTS:
4899 Algorithm will request following information during its operation:
4900 * function vector F[] and Jacobian matrix at given point X
4901 * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X
4902
4903
4904 USAGE:
4905 1. User initializes algorithm state with NLEQCreateLM() call
4906 2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and
4907 other functions
4908 3. User calls NLEQSolve() function which takes algorithm state and
4909 pointers (delegates, etc.) to callback functions which calculate merit
4910 function value and Jacobian.
4911 4. User calls NLEQResults() to get solution
4912 5. Optionally, user may call NLEQRestartFrom() to solve another problem
4913 with same parameters (N/M) but another starting point and/or another
4914 function vector. NLEQRestartFrom() allows to reuse already initialized
4915 structure.
4916
4917
4918 INPUT PARAMETERS:
4919 N - space dimension, N>1:
4920 * if provided, only leading N elements of X are used
4921 * if not provided, determined automatically from size of X
4922 M - system size
4923 X - starting point
4924
4925
4926 OUTPUT PARAMETERS:
4927 State - structure which stores algorithm state
4928
4929
4930 NOTES:
4931 1. you may tune stopping conditions with NLEQSetCond() function
4932 2. if target function contains exp() or other fast growing functions, and
4933 optimization algorithm makes too large steps which leads to overflow,
4934 use NLEQSetStpMax() function to bound algorithm's steps.
4935 3. this algorithm is a slightly modified implementation of the method
4936 described in 'Levenberg-Marquardt method for constrained nonlinear
4937 equations with strong local convergence properties' by Christian Kanzow
4938 Nobuo Yamashita and Masao Fukushima and further developed in 'On the
4939 convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and
4940 Ya-Xiang Yuan.
4941
4942
4943 -- ALGLIB --
4944 Copyright 20.08.2009 by Bochkanov Sergey
4945 *************************************************************************/
nleqcreatelm(ae_int_t n,ae_int_t m,ae_vector * x,nleqstate * state,ae_state * _state)4946 void nleqcreatelm(ae_int_t n,
4947 ae_int_t m,
4948 /* Real */ ae_vector* x,
4949 nleqstate* state,
4950 ae_state *_state)
4951 {
4952
4953 _nleqstate_clear(state);
4954
4955 ae_assert(n>=1, "NLEQCreateLM: N<1!", _state);
4956 ae_assert(m>=1, "NLEQCreateLM: M<1!", _state);
4957 ae_assert(x->cnt>=n, "NLEQCreateLM: Length(X)<N!", _state);
4958 ae_assert(isfinitevector(x, n, _state), "NLEQCreateLM: X contains infinite or NaN values!", _state);
4959
4960 /*
4961 * Initialize
4962 */
4963 state->n = n;
4964 state->m = m;
4965 nleqsetcond(state, 0, 0, _state);
4966 nleqsetxrep(state, ae_false, _state);
4967 nleqsetstpmax(state, 0, _state);
4968 ae_vector_set_length(&state->x, n, _state);
4969 ae_vector_set_length(&state->xbase, n, _state);
4970 ae_matrix_set_length(&state->j, m, n, _state);
4971 ae_vector_set_length(&state->fi, m, _state);
4972 ae_vector_set_length(&state->rightpart, n, _state);
4973 ae_vector_set_length(&state->candstep, n, _state);
4974 nleqrestartfrom(state, x, _state);
4975 }
4976
4977
4978 /*************************************************************************
4979 This function sets stopping conditions for the nonlinear solver
4980
4981 INPUT PARAMETERS:
4982 State - structure which stores algorithm state
4983 EpsF - >=0
4984 The subroutine finishes its work if on k+1-th iteration
4985 the condition ||F||<=EpsF is satisfied
4986 MaxIts - maximum number of iterations. If MaxIts=0, the number of
4987 iterations is unlimited.
4988
4989 Passing EpsF=0 and MaxIts=0 simultaneously will lead to automatic
4990 stopping criterion selection (small EpsF).
4991
4992 NOTES:
4993
4994 -- ALGLIB --
4995 Copyright 20.08.2010 by Bochkanov Sergey
4996 *************************************************************************/
nleqsetcond(nleqstate * state,double epsf,ae_int_t maxits,ae_state * _state)4997 void nleqsetcond(nleqstate* state,
4998 double epsf,
4999 ae_int_t maxits,
5000 ae_state *_state)
5001 {
5002
5003
5004 ae_assert(ae_isfinite(epsf, _state), "NLEQSetCond: EpsF is not finite number!", _state);
5005 ae_assert(ae_fp_greater_eq(epsf,0), "NLEQSetCond: negative EpsF!", _state);
5006 ae_assert(maxits>=0, "NLEQSetCond: negative MaxIts!", _state);
5007 if( ae_fp_eq(epsf,0)&&maxits==0 )
5008 {
5009 epsf = 1.0E-6;
5010 }
5011 state->epsf = epsf;
5012 state->maxits = maxits;
5013 }
5014
5015
5016 /*************************************************************************
5017 This function turns on/off reporting.
5018
5019 INPUT PARAMETERS:
5020 State - structure which stores algorithm state
5021 NeedXRep- whether iteration reports are needed or not
5022
5023 If NeedXRep is True, algorithm will call rep() callback function if it is
5024 provided to NLEQSolve().
5025
5026 -- ALGLIB --
5027 Copyright 20.08.2010 by Bochkanov Sergey
5028 *************************************************************************/
nleqsetxrep(nleqstate * state,ae_bool needxrep,ae_state * _state)5029 void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state)
5030 {
5031
5032
5033 state->xrep = needxrep;
5034 }
5035
5036
5037 /*************************************************************************
5038 This function sets maximum step length
5039
5040 INPUT PARAMETERS:
5041 State - structure which stores algorithm state
5042 StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't
5043 want to limit step length.
5044
5045 Use this subroutine when target function contains exp() or other fast
5046 growing functions, and algorithm makes too large steps which lead to
5047 overflow. This function allows us to reject steps that are too large (and
5048 therefore expose us to the possible overflow) without actually calculating
5049 function value at the x+stp*d.
5050
5051 -- ALGLIB --
5052 Copyright 20.08.2010 by Bochkanov Sergey
5053 *************************************************************************/
nleqsetstpmax(nleqstate * state,double stpmax,ae_state * _state)5054 void nleqsetstpmax(nleqstate* state, double stpmax, ae_state *_state)
5055 {
5056
5057
5058 ae_assert(ae_isfinite(stpmax, _state), "NLEQSetStpMax: StpMax is not finite!", _state);
5059 ae_assert(ae_fp_greater_eq(stpmax,0), "NLEQSetStpMax: StpMax<0!", _state);
5060 state->stpmax = stpmax;
5061 }
5062
5063
5064 /*************************************************************************
5065
5066 -- ALGLIB --
5067 Copyright 20.03.2009 by Bochkanov Sergey
5068 *************************************************************************/
nleqiteration(nleqstate * state,ae_state * _state)5069 ae_bool nleqiteration(nleqstate* state, ae_state *_state)
5070 {
5071 ae_int_t n;
5072 ae_int_t m;
5073 ae_int_t i;
5074 double lambdaup;
5075 double lambdadown;
5076 double lambdav;
5077 double rho;
5078 double mu;
5079 double stepnorm;
5080 ae_bool b;
5081 ae_bool result;
5082
5083
5084
5085 /*
5086 * Reverse communication preparations
5087 * I know it looks ugly, but it works the same way
5088 * anywhere from C++ to Python.
5089 *
5090 * This code initializes locals by:
5091 * * random values determined during code
5092 * generation - on first subroutine call
5093 * * values from previous call - on subsequent calls
5094 */
5095 if( state->rstate.stage>=0 )
5096 {
5097 n = state->rstate.ia.ptr.p_int[0];
5098 m = state->rstate.ia.ptr.p_int[1];
5099 i = state->rstate.ia.ptr.p_int[2];
5100 b = state->rstate.ba.ptr.p_bool[0];
5101 lambdaup = state->rstate.ra.ptr.p_double[0];
5102 lambdadown = state->rstate.ra.ptr.p_double[1];
5103 lambdav = state->rstate.ra.ptr.p_double[2];
5104 rho = state->rstate.ra.ptr.p_double[3];
5105 mu = state->rstate.ra.ptr.p_double[4];
5106 stepnorm = state->rstate.ra.ptr.p_double[5];
5107 }
5108 else
5109 {
5110 n = -983;
5111 m = -989;
5112 i = -834;
5113 b = ae_false;
5114 lambdaup = -287;
5115 lambdadown = 364;
5116 lambdav = 214;
5117 rho = -338;
5118 mu = -686;
5119 stepnorm = 912;
5120 }
5121 if( state->rstate.stage==0 )
5122 {
5123 goto lbl_0;
5124 }
5125 if( state->rstate.stage==1 )
5126 {
5127 goto lbl_1;
5128 }
5129 if( state->rstate.stage==2 )
5130 {
5131 goto lbl_2;
5132 }
5133 if( state->rstate.stage==3 )
5134 {
5135 goto lbl_3;
5136 }
5137 if( state->rstate.stage==4 )
5138 {
5139 goto lbl_4;
5140 }
5141
5142 /*
5143 * Routine body
5144 */
5145
5146 /*
5147 * Prepare
5148 */
5149 n = state->n;
5150 m = state->m;
5151 state->repterminationtype = 0;
5152 state->repiterationscount = 0;
5153 state->repnfunc = 0;
5154 state->repnjac = 0;
5155
5156 /*
5157 * Calculate F/G, initialize algorithm
5158 */
5159 nleq_clearrequestfields(state, _state);
5160 state->needf = ae_true;
5161 state->rstate.stage = 0;
5162 goto lbl_rcomm;
5163 lbl_0:
5164 state->needf = ae_false;
5165 state->repnfunc = state->repnfunc+1;
5166 ae_v_move(&state->xbase.ptr.p_double[0], 1, &state->x.ptr.p_double[0], 1, ae_v_len(0,n-1));
5167 state->fbase = state->f;
5168 state->fprev = ae_maxrealnumber;
5169 if( !state->xrep )
5170 {
5171 goto lbl_5;
5172 }
5173
5174 /*
5175 * progress report
5176 */
5177 nleq_clearrequestfields(state, _state);
5178 state->xupdated = ae_true;
5179 state->rstate.stage = 1;
5180 goto lbl_rcomm;
5181 lbl_1:
5182 state->xupdated = ae_false;
5183 lbl_5:
5184 if( ae_fp_less_eq(state->f,ae_sqr(state->epsf, _state)) )
5185 {
5186 state->repterminationtype = 1;
5187 result = ae_false;
5188 return result;
5189 }
5190
5191 /*
5192 * Main cycle
5193 */
5194 lambdaup = 10;
5195 lambdadown = 0.3;
5196 lambdav = 0.001;
5197 rho = 1;
5198 lbl_7:
5199 if( ae_false )
5200 {
5201 goto lbl_8;
5202 }
5203
5204 /*
5205 * Get Jacobian;
5206 * before we get to this point we already have State.XBase filled
5207 * with current point and State.FBase filled with function value
5208 * at XBase
5209 */
5210 nleq_clearrequestfields(state, _state);
5211 state->needfij = ae_true;
5212 ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
5213 state->rstate.stage = 2;
5214 goto lbl_rcomm;
5215 lbl_2:
5216 state->needfij = ae_false;
5217 state->repnfunc = state->repnfunc+1;
5218 state->repnjac = state->repnjac+1;
5219 rmatrixmv(n, m, &state->j, 0, 0, 1, &state->fi, 0, &state->rightpart, 0, _state);
5220 ae_v_muld(&state->rightpart.ptr.p_double[0], 1, ae_v_len(0,n-1), -1);
5221
5222 /*
5223 * Inner cycle: find good lambda
5224 */
5225 lbl_9:
5226 if( ae_false )
5227 {
5228 goto lbl_10;
5229 }
5230
5231 /*
5232 * Solve (J^T*J + (Lambda+Mu)*I)*y = J^T*F
5233 * to get step d=-y where:
5234 * * Mu=||F|| - is damping parameter for nonlinear system
5235 * * Lambda - is additional Levenberg-Marquardt parameter
5236 * for better convergence when far away from minimum
5237 */
5238 for(i=0; i<=n-1; i++)
5239 {
5240 state->candstep.ptr.p_double[i] = 0;
5241 }
5242 fblssolvecgx(&state->j, m, n, lambdav, &state->rightpart, &state->candstep, &state->cgbuf, _state);
5243
5244 /*
5245 * Normalize step (it must be no more than StpMax)
5246 */
5247 stepnorm = 0;
5248 for(i=0; i<=n-1; i++)
5249 {
5250 if( ae_fp_neq(state->candstep.ptr.p_double[i],0) )
5251 {
5252 stepnorm = 1;
5253 break;
5254 }
5255 }
5256 linminnormalized(&state->candstep, &stepnorm, n, _state);
5257 if( ae_fp_neq(state->stpmax,0) )
5258 {
5259 stepnorm = ae_minreal(stepnorm, state->stpmax, _state);
5260 }
5261
5262 /*
5263 * Test new step - is it good enough?
5264 * * if not, Lambda is increased and we try again.
5265 * * if step is good, we decrease Lambda and move on.
5266 *
5267 * We can break this cycle on two occasions:
5268 * * step is so small that x+step==x (in floating point arithmetics)
5269 * * lambda is so large
5270 */
5271 ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
5272 ae_v_addd(&state->x.ptr.p_double[0], 1, &state->candstep.ptr.p_double[0], 1, ae_v_len(0,n-1), stepnorm);
5273 b = ae_true;
5274 for(i=0; i<=n-1; i++)
5275 {
5276 if( ae_fp_neq(state->x.ptr.p_double[i],state->xbase.ptr.p_double[i]) )
5277 {
5278 b = ae_false;
5279 break;
5280 }
5281 }
5282 if( b )
5283 {
5284
5285 /*
5286 * Step is too small, force zero step and break
5287 */
5288 stepnorm = 0;
5289 ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
5290 state->f = state->fbase;
5291 goto lbl_10;
5292 }
5293 nleq_clearrequestfields(state, _state);
5294 state->needf = ae_true;
5295 state->rstate.stage = 3;
5296 goto lbl_rcomm;
5297 lbl_3:
5298 state->needf = ae_false;
5299 state->repnfunc = state->repnfunc+1;
5300 if( ae_fp_less(state->f,state->fbase) )
5301 {
5302
5303 /*
5304 * function value decreased, move on
5305 */
5306 nleq_decreaselambda(&lambdav, &rho, lambdadown, _state);
5307 goto lbl_10;
5308 }
5309 if( !nleq_increaselambda(&lambdav, &rho, lambdaup, _state) )
5310 {
5311
5312 /*
5313 * Lambda is too large (near overflow), force zero step and break
5314 */
5315 stepnorm = 0;
5316 ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
5317 state->f = state->fbase;
5318 goto lbl_10;
5319 }
5320 goto lbl_9;
5321 lbl_10:
5322
5323 /*
5324 * Accept step:
5325 * * new position
5326 * * new function value
5327 */
5328 state->fbase = state->f;
5329 ae_v_addd(&state->xbase.ptr.p_double[0], 1, &state->candstep.ptr.p_double[0], 1, ae_v_len(0,n-1), stepnorm);
5330 state->repiterationscount = state->repiterationscount+1;
5331
5332 /*
5333 * Report new iteration
5334 */
5335 if( !state->xrep )
5336 {
5337 goto lbl_11;
5338 }
5339 nleq_clearrequestfields(state, _state);
5340 state->xupdated = ae_true;
5341 state->f = state->fbase;
5342 ae_v_move(&state->x.ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,n-1));
5343 state->rstate.stage = 4;
5344 goto lbl_rcomm;
5345 lbl_4:
5346 state->xupdated = ae_false;
5347 lbl_11:
5348
5349 /*
5350 * Test stopping conditions on F, step (zero/non-zero) and MaxIts;
5351 * If one of the conditions is met, RepTerminationType is changed.
5352 */
5353 if( ae_fp_less_eq(ae_sqrt(state->f, _state),state->epsf) )
5354 {
5355 state->repterminationtype = 1;
5356 }
5357 if( ae_fp_eq(stepnorm,0)&&state->repterminationtype==0 )
5358 {
5359 state->repterminationtype = -4;
5360 }
5361 if( state->repiterationscount>=state->maxits&&state->maxits>0 )
5362 {
5363 state->repterminationtype = 5;
5364 }
5365 if( state->repterminationtype!=0 )
5366 {
5367 goto lbl_8;
5368 }
5369
5370 /*
5371 * Now, iteration is finally over
5372 */
5373 goto lbl_7;
5374 lbl_8:
5375 result = ae_false;
5376 return result;
5377
5378 /*
5379 * Saving state
5380 */
5381 lbl_rcomm:
5382 result = ae_true;
5383 state->rstate.ia.ptr.p_int[0] = n;
5384 state->rstate.ia.ptr.p_int[1] = m;
5385 state->rstate.ia.ptr.p_int[2] = i;
5386 state->rstate.ba.ptr.p_bool[0] = b;
5387 state->rstate.ra.ptr.p_double[0] = lambdaup;
5388 state->rstate.ra.ptr.p_double[1] = lambdadown;
5389 state->rstate.ra.ptr.p_double[2] = lambdav;
5390 state->rstate.ra.ptr.p_double[3] = rho;
5391 state->rstate.ra.ptr.p_double[4] = mu;
5392 state->rstate.ra.ptr.p_double[5] = stepnorm;
5393 return result;
5394 }
5395
5396
5397 /*************************************************************************
5398 NLEQ solver results
5399
5400 INPUT PARAMETERS:
5401 State - algorithm state.
5402
5403 OUTPUT PARAMETERS:
5404 X - array[0..N-1], solution
5405 Rep - optimization report:
5406 * Rep.TerminationType completetion code:
5407 * -4 ERROR: algorithm has converged to the
5408 stationary point Xf which is local minimum of
5409 f=F[0]^2+...+F[m-1]^2, but is not solution of
5410 nonlinear system.
5411 * 1 sqrt(f)<=EpsF.
5412 * 5 MaxIts steps was taken
5413 * 7 stopping conditions are too stringent,
5414 further improvement is impossible
5415 * Rep.IterationsCount contains iterations count
5416 * NFEV countains number of function calculations
5417 * ActiveConstraints contains number of active constraints
5418
5419 -- ALGLIB --
5420 Copyright 20.08.2009 by Bochkanov Sergey
5421 *************************************************************************/
nleqresults(nleqstate * state,ae_vector * x,nleqreport * rep,ae_state * _state)5422 void nleqresults(nleqstate* state,
5423 /* Real */ ae_vector* x,
5424 nleqreport* rep,
5425 ae_state *_state)
5426 {
5427
5428 ae_vector_clear(x);
5429 _nleqreport_clear(rep);
5430
5431 nleqresultsbuf(state, x, rep, _state);
5432 }
5433
5434
5435 /*************************************************************************
5436 NLEQ solver results
5437
5438 Buffered implementation of NLEQResults(), which uses pre-allocated buffer
5439 to store X[]. If buffer size is too small, it resizes buffer. It is
5440 intended to be used in the inner cycles of performance critical algorithms
5441 where array reallocation penalty is too large to be ignored.
5442
5443 -- ALGLIB --
5444 Copyright 20.08.2009 by Bochkanov Sergey
5445 *************************************************************************/
nleqresultsbuf(nleqstate * state,ae_vector * x,nleqreport * rep,ae_state * _state)5446 void nleqresultsbuf(nleqstate* state,
5447 /* Real */ ae_vector* x,
5448 nleqreport* rep,
5449 ae_state *_state)
5450 {
5451
5452
5453 if( x->cnt<state->n )
5454 {
5455 ae_vector_set_length(x, state->n, _state);
5456 }
5457 ae_v_move(&x->ptr.p_double[0], 1, &state->xbase.ptr.p_double[0], 1, ae_v_len(0,state->n-1));
5458 rep->iterationscount = state->repiterationscount;
5459 rep->nfunc = state->repnfunc;
5460 rep->njac = state->repnjac;
5461 rep->terminationtype = state->repterminationtype;
5462 }
5463
5464
5465 /*************************************************************************
5466 This subroutine restarts CG algorithm from new point. All optimization
5467 parameters are left unchanged.
5468
5469 This function allows to solve multiple optimization problems (which
5470 must have same number of dimensions) without object reallocation penalty.
5471
5472 INPUT PARAMETERS:
5473 State - structure used for reverse communication previously
5474 allocated with MinCGCreate call.
5475 X - new starting point.
5476 BndL - new lower bounds
5477 BndU - new upper bounds
5478
5479 -- ALGLIB --
5480 Copyright 30.07.2010 by Bochkanov Sergey
5481 *************************************************************************/
nleqrestartfrom(nleqstate * state,ae_vector * x,ae_state * _state)5482 void nleqrestartfrom(nleqstate* state,
5483 /* Real */ ae_vector* x,
5484 ae_state *_state)
5485 {
5486
5487
5488 ae_assert(x->cnt>=state->n, "NLEQRestartFrom: Length(X)<N!", _state);
5489 ae_assert(isfinitevector(x, state->n, _state), "NLEQRestartFrom: X contains infinite or NaN values!", _state);
5490 ae_v_move(&state->x.ptr.p_double[0], 1, &x->ptr.p_double[0], 1, ae_v_len(0,state->n-1));
5491 ae_vector_set_length(&state->rstate.ia, 2+1, _state);
5492 ae_vector_set_length(&state->rstate.ba, 0+1, _state);
5493 ae_vector_set_length(&state->rstate.ra, 5+1, _state);
5494 state->rstate.stage = -1;
5495 nleq_clearrequestfields(state, _state);
5496 }
5497
5498
5499 /*************************************************************************
5500 Clears request fileds (to be sure that we don't forgot to clear something)
5501 *************************************************************************/
nleq_clearrequestfields(nleqstate * state,ae_state * _state)5502 static void nleq_clearrequestfields(nleqstate* state, ae_state *_state)
5503 {
5504
5505
5506 state->needf = ae_false;
5507 state->needfij = ae_false;
5508 state->xupdated = ae_false;
5509 }
5510
5511
5512 /*************************************************************************
5513 Increases lambda, returns False when there is a danger of overflow
5514 *************************************************************************/
nleq_increaselambda(double * lambdav,double * nu,double lambdaup,ae_state * _state)5515 static ae_bool nleq_increaselambda(double* lambdav,
5516 double* nu,
5517 double lambdaup,
5518 ae_state *_state)
5519 {
5520 double lnlambda;
5521 double lnnu;
5522 double lnlambdaup;
5523 double lnmax;
5524 ae_bool result;
5525
5526
5527 result = ae_false;
5528 lnlambda = ae_log(*lambdav, _state);
5529 lnlambdaup = ae_log(lambdaup, _state);
5530 lnnu = ae_log(*nu, _state);
5531 lnmax = 0.5*ae_log(ae_maxrealnumber, _state);
5532 if( ae_fp_greater(lnlambda+lnlambdaup+lnnu,lnmax) )
5533 {
5534 return result;
5535 }
5536 if( ae_fp_greater(lnnu+ae_log(2, _state),lnmax) )
5537 {
5538 return result;
5539 }
5540 *lambdav = *lambdav*lambdaup*(*nu);
5541 *nu = *nu*2;
5542 result = ae_true;
5543 return result;
5544 }
5545
5546
5547 /*************************************************************************
5548 Decreases lambda, but leaves it unchanged when there is danger of underflow.
5549 *************************************************************************/
nleq_decreaselambda(double * lambdav,double * nu,double lambdadown,ae_state * _state)5550 static void nleq_decreaselambda(double* lambdav,
5551 double* nu,
5552 double lambdadown,
5553 ae_state *_state)
5554 {
5555
5556
5557 *nu = 1;
5558 if( ae_fp_less(ae_log(*lambdav, _state)+ae_log(lambdadown, _state),ae_log(ae_minrealnumber, _state)) )
5559 {
5560 *lambdav = ae_minrealnumber;
5561 }
5562 else
5563 {
5564 *lambdav = *lambdav*lambdadown;
5565 }
5566 }
5567
5568
_nleqstate_init(nleqstate * p,ae_state * _state,ae_bool make_automatic)5569 ae_bool _nleqstate_init(nleqstate* p, ae_state *_state, ae_bool make_automatic)
5570 {
5571 if( !ae_vector_init(&p->x, 0, DT_REAL, _state, make_automatic) )
5572 return ae_false;
5573 if( !ae_vector_init(&p->fi, 0, DT_REAL, _state, make_automatic) )
5574 return ae_false;
5575 if( !ae_matrix_init(&p->j, 0, 0, DT_REAL, _state, make_automatic) )
5576 return ae_false;
5577 if( !_rcommstate_init(&p->rstate, _state, make_automatic) )
5578 return ae_false;
5579 if( !ae_vector_init(&p->xbase, 0, DT_REAL, _state, make_automatic) )
5580 return ae_false;
5581 if( !ae_vector_init(&p->candstep, 0, DT_REAL, _state, make_automatic) )
5582 return ae_false;
5583 if( !ae_vector_init(&p->rightpart, 0, DT_REAL, _state, make_automatic) )
5584 return ae_false;
5585 if( !ae_vector_init(&p->cgbuf, 0, DT_REAL, _state, make_automatic) )
5586 return ae_false;
5587 return ae_true;
5588 }
5589
5590
_nleqstate_init_copy(nleqstate * dst,nleqstate * src,ae_state * _state,ae_bool make_automatic)5591 ae_bool _nleqstate_init_copy(nleqstate* dst, nleqstate* src, ae_state *_state, ae_bool make_automatic)
5592 {
5593 dst->n = src->n;
5594 dst->m = src->m;
5595 dst->epsf = src->epsf;
5596 dst->maxits = src->maxits;
5597 dst->xrep = src->xrep;
5598 dst->stpmax = src->stpmax;
5599 if( !ae_vector_init_copy(&dst->x, &src->x, _state, make_automatic) )
5600 return ae_false;
5601 dst->f = src->f;
5602 if( !ae_vector_init_copy(&dst->fi, &src->fi, _state, make_automatic) )
5603 return ae_false;
5604 if( !ae_matrix_init_copy(&dst->j, &src->j, _state, make_automatic) )
5605 return ae_false;
5606 dst->needf = src->needf;
5607 dst->needfij = src->needfij;
5608 dst->xupdated = src->xupdated;
5609 if( !_rcommstate_init_copy(&dst->rstate, &src->rstate, _state, make_automatic) )
5610 return ae_false;
5611 dst->repiterationscount = src->repiterationscount;
5612 dst->repnfunc = src->repnfunc;
5613 dst->repnjac = src->repnjac;
5614 dst->repterminationtype = src->repterminationtype;
5615 if( !ae_vector_init_copy(&dst->xbase, &src->xbase, _state, make_automatic) )
5616 return ae_false;
5617 dst->fbase = src->fbase;
5618 dst->fprev = src->fprev;
5619 if( !ae_vector_init_copy(&dst->candstep, &src->candstep, _state, make_automatic) )
5620 return ae_false;
5621 if( !ae_vector_init_copy(&dst->rightpart, &src->rightpart, _state, make_automatic) )
5622 return ae_false;
5623 if( !ae_vector_init_copy(&dst->cgbuf, &src->cgbuf, _state, make_automatic) )
5624 return ae_false;
5625 return ae_true;
5626 }
5627
5628
_nleqstate_clear(nleqstate * p)5629 void _nleqstate_clear(nleqstate* p)
5630 {
5631 ae_vector_clear(&p->x);
5632 ae_vector_clear(&p->fi);
5633 ae_matrix_clear(&p->j);
5634 _rcommstate_clear(&p->rstate);
5635 ae_vector_clear(&p->xbase);
5636 ae_vector_clear(&p->candstep);
5637 ae_vector_clear(&p->rightpart);
5638 ae_vector_clear(&p->cgbuf);
5639 }
5640
5641
_nleqreport_init(nleqreport * p,ae_state * _state,ae_bool make_automatic)5642 ae_bool _nleqreport_init(nleqreport* p, ae_state *_state, ae_bool make_automatic)
5643 {
5644 return ae_true;
5645 }
5646
5647
_nleqreport_init_copy(nleqreport * dst,nleqreport * src,ae_state * _state,ae_bool make_automatic)5648 ae_bool _nleqreport_init_copy(nleqreport* dst, nleqreport* src, ae_state *_state, ae_bool make_automatic)
5649 {
5650 dst->iterationscount = src->iterationscount;
5651 dst->nfunc = src->nfunc;
5652 dst->njac = src->njac;
5653 dst->terminationtype = src->terminationtype;
5654 return ae_true;
5655 }
5656
5657
_nleqreport_clear(nleqreport * p)5658 void _nleqreport_clear(nleqreport* p)
5659 {
5660 }
5661
5662
5663
5664 }
5665
5666