1 /************************************************************************* 2 ALGLIB 3.15.0 (source code generated 2019-02-20) 3 Copyright (c) Sergey Bochkanov (ALGLIB project). 4 5 >>> SOURCE LICENSE >>> 6 This program is free software; you can redistribute it and/or modify 7 it under the terms of the GNU General Public License as published by 8 the Free Software Foundation (www.fsf.org); either version 2 of the 9 License, or (at your option) any later version. 10 11 This program is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 GNU General Public License for more details. 15 16 A copy of the GNU General Public License is available at 17 http://www.fsf.org/licensing/licenses 18 >>> END OF LICENSE >>> 19 *************************************************************************/ 20 #ifndef _solvers_pkg_h 21 #define _solvers_pkg_h 22 #include "ap.h" 23 #include "alglibinternal.h" 24 #include "alglibmisc.h" 25 #include "linalg.h" 26 27 ///////////////////////////////////////////////////////////////////////// 28 // 29 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (DATATYPES) 30 // 31 ///////////////////////////////////////////////////////////////////////// 32 namespace alglib_impl 33 { 34 #if defined(AE_COMPILE_DIRECTDENSESOLVERS) || !defined(AE_PARTIAL_BUILD) 35 typedef struct 36 { 37 double r1; 38 double rinf; 39 } densesolverreport; 40 typedef struct 41 { 42 double r2; 43 ae_matrix cx; 44 ae_int_t n; 45 ae_int_t k; 46 } densesolverlsreport; 47 #endif 48 #if defined(AE_COMPILE_LINLSQR) || !defined(AE_PARTIAL_BUILD) 49 typedef struct 50 { 51 normestimatorstate nes; 52 ae_vector rx; 53 ae_vector b; 54 ae_int_t n; 55 ae_int_t m; 56 ae_int_t prectype; 57 ae_vector ui; 58 ae_vector uip1; 59 ae_vector vi; 60 ae_vector vip1; 61 ae_vector omegai; 62 ae_vector omegaip1; 63 double alphai; 64 double alphaip1; 65 double betai; 66 double betaip1; 67 double phibari; 68 double phibarip1; 69 double phii; 70 double rhobari; 71 double rhobarip1; 72 double rhoi; 73 double ci; 74 double si; 75 double theta; 76 double lambdai; 77 ae_vector d; 78 double anorm; 79 double bnorm2; 80 double dnorm; 81 double r2; 82 ae_vector x; 83 ae_vector mv; 84 ae_vector mtv; 85 double epsa; 86 double epsb; 87 double epsc; 88 ae_int_t maxits; 89 ae_bool xrep; 90 ae_bool xupdated; 91 ae_bool needmv; 92 ae_bool needmtv; 93 ae_bool needmv2; 94 ae_bool needvmv; 95 ae_bool needprec; 96 ae_int_t repiterationscount; 97 ae_int_t repnmv; 98 ae_int_t repterminationtype; 99 ae_bool running; 100 ae_bool userterminationneeded; 101 ae_vector tmpd; 102 ae_vector tmpx; 103 rcommstate rstate; 104 } linlsqrstate; 105 typedef struct 106 { 107 ae_int_t iterationscount; 108 ae_int_t nmv; 109 ae_int_t terminationtype; 110 } linlsqrreport; 111 #endif 112 #if defined(AE_COMPILE_POLYNOMIALSOLVER) || !defined(AE_PARTIAL_BUILD) 113 typedef struct 114 { 115 double maxerr; 116 } polynomialsolverreport; 117 #endif 118 #if defined(AE_COMPILE_NLEQ) || !defined(AE_PARTIAL_BUILD) 119 typedef struct 120 { 121 ae_int_t n; 122 ae_int_t m; 123 double epsf; 124 ae_int_t maxits; 125 ae_bool xrep; 126 double stpmax; 127 ae_vector x; 128 double f; 129 ae_vector fi; 130 ae_matrix j; 131 ae_bool needf; 132 ae_bool needfij; 133 ae_bool xupdated; 134 rcommstate rstate; 135 ae_int_t repiterationscount; 136 ae_int_t repnfunc; 137 ae_int_t repnjac; 138 ae_int_t repterminationtype; 139 ae_vector xbase; 140 double fbase; 141 double fprev; 142 ae_vector candstep; 143 ae_vector rightpart; 144 ae_vector cgbuf; 145 } nleqstate; 146 typedef struct 147 { 148 ae_int_t iterationscount; 149 ae_int_t nfunc; 150 ae_int_t njac; 151 ae_int_t terminationtype; 152 } nleqreport; 153 #endif 154 #if defined(AE_COMPILE_DIRECTSPARSESOLVERS) || !defined(AE_PARTIAL_BUILD) 155 typedef struct 156 { 157 ae_int_t terminationtype; 158 } sparsesolverreport; 159 #endif 160 #if defined(AE_COMPILE_LINCG) || !defined(AE_PARTIAL_BUILD) 161 typedef struct 162 { 163 ae_vector rx; 164 ae_vector b; 165 ae_int_t n; 166 ae_int_t prectype; 167 ae_vector cx; 168 ae_vector cr; 169 ae_vector cz; 170 ae_vector p; 171 ae_vector r; 172 ae_vector z; 173 double alpha; 174 double beta; 175 double r2; 176 double meritfunction; 177 ae_vector x; 178 ae_vector mv; 179 ae_vector pv; 180 double vmv; 181 ae_vector startx; 182 double epsf; 183 ae_int_t maxits; 184 ae_int_t itsbeforerestart; 185 ae_int_t itsbeforerupdate; 186 ae_bool xrep; 187 ae_bool xupdated; 188 ae_bool needmv; 189 ae_bool needmtv; 190 ae_bool needmv2; 191 ae_bool needvmv; 192 ae_bool needprec; 193 ae_int_t repiterationscount; 194 ae_int_t repnmv; 195 ae_int_t repterminationtype; 196 ae_bool running; 197 ae_vector tmpd; 198 rcommstate rstate; 199 } lincgstate; 200 typedef struct 201 { 202 ae_int_t iterationscount; 203 ae_int_t nmv; 204 ae_int_t terminationtype; 205 double r2; 206 } lincgreport; 207 #endif 208 209 } 210 211 ///////////////////////////////////////////////////////////////////////// 212 // 213 // THIS SECTION CONTAINS C++ INTERFACE 214 // 215 ///////////////////////////////////////////////////////////////////////// 216 namespace alglib 217 { 218 219 #if defined(AE_COMPILE_DIRECTDENSESOLVERS) || !defined(AE_PARTIAL_BUILD) 220 /************************************************************************* 221 222 *************************************************************************/ 223 class _densesolverreport_owner 224 { 225 public: 226 _densesolverreport_owner(); 227 _densesolverreport_owner(const _densesolverreport_owner &rhs); 228 _densesolverreport_owner& operator=(const _densesolverreport_owner &rhs); 229 virtual ~_densesolverreport_owner(); 230 alglib_impl::densesolverreport* c_ptr(); 231 alglib_impl::densesolverreport* c_ptr() const; 232 protected: 233 alglib_impl::densesolverreport *p_struct; 234 }; 235 class densesolverreport : public _densesolverreport_owner 236 { 237 public: 238 densesolverreport(); 239 densesolverreport(const densesolverreport &rhs); 240 densesolverreport& operator=(const densesolverreport &rhs); 241 virtual ~densesolverreport(); 242 double &r1; 243 double &rinf; 244 245 }; 246 247 248 /************************************************************************* 249 250 *************************************************************************/ 251 class _densesolverlsreport_owner 252 { 253 public: 254 _densesolverlsreport_owner(); 255 _densesolverlsreport_owner(const _densesolverlsreport_owner &rhs); 256 _densesolverlsreport_owner& operator=(const _densesolverlsreport_owner &rhs); 257 virtual ~_densesolverlsreport_owner(); 258 alglib_impl::densesolverlsreport* c_ptr(); 259 alglib_impl::densesolverlsreport* c_ptr() const; 260 protected: 261 alglib_impl::densesolverlsreport *p_struct; 262 }; 263 class densesolverlsreport : public _densesolverlsreport_owner 264 { 265 public: 266 densesolverlsreport(); 267 densesolverlsreport(const densesolverlsreport &rhs); 268 densesolverlsreport& operator=(const densesolverlsreport &rhs); 269 virtual ~densesolverlsreport(); 270 double &r2; 271 real_2d_array cx; 272 ae_int_t &n; 273 ae_int_t &k; 274 275 }; 276 #endif 277 278 #if defined(AE_COMPILE_LINLSQR) || !defined(AE_PARTIAL_BUILD) 279 /************************************************************************* 280 This object stores state of the LinLSQR method. 281 282 You should use ALGLIB functions to work with this object. 283 *************************************************************************/ 284 class _linlsqrstate_owner 285 { 286 public: 287 _linlsqrstate_owner(); 288 _linlsqrstate_owner(const _linlsqrstate_owner &rhs); 289 _linlsqrstate_owner& operator=(const _linlsqrstate_owner &rhs); 290 virtual ~_linlsqrstate_owner(); 291 alglib_impl::linlsqrstate* c_ptr(); 292 alglib_impl::linlsqrstate* c_ptr() const; 293 protected: 294 alglib_impl::linlsqrstate *p_struct; 295 }; 296 class linlsqrstate : public _linlsqrstate_owner 297 { 298 public: 299 linlsqrstate(); 300 linlsqrstate(const linlsqrstate &rhs); 301 linlsqrstate& operator=(const linlsqrstate &rhs); 302 virtual ~linlsqrstate(); 303 304 }; 305 306 307 /************************************************************************* 308 309 *************************************************************************/ 310 class _linlsqrreport_owner 311 { 312 public: 313 _linlsqrreport_owner(); 314 _linlsqrreport_owner(const _linlsqrreport_owner &rhs); 315 _linlsqrreport_owner& operator=(const _linlsqrreport_owner &rhs); 316 virtual ~_linlsqrreport_owner(); 317 alglib_impl::linlsqrreport* c_ptr(); 318 alglib_impl::linlsqrreport* c_ptr() const; 319 protected: 320 alglib_impl::linlsqrreport *p_struct; 321 }; 322 class linlsqrreport : public _linlsqrreport_owner 323 { 324 public: 325 linlsqrreport(); 326 linlsqrreport(const linlsqrreport &rhs); 327 linlsqrreport& operator=(const linlsqrreport &rhs); 328 virtual ~linlsqrreport(); 329 ae_int_t &iterationscount; 330 ae_int_t &nmv; 331 ae_int_t &terminationtype; 332 333 }; 334 #endif 335 336 #if defined(AE_COMPILE_POLYNOMIALSOLVER) || !defined(AE_PARTIAL_BUILD) 337 /************************************************************************* 338 339 *************************************************************************/ 340 class _polynomialsolverreport_owner 341 { 342 public: 343 _polynomialsolverreport_owner(); 344 _polynomialsolverreport_owner(const _polynomialsolverreport_owner &rhs); 345 _polynomialsolverreport_owner& operator=(const _polynomialsolverreport_owner &rhs); 346 virtual ~_polynomialsolverreport_owner(); 347 alglib_impl::polynomialsolverreport* c_ptr(); 348 alglib_impl::polynomialsolverreport* c_ptr() const; 349 protected: 350 alglib_impl::polynomialsolverreport *p_struct; 351 }; 352 class polynomialsolverreport : public _polynomialsolverreport_owner 353 { 354 public: 355 polynomialsolverreport(); 356 polynomialsolverreport(const polynomialsolverreport &rhs); 357 polynomialsolverreport& operator=(const polynomialsolverreport &rhs); 358 virtual ~polynomialsolverreport(); 359 double &maxerr; 360 361 }; 362 #endif 363 364 #if defined(AE_COMPILE_NLEQ) || !defined(AE_PARTIAL_BUILD) 365 /************************************************************************* 366 367 *************************************************************************/ 368 class _nleqstate_owner 369 { 370 public: 371 _nleqstate_owner(); 372 _nleqstate_owner(const _nleqstate_owner &rhs); 373 _nleqstate_owner& operator=(const _nleqstate_owner &rhs); 374 virtual ~_nleqstate_owner(); 375 alglib_impl::nleqstate* c_ptr(); 376 alglib_impl::nleqstate* c_ptr() const; 377 protected: 378 alglib_impl::nleqstate *p_struct; 379 }; 380 class nleqstate : public _nleqstate_owner 381 { 382 public: 383 nleqstate(); 384 nleqstate(const nleqstate &rhs); 385 nleqstate& operator=(const nleqstate &rhs); 386 virtual ~nleqstate(); 387 ae_bool &needf; 388 ae_bool &needfij; 389 ae_bool &xupdated; 390 double &f; 391 real_1d_array fi; 392 real_2d_array j; 393 real_1d_array x; 394 395 }; 396 397 398 /************************************************************************* 399 400 *************************************************************************/ 401 class _nleqreport_owner 402 { 403 public: 404 _nleqreport_owner(); 405 _nleqreport_owner(const _nleqreport_owner &rhs); 406 _nleqreport_owner& operator=(const _nleqreport_owner &rhs); 407 virtual ~_nleqreport_owner(); 408 alglib_impl::nleqreport* c_ptr(); 409 alglib_impl::nleqreport* c_ptr() const; 410 protected: 411 alglib_impl::nleqreport *p_struct; 412 }; 413 class nleqreport : public _nleqreport_owner 414 { 415 public: 416 nleqreport(); 417 nleqreport(const nleqreport &rhs); 418 nleqreport& operator=(const nleqreport &rhs); 419 virtual ~nleqreport(); 420 ae_int_t &iterationscount; 421 ae_int_t &nfunc; 422 ae_int_t &njac; 423 ae_int_t &terminationtype; 424 425 }; 426 #endif 427 428 #if defined(AE_COMPILE_DIRECTSPARSESOLVERS) || !defined(AE_PARTIAL_BUILD) 429 /************************************************************************* 430 This structure is a sparse solver report. 431 432 Following fields can be accessed by users: 433 *************************************************************************/ 434 class _sparsesolverreport_owner 435 { 436 public: 437 _sparsesolverreport_owner(); 438 _sparsesolverreport_owner(const _sparsesolverreport_owner &rhs); 439 _sparsesolverreport_owner& operator=(const _sparsesolverreport_owner &rhs); 440 virtual ~_sparsesolverreport_owner(); 441 alglib_impl::sparsesolverreport* c_ptr(); 442 alglib_impl::sparsesolverreport* c_ptr() const; 443 protected: 444 alglib_impl::sparsesolverreport *p_struct; 445 }; 446 class sparsesolverreport : public _sparsesolverreport_owner 447 { 448 public: 449 sparsesolverreport(); 450 sparsesolverreport(const sparsesolverreport &rhs); 451 sparsesolverreport& operator=(const sparsesolverreport &rhs); 452 virtual ~sparsesolverreport(); 453 ae_int_t &terminationtype; 454 455 }; 456 #endif 457 458 #if defined(AE_COMPILE_LINCG) || !defined(AE_PARTIAL_BUILD) 459 /************************************************************************* 460 This object stores state of the linear CG method. 461 462 You should use ALGLIB functions to work with this object. 463 Never try to access its fields directly! 464 *************************************************************************/ 465 class _lincgstate_owner 466 { 467 public: 468 _lincgstate_owner(); 469 _lincgstate_owner(const _lincgstate_owner &rhs); 470 _lincgstate_owner& operator=(const _lincgstate_owner &rhs); 471 virtual ~_lincgstate_owner(); 472 alglib_impl::lincgstate* c_ptr(); 473 alglib_impl::lincgstate* c_ptr() const; 474 protected: 475 alglib_impl::lincgstate *p_struct; 476 }; 477 class lincgstate : public _lincgstate_owner 478 { 479 public: 480 lincgstate(); 481 lincgstate(const lincgstate &rhs); 482 lincgstate& operator=(const lincgstate &rhs); 483 virtual ~lincgstate(); 484 485 }; 486 487 488 /************************************************************************* 489 490 *************************************************************************/ 491 class _lincgreport_owner 492 { 493 public: 494 _lincgreport_owner(); 495 _lincgreport_owner(const _lincgreport_owner &rhs); 496 _lincgreport_owner& operator=(const _lincgreport_owner &rhs); 497 virtual ~_lincgreport_owner(); 498 alglib_impl::lincgreport* c_ptr(); 499 alglib_impl::lincgreport* c_ptr() const; 500 protected: 501 alglib_impl::lincgreport *p_struct; 502 }; 503 class lincgreport : public _lincgreport_owner 504 { 505 public: 506 lincgreport(); 507 lincgreport(const lincgreport &rhs); 508 lincgreport& operator=(const lincgreport &rhs); 509 virtual ~lincgreport(); 510 ae_int_t &iterationscount; 511 ae_int_t &nmv; 512 ae_int_t &terminationtype; 513 double &r2; 514 515 }; 516 #endif 517 518 #if defined(AE_COMPILE_DIRECTDENSESOLVERS) || !defined(AE_PARTIAL_BUILD) 519 /************************************************************************* 520 Dense solver for A*x=b with N*N real matrix A and N*1 real vectorx x and 521 b. This is "slow-but-feature rich" version of the linear solver. Faster 522 version is RMatrixSolveFast() function. 523 524 Algorithm features: 525 * automatic detection of degenerate cases 526 * condition number estimation 527 * iterative refinement 528 * O(N^3) complexity 529 530 IMPORTANT: ! this function is NOT the most efficient linear solver provided 531 ! by ALGLIB. It estimates condition number of linear system 532 ! and performs iterative refinement, which results in 533 ! significant performance penalty when compared with "fast" 534 ! version which just performs LU decomposition and calls 535 ! triangular solver. 536 ! 537 ! This performance penalty is especially visible in the 538 ! multithreaded mode, because both condition number estimation 539 ! and iterative refinement are inherently sequential 540 ! calculations. It is also very significant on small matrices. 541 ! 542 ! Thus, if you need high performance and if you are pretty sure 543 ! that your system is well conditioned, we strongly recommend 544 ! you to use faster solver, RMatrixSolveFast() function. 545 546 ! COMMERCIAL EDITION OF ALGLIB: 547 ! 548 ! Commercial Edition of ALGLIB includes following important improvements 549 ! of this function: 550 ! * high-performance native backend with same C# interface (C# version) 551 ! * multithreading support (C++ and C# versions) 552 ! * hardware vendor (Intel) implementations of linear algebra primitives 553 ! (C++ and C# versions, x86/x64 platform) 554 ! 555 ! We recommend you to read 'Working with commercial version' section of 556 ! ALGLIB Reference Manual in order to find out how to use performance- 557 ! related features provided by commercial edition of ALGLIB. 558 559 INPUT PARAMETERS 560 A - array[0..N-1,0..N-1], system matrix 561 N - size of A 562 B - array[0..N-1], right part 563 564 OUTPUT PARAMETERS 565 Info - return code: 566 * -3 matrix is very badly conditioned or exactly singular. 567 * -1 N<=0 was passed 568 * 1 task is solved (but matrix A may be ill-conditioned, 569 check R1/RInf parameters for condition numbers). 570 Rep - additional report, following fields are set: 571 * rep.r1 condition number in 1-norm 572 * rep.rinf condition number in inf-norm 573 X - array[N], it contains: 574 * info>0 => solution 575 * info=-3 => filled by zeros 576 577 -- ALGLIB -- 578 Copyright 27.01.2010 by Bochkanov Sergey 579 *************************************************************************/ 580 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, const xparams _xparams = alglib::xdefault); 581 582 583 /************************************************************************* 584 Dense solver. 585 586 This subroutine solves a system A*x=b, where A is NxN non-denegerate 587 real matrix, x and b are vectors. This is a "fast" version of linear 588 solver which does NOT provide any additional functions like condition 589 number estimation or iterative refinement. 590 591 Algorithm features: 592 * efficient algorithm O(N^3) complexity 593 * no performance overhead from additional functionality 594 595 If you need condition number estimation or iterative refinement, use more 596 feature-rich version - RMatrixSolve(). 597 598 ! COMMERCIAL EDITION OF ALGLIB: 599 ! 600 ! Commercial Edition of ALGLIB includes following important improvements 601 ! of this function: 602 ! * high-performance native backend with same C# interface (C# version) 603 ! * multithreading support (C++ and C# versions) 604 ! * hardware vendor (Intel) implementations of linear algebra primitives 605 ! (C++ and C# versions, x86/x64 platform) 606 ! 607 ! We recommend you to read 'Working with commercial version' section of 608 ! ALGLIB Reference Manual in order to find out how to use performance- 609 ! related features provided by commercial edition of ALGLIB. 610 611 INPUT PARAMETERS 612 A - array[0..N-1,0..N-1], system matrix 613 N - size of A 614 B - array[0..N-1], right part 615 616 OUTPUT PARAMETERS 617 Info - return code: 618 * -3 matrix is exactly singular (ill conditioned matrices 619 are not recognized). 620 * -1 N<=0 was passed 621 * 1 task is solved 622 B - array[N]: 623 * info>0 => overwritten by solution 624 * info=-3 => filled by zeros 625 626 -- ALGLIB -- 627 Copyright 16.03.2015 by Bochkanov Sergey 628 *************************************************************************/ 629 void rmatrixsolvefast(const real_2d_array &a, const ae_int_t n, const real_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 630 631 632 /************************************************************************* 633 Dense solver. 634 635 Similar to RMatrixSolve() but solves task with multiple right parts (where 636 b and x are NxM matrices). This is "slow-but-robust" version of linear 637 solver with additional functionality like condition number estimation. 638 There also exists faster version - RMatrixSolveMFast(). 639 640 Algorithm features: 641 * automatic detection of degenerate cases 642 * condition number estimation 643 * optional iterative refinement 644 * O(N^3+M*N^2) complexity 645 646 IMPORTANT: ! this function is NOT the most efficient linear solver provided 647 ! by ALGLIB. It estimates condition number of linear system 648 ! and performs iterative refinement, which results in 649 ! significant performance penalty when compared with "fast" 650 ! version which just performs LU decomposition and calls 651 ! triangular solver. 652 ! 653 ! This performance penalty is especially visible in the 654 ! multithreaded mode, because both condition number estimation 655 ! and iterative refinement are inherently sequential 656 ! calculations. It also very significant on small matrices. 657 ! 658 ! Thus, if you need high performance and if you are pretty sure 659 ! that your system is well conditioned, we strongly recommend 660 ! you to use faster solver, RMatrixSolveMFast() function. 661 662 ! COMMERCIAL EDITION OF ALGLIB: 663 ! 664 ! Commercial Edition of ALGLIB includes following important improvements 665 ! of this function: 666 ! * high-performance native backend with same C# interface (C# version) 667 ! * multithreading support (C++ and C# versions) 668 ! * hardware vendor (Intel) implementations of linear algebra primitives 669 ! (C++ and C# versions, x86/x64 platform) 670 ! 671 ! We recommend you to read 'Working with commercial version' section of 672 ! ALGLIB Reference Manual in order to find out how to use performance- 673 ! related features provided by commercial edition of ALGLIB. 674 675 INPUT PARAMETERS 676 A - array[0..N-1,0..N-1], system matrix 677 N - size of A 678 B - array[0..N-1,0..M-1], right part 679 M - right part size 680 RFS - iterative refinement switch: 681 * True - refinement is used. 682 Less performance, more precision. 683 * False - refinement is not used. 684 More performance, less precision. 685 686 OUTPUT PARAMETERS 687 Info - return code: 688 * -3 A is ill conditioned or singular. 689 X is filled by zeros in such cases. 690 * -1 N<=0 was passed 691 * 1 task is solved (but matrix A may be ill-conditioned, 692 check R1/RInf parameters for condition numbers). 693 Rep - additional report, following fields are set: 694 * rep.r1 condition number in 1-norm 695 * rep.rinf condition number in inf-norm 696 X - array[N], it contains: 697 * info>0 => solution 698 * info=-3 => filled by zeros 699 700 701 -- ALGLIB -- 702 Copyright 27.01.2010 by Bochkanov Sergey 703 *************************************************************************/ 704 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, const xparams _xparams = alglib::xdefault); 705 706 707 /************************************************************************* 708 Dense solver. 709 710 Similar to RMatrixSolve() but solves task with multiple right parts (where 711 b and x are NxM matrices). This is "fast" version of linear solver which 712 does NOT offer additional functions like condition number estimation or 713 iterative refinement. 714 715 Algorithm features: 716 * O(N^3+M*N^2) complexity 717 * no additional functionality, highest performance 718 719 ! COMMERCIAL EDITION OF ALGLIB: 720 ! 721 ! Commercial Edition of ALGLIB includes following important improvements 722 ! of this function: 723 ! * high-performance native backend with same C# interface (C# version) 724 ! * multithreading support (C++ and C# versions) 725 ! * hardware vendor (Intel) implementations of linear algebra primitives 726 ! (C++ and C# versions, x86/x64 platform) 727 ! 728 ! We recommend you to read 'Working with commercial version' section of 729 ! ALGLIB Reference Manual in order to find out how to use performance- 730 ! related features provided by commercial edition of ALGLIB. 731 732 INPUT PARAMETERS 733 A - array[0..N-1,0..N-1], system matrix 734 N - size of A 735 B - array[0..N-1,0..M-1], right part 736 M - right part size 737 RFS - iterative refinement switch: 738 * True - refinement is used. 739 Less performance, more precision. 740 * False - refinement is not used. 741 More performance, less precision. 742 743 OUTPUT PARAMETERS 744 Info - return code: 745 * -3 matrix is exactly singular (ill conditioned matrices 746 are not recognized). 747 X is filled by zeros in such cases. 748 * -1 N<=0 was passed 749 * 1 task is solved 750 Rep - additional report, following fields are set: 751 * rep.r1 condition number in 1-norm 752 * rep.rinf condition number in inf-norm 753 B - array[N]: 754 * info>0 => overwritten by solution 755 * info=-3 => filled by zeros 756 757 758 -- ALGLIB -- 759 Copyright 27.01.2010 by Bochkanov Sergey 760 *************************************************************************/ 761 void rmatrixsolvemfast(const real_2d_array &a, const ae_int_t n, const real_2d_array &b, const ae_int_t m, ae_int_t &info, const xparams _xparams = alglib::xdefault); 762 763 764 /************************************************************************* 765 Dense solver. 766 767 This subroutine solves a system A*x=b, where A is NxN non-denegerate 768 real matrix given by its LU decomposition, x and b are real vectors. This 769 is "slow-but-robust" version of the linear LU-based solver. Faster version 770 is RMatrixLUSolveFast() function. 771 772 Algorithm features: 773 * automatic detection of degenerate cases 774 * O(N^2) complexity 775 * condition number estimation 776 777 No iterative refinement is provided because exact form of original matrix 778 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you 779 need iterative refinement. 780 781 IMPORTANT: ! this function is NOT the most efficient linear solver provided 782 ! by ALGLIB. It estimates condition number of linear system, 783 ! which results in 10-15x performance penalty when compared 784 ! with "fast" version which just calls triangular solver. 785 ! 786 ! This performance penalty is insignificant when compared with 787 ! cost of large LU decomposition. However, if you call this 788 ! function many times for the same left side, this overhead 789 ! BECOMES significant. It also becomes significant for small- 790 ! scale problems. 791 ! 792 ! In such cases we strongly recommend you to use faster solver, 793 ! RMatrixLUSolveFast() function. 794 795 INPUT PARAMETERS 796 LUA - array[N,N], LU decomposition, RMatrixLU result 797 P - array[N], pivots array, RMatrixLU result 798 N - size of A 799 B - array[N], right part 800 801 OUTPUT PARAMETERS 802 Info - return code: 803 * -3 matrix is very badly conditioned or exactly singular. 804 * -1 N<=0 was passed 805 * 1 task is solved (but matrix A may be ill-conditioned, 806 check R1/RInf parameters for condition numbers). 807 Rep - additional report, following fields are set: 808 * rep.r1 condition number in 1-norm 809 * rep.rinf condition number in inf-norm 810 X - array[N], it contains: 811 * info>0 => solution 812 * info=-3 => filled by zeros 813 814 815 -- ALGLIB -- 816 Copyright 27.01.2010 by Bochkanov Sergey 817 *************************************************************************/ 818 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, const xparams _xparams = alglib::xdefault); 819 820 821 /************************************************************************* 822 Dense solver. 823 824 This subroutine solves a system A*x=b, where A is NxN non-denegerate 825 real matrix given by its LU decomposition, x and b are real vectors. This 826 is "fast-without-any-checks" version of the linear LU-based solver. Slower 827 but more robust version is RMatrixLUSolve() function. 828 829 Algorithm features: 830 * O(N^2) complexity 831 * fast algorithm without ANY additional checks, just triangular solver 832 833 INPUT PARAMETERS 834 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 835 P - array[0..N-1], pivots array, RMatrixLU result 836 N - size of A 837 B - array[0..N-1], right part 838 839 OUTPUT PARAMETERS 840 Info - return code: 841 * -3 matrix is exactly singular (ill conditioned matrices 842 are not recognized). 843 X is filled by zeros in such cases. 844 * -1 N<=0 was passed 845 * 1 task is solved 846 B - array[N]: 847 * info>0 => overwritten by solution 848 * info=-3 => filled by zeros 849 850 -- ALGLIB -- 851 Copyright 18.03.2015 by Bochkanov Sergey 852 *************************************************************************/ 853 void rmatrixlusolvefast(const real_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const real_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 854 855 856 /************************************************************************* 857 Dense solver. 858 859 Similar to RMatrixLUSolve() but solves task with multiple right parts 860 (where b and x are NxM matrices). This is "robust-but-slow" version of 861 LU-based solver which performs additional checks for non-degeneracy of 862 inputs (condition number estimation). If you need best performance, use 863 "fast-without-any-checks" version, RMatrixLUSolveMFast(). 864 865 Algorithm features: 866 * automatic detection of degenerate cases 867 * O(M*N^2) complexity 868 * condition number estimation 869 870 No iterative refinement is provided because exact form of original matrix 871 is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you 872 need iterative refinement. 873 874 IMPORTANT: ! this function is NOT the most efficient linear solver provided 875 ! by ALGLIB. It estimates condition number of linear system, 876 ! which results in significant performance penalty when 877 ! compared with "fast" version which just calls triangular 878 ! solver. 879 ! 880 ! This performance penalty is especially apparent when you use 881 ! ALGLIB parallel capabilities (condition number estimation is 882 ! inherently sequential). It also becomes significant for 883 ! small-scale problems. 884 ! 885 ! In such cases we strongly recommend you to use faster solver, 886 ! RMatrixLUSolveMFast() function. 887 888 ! COMMERCIAL EDITION OF ALGLIB: 889 ! 890 ! Commercial Edition of ALGLIB includes following important improvements 891 ! of this function: 892 ! * high-performance native backend with same C# interface (C# version) 893 ! * multithreading support (C++ and C# versions) 894 ! * hardware vendor (Intel) implementations of linear algebra primitives 895 ! (C++ and C# versions, x86/x64 platform) 896 ! 897 ! We recommend you to read 'Working with commercial version' section of 898 ! ALGLIB Reference Manual in order to find out how to use performance- 899 ! related features provided by commercial edition of ALGLIB. 900 901 INPUT PARAMETERS 902 LUA - array[N,N], LU decomposition, RMatrixLU result 903 P - array[N], pivots array, RMatrixLU result 904 N - size of A 905 B - array[0..N-1,0..M-1], right part 906 M - right part size 907 908 OUTPUT PARAMETERS 909 Info - return code: 910 * -3 matrix is very badly conditioned or exactly singular. 911 X is filled by zeros in such cases. 912 * -1 N<=0 was passed 913 * 1 task is solved (but matrix A may be ill-conditioned, 914 check R1/RInf parameters for condition numbers). 915 Rep - additional report, following fields are set: 916 * rep.r1 condition number in 1-norm 917 * rep.rinf condition number in inf-norm 918 X - array[N,M], it contains: 919 * info>0 => solution 920 * info=-3 => filled by zeros 921 922 923 -- ALGLIB -- 924 Copyright 27.01.2010 by Bochkanov Sergey 925 *************************************************************************/ 926 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, const xparams _xparams = alglib::xdefault); 927 928 929 /************************************************************************* 930 Dense solver. 931 932 Similar to RMatrixLUSolve() but solves task with multiple right parts, 933 where b and x are NxM matrices. This is "fast-without-any-checks" version 934 of LU-based solver. It does not estimate condition number of a system, 935 so it is extremely fast. If you need better detection of near-degenerate 936 cases, use RMatrixLUSolveM() function. 937 938 Algorithm features: 939 * O(M*N^2) complexity 940 * fast algorithm without ANY additional checks, just triangular solver 941 942 ! COMMERCIAL EDITION OF ALGLIB: 943 ! 944 ! Commercial Edition of ALGLIB includes following important improvements 945 ! of this function: 946 ! * high-performance native backend with same C# interface (C# version) 947 ! * multithreading support (C++ and C# versions) 948 ! * hardware vendor (Intel) implementations of linear algebra primitives 949 ! (C++ and C# versions, x86/x64 platform) 950 ! 951 ! We recommend you to read 'Working with commercial version' section of 952 ! ALGLIB Reference Manual in order to find out how to use performance- 953 ! related features provided by commercial edition of ALGLIB. 954 955 INPUT PARAMETERS: 956 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 957 P - array[0..N-1], pivots array, RMatrixLU result 958 N - size of A 959 B - array[0..N-1,0..M-1], right part 960 M - right part size 961 962 OUTPUT PARAMETERS: 963 Info - return code: 964 * -3 matrix is exactly singular (ill conditioned matrices 965 are not recognized). 966 * -1 N<=0 was passed 967 * 1 task is solved 968 B - array[N,M]: 969 * info>0 => overwritten by solution 970 * info=-3 => filled by zeros 971 972 -- ALGLIB -- 973 Copyright 18.03.2015 by Bochkanov Sergey 974 *************************************************************************/ 975 void rmatrixlusolvemfast(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, const xparams _xparams = alglib::xdefault); 976 977 978 /************************************************************************* 979 Dense solver. 980 981 This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS 982 LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have 983 both A and its LU decomposition. 984 985 Algorithm features: 986 * automatic detection of degenerate cases 987 * condition number estimation 988 * iterative refinement 989 * O(N^2) complexity 990 991 INPUT PARAMETERS 992 A - array[0..N-1,0..N-1], system matrix 993 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 994 P - array[0..N-1], pivots array, RMatrixLU result 995 N - size of A 996 B - array[0..N-1], right part 997 998 OUTPUT PARAMETERS 999 Info - return code: 1000 * -3 matrix is very badly conditioned or exactly singular. 1001 * -1 N<=0 was passed 1002 * 1 task is solved (but matrix A may be ill-conditioned, 1003 check R1/RInf parameters for condition numbers). 1004 Rep - additional report, following fields are set: 1005 * rep.r1 condition number in 1-norm 1006 * rep.rinf condition number in inf-norm 1007 X - array[N], it contains: 1008 * info>0 => solution 1009 * info=-3 => filled by zeros 1010 1011 -- ALGLIB -- 1012 Copyright 27.01.2010 by Bochkanov Sergey 1013 *************************************************************************/ 1014 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, const xparams _xparams = alglib::xdefault); 1015 1016 1017 /************************************************************************* 1018 Dense solver. 1019 1020 Similar to RMatrixMixedSolve() but solves task with multiple right parts 1021 (where b and x are NxM matrices). 1022 1023 Algorithm features: 1024 * automatic detection of degenerate cases 1025 * condition number estimation 1026 * iterative refinement 1027 * O(M*N^2) complexity 1028 1029 INPUT PARAMETERS 1030 A - array[0..N-1,0..N-1], system matrix 1031 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 1032 P - array[0..N-1], pivots array, RMatrixLU result 1033 N - size of A 1034 B - array[0..N-1,0..M-1], right part 1035 M - right part size 1036 1037 OUTPUT PARAMETERS 1038 Info - return code: 1039 * -3 matrix is very badly conditioned or exactly singular. 1040 * -1 N<=0 was passed 1041 * 1 task is solved (but matrix A may be ill-conditioned, 1042 check R1/RInf parameters for condition numbers). 1043 Rep - additional report, following fields are set: 1044 * rep.r1 condition number in 1-norm 1045 * rep.rinf condition number in inf-norm 1046 X - array[N,M], it contains: 1047 * info>0 => solution 1048 * info=-3 => filled by zeros 1049 1050 -- ALGLIB -- 1051 Copyright 27.01.2010 by Bochkanov Sergey 1052 *************************************************************************/ 1053 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, const xparams _xparams = alglib::xdefault); 1054 1055 1056 /************************************************************************* 1057 Complex dense solver for A*X=B with N*N complex matrix A, N*M complex 1058 matrices X and B. "Slow-but-feature-rich" version which provides 1059 additional functions, at the cost of slower performance. Faster version 1060 may be invoked with CMatrixSolveMFast() function. 1061 1062 Algorithm features: 1063 * automatic detection of degenerate cases 1064 * condition number estimation 1065 * iterative refinement 1066 * O(N^3+M*N^2) complexity 1067 1068 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1069 ! by ALGLIB. It estimates condition number of linear system 1070 ! and performs iterative refinement, which results in 1071 ! significant performance penalty when compared with "fast" 1072 ! version which just performs LU decomposition and calls 1073 ! triangular solver. 1074 ! 1075 ! This performance penalty is especially visible in the 1076 ! multithreaded mode, because both condition number estimation 1077 ! and iterative refinement are inherently sequential 1078 ! calculations. 1079 ! 1080 ! Thus, if you need high performance and if you are pretty sure 1081 ! that your system is well conditioned, we strongly recommend 1082 ! you to use faster solver, CMatrixSolveMFast() function. 1083 1084 ! COMMERCIAL EDITION OF ALGLIB: 1085 ! 1086 ! Commercial Edition of ALGLIB includes following important improvements 1087 ! of this function: 1088 ! * high-performance native backend with same C# interface (C# version) 1089 ! * multithreading support (C++ and C# versions) 1090 ! * hardware vendor (Intel) implementations of linear algebra primitives 1091 ! (C++ and C# versions, x86/x64 platform) 1092 ! 1093 ! We recommend you to read 'Working with commercial version' section of 1094 ! ALGLIB Reference Manual in order to find out how to use performance- 1095 ! related features provided by commercial edition of ALGLIB. 1096 1097 INPUT PARAMETERS 1098 A - array[0..N-1,0..N-1], system matrix 1099 N - size of A 1100 B - array[0..N-1,0..M-1], right part 1101 M - right part size 1102 RFS - iterative refinement switch: 1103 * True - refinement is used. 1104 Less performance, more precision. 1105 * False - refinement is not used. 1106 More performance, less precision. 1107 1108 OUTPUT PARAMETERS 1109 Info - return code: 1110 * -3 matrix is very badly conditioned or exactly singular. 1111 X is filled by zeros in such cases. 1112 * -1 N<=0 was passed 1113 * 1 task is solved (but matrix A may be ill-conditioned, 1114 check R1/RInf parameters for condition numbers). 1115 Rep - additional report, following fields are set: 1116 * rep.r1 condition number in 1-norm 1117 * rep.rinf condition number in inf-norm 1118 X - array[N,M], it contains: 1119 * info>0 => solution 1120 * info=-3 => filled by zeros 1121 1122 -- ALGLIB -- 1123 Copyright 27.01.2010 by Bochkanov Sergey 1124 *************************************************************************/ 1125 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, const xparams _xparams = alglib::xdefault); 1126 1127 1128 /************************************************************************* 1129 Complex dense solver for A*X=B with N*N complex matrix A, N*M complex 1130 matrices X and B. "Fast-but-lightweight" version which provides just 1131 triangular solver - and no additional functions like iterative refinement 1132 or condition number estimation. 1133 1134 Algorithm features: 1135 * O(N^3+M*N^2) complexity 1136 * no additional time consuming functions 1137 1138 ! COMMERCIAL EDITION OF ALGLIB: 1139 ! 1140 ! Commercial Edition of ALGLIB includes following important improvements 1141 ! of this function: 1142 ! * high-performance native backend with same C# interface (C# version) 1143 ! * multithreading support (C++ and C# versions) 1144 ! * hardware vendor (Intel) implementations of linear algebra primitives 1145 ! (C++ and C# versions, x86/x64 platform) 1146 ! 1147 ! We recommend you to read 'Working with commercial version' section of 1148 ! ALGLIB Reference Manual in order to find out how to use performance- 1149 ! related features provided by commercial edition of ALGLIB. 1150 1151 INPUT PARAMETERS 1152 A - array[0..N-1,0..N-1], system matrix 1153 N - size of A 1154 B - array[0..N-1,0..M-1], right part 1155 M - right part size 1156 1157 OUTPUT PARAMETERS: 1158 Info - return code: 1159 * -3 matrix is exactly singular (ill conditioned matrices 1160 are not recognized). 1161 * -1 N<=0 was passed 1162 * 1 task is solved 1163 B - array[N,M]: 1164 * info>0 => overwritten by solution 1165 * info=-3 => filled by zeros 1166 1167 -- ALGLIB -- 1168 Copyright 16.03.2015 by Bochkanov Sergey 1169 *************************************************************************/ 1170 void cmatrixsolvemfast(const complex_2d_array &a, const ae_int_t n, const complex_2d_array &b, const ae_int_t m, ae_int_t &info, const xparams _xparams = alglib::xdefault); 1171 1172 1173 /************************************************************************* 1174 Complex dense solver for A*x=B with N*N complex matrix A and N*1 complex 1175 vectors x and b. "Slow-but-feature-rich" version of the solver. 1176 1177 Algorithm features: 1178 * automatic detection of degenerate cases 1179 * condition number estimation 1180 * iterative refinement 1181 * O(N^3) complexity 1182 1183 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1184 ! by ALGLIB. It estimates condition number of linear system 1185 ! and performs iterative refinement, which results in 1186 ! significant performance penalty when compared with "fast" 1187 ! version which just performs LU decomposition and calls 1188 ! triangular solver. 1189 ! 1190 ! This performance penalty is especially visible in the 1191 ! multithreaded mode, because both condition number estimation 1192 ! and iterative refinement are inherently sequential 1193 ! calculations. 1194 ! 1195 ! Thus, if you need high performance and if you are pretty sure 1196 ! that your system is well conditioned, we strongly recommend 1197 ! you to use faster solver, CMatrixSolveFast() function. 1198 1199 ! COMMERCIAL EDITION OF ALGLIB: 1200 ! 1201 ! Commercial Edition of ALGLIB includes following important improvements 1202 ! of this function: 1203 ! * high-performance native backend with same C# interface (C# version) 1204 ! * multithreading support (C++ and C# versions) 1205 ! * hardware vendor (Intel) implementations of linear algebra primitives 1206 ! (C++ and C# versions, x86/x64 platform) 1207 ! 1208 ! We recommend you to read 'Working with commercial version' section of 1209 ! ALGLIB Reference Manual in order to find out how to use performance- 1210 ! related features provided by commercial edition of ALGLIB. 1211 1212 INPUT PARAMETERS 1213 A - array[0..N-1,0..N-1], system matrix 1214 N - size of A 1215 B - array[0..N-1], right part 1216 1217 OUTPUT PARAMETERS 1218 Info - return code: 1219 * -3 matrix is very badly conditioned or exactly singular. 1220 * -1 N<=0 was passed 1221 * 1 task is solved (but matrix A may be ill-conditioned, 1222 check R1/RInf parameters for condition numbers). 1223 Rep - additional report, following fields are set: 1224 * rep.r1 condition number in 1-norm 1225 * rep.rinf condition number in inf-norm 1226 X - array[N], it contains: 1227 * info>0 => solution 1228 * info=-3 => filled by zeros 1229 1230 -- ALGLIB -- 1231 Copyright 27.01.2010 by Bochkanov Sergey 1232 *************************************************************************/ 1233 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, const xparams _xparams = alglib::xdefault); 1234 1235 1236 /************************************************************************* 1237 Complex dense solver for A*x=B with N*N complex matrix A and N*1 complex 1238 vectors x and b. "Fast-but-lightweight" version of the solver. 1239 1240 Algorithm features: 1241 * O(N^3) complexity 1242 * no additional time consuming features, just triangular solver 1243 1244 ! COMMERCIAL EDITION OF ALGLIB: 1245 ! 1246 ! Commercial Edition of ALGLIB includes following important improvements 1247 ! of this function: 1248 ! * high-performance native backend with same C# interface (C# version) 1249 ! * multithreading support (C++ and C# versions) 1250 ! * hardware vendor (Intel) implementations of linear algebra primitives 1251 ! (C++ and C# versions, x86/x64 platform) 1252 ! 1253 ! We recommend you to read 'Working with commercial version' section of 1254 ! ALGLIB Reference Manual in order to find out how to use performance- 1255 ! related features provided by commercial edition of ALGLIB. 1256 1257 INPUT PARAMETERS: 1258 A - array[0..N-1,0..N-1], system matrix 1259 N - size of A 1260 B - array[0..N-1], right part 1261 1262 OUTPUT PARAMETERS: 1263 Info - return code: 1264 * -3 matrix is exactly singular (ill conditioned matrices 1265 are not recognized). 1266 * -1 N<=0 was passed 1267 * 1 task is solved 1268 B - array[N]: 1269 * info>0 => overwritten by solution 1270 * info=-3 => filled by zeros 1271 1272 -- ALGLIB -- 1273 Copyright 27.01.2010 by Bochkanov Sergey 1274 *************************************************************************/ 1275 void cmatrixsolvefast(const complex_2d_array &a, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 1276 1277 1278 /************************************************************************* 1279 Dense solver for A*X=B with N*N complex A given by its LU decomposition, 1280 and N*M matrices X and B (multiple right sides). "Slow-but-feature-rich" 1281 version of the solver. 1282 1283 Algorithm features: 1284 * automatic detection of degenerate cases 1285 * O(M*N^2) complexity 1286 * condition number estimation 1287 1288 No iterative refinement is provided because exact form of original matrix 1289 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you 1290 need iterative refinement. 1291 1292 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1293 ! by ALGLIB. It estimates condition number of linear system, 1294 ! which results in significant performance penalty when 1295 ! compared with "fast" version which just calls triangular 1296 ! solver. 1297 ! 1298 ! This performance penalty is especially apparent when you use 1299 ! ALGLIB parallel capabilities (condition number estimation is 1300 ! inherently sequential). It also becomes significant for 1301 ! small-scale problems. 1302 ! 1303 ! In such cases we strongly recommend you to use faster solver, 1304 ! CMatrixLUSolveMFast() function. 1305 1306 ! COMMERCIAL EDITION OF ALGLIB: 1307 ! 1308 ! Commercial Edition of ALGLIB includes following important improvements 1309 ! of this function: 1310 ! * high-performance native backend with same C# interface (C# version) 1311 ! * multithreading support (C++ and C# versions) 1312 ! * hardware vendor (Intel) implementations of linear algebra primitives 1313 ! (C++ and C# versions, x86/x64 platform) 1314 ! 1315 ! We recommend you to read 'Working with commercial version' section of 1316 ! ALGLIB Reference Manual in order to find out how to use performance- 1317 ! related features provided by commercial edition of ALGLIB. 1318 1319 INPUT PARAMETERS 1320 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 1321 P - array[0..N-1], pivots array, RMatrixLU result 1322 N - size of A 1323 B - array[0..N-1,0..M-1], right part 1324 M - right part size 1325 1326 OUTPUT PARAMETERS 1327 Info - return code: 1328 * -3 matrix is very badly conditioned or exactly singular. 1329 * -1 N<=0 was passed 1330 * 1 task is solved (but matrix A may be ill-conditioned, 1331 check R1/RInf parameters for condition numbers). 1332 Rep - additional report, following fields are set: 1333 * rep.r1 condition number in 1-norm 1334 * rep.rinf condition number in inf-norm 1335 X - array[N,M], it contains: 1336 * info>0 => solution 1337 * info=-3 => filled by zeros 1338 1339 -- ALGLIB -- 1340 Copyright 27.01.2010 by Bochkanov Sergey 1341 *************************************************************************/ 1342 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, const xparams _xparams = alglib::xdefault); 1343 1344 1345 /************************************************************************* 1346 Dense solver for A*X=B with N*N complex A given by its LU decomposition, 1347 and N*M matrices X and B (multiple right sides). "Fast-but-lightweight" 1348 version of the solver. 1349 1350 Algorithm features: 1351 * O(M*N^2) complexity 1352 * no additional time-consuming features 1353 1354 ! COMMERCIAL EDITION OF ALGLIB: 1355 ! 1356 ! Commercial Edition of ALGLIB includes following important improvements 1357 ! of this function: 1358 ! * high-performance native backend with same C# interface (C# version) 1359 ! * multithreading support (C++ and C# versions) 1360 ! * hardware vendor (Intel) implementations of linear algebra primitives 1361 ! (C++ and C# versions, x86/x64 platform) 1362 ! 1363 ! We recommend you to read 'Working with commercial version' section of 1364 ! ALGLIB Reference Manual in order to find out how to use performance- 1365 ! related features provided by commercial edition of ALGLIB. 1366 1367 INPUT PARAMETERS 1368 LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result 1369 P - array[0..N-1], pivots array, RMatrixLU result 1370 N - size of A 1371 B - array[0..N-1,0..M-1], right part 1372 M - right part size 1373 1374 OUTPUT PARAMETERS 1375 Info - return code: 1376 * -3 matrix is exactly singular (ill conditioned matrices 1377 are not recognized). 1378 * -1 N<=0 was passed 1379 * 1 task is solved 1380 B - array[N,M]: 1381 * info>0 => overwritten by solution 1382 * info=-3 => filled by zeros 1383 1384 1385 -- ALGLIB -- 1386 Copyright 27.01.2010 by Bochkanov Sergey 1387 *************************************************************************/ 1388 void cmatrixlusolvemfast(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, const xparams _xparams = alglib::xdefault); 1389 1390 1391 /************************************************************************* 1392 Complex dense linear solver for A*x=b with complex N*N A given by its LU 1393 decomposition and N*1 vectors x and b. This is "slow-but-robust" version 1394 of the complex linear solver with additional features which add 1395 significant performance overhead. Faster version is CMatrixLUSolveFast() 1396 function. 1397 1398 Algorithm features: 1399 * automatic detection of degenerate cases 1400 * O(N^2) complexity 1401 * condition number estimation 1402 1403 No iterative refinement is provided because exact form of original matrix 1404 is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you 1405 need iterative refinement. 1406 1407 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1408 ! by ALGLIB. It estimates condition number of linear system, 1409 ! which results in 10-15x performance penalty when compared 1410 ! with "fast" version which just calls triangular solver. 1411 ! 1412 ! This performance penalty is insignificant when compared with 1413 ! cost of large LU decomposition. However, if you call this 1414 ! function many times for the same left side, this overhead 1415 ! BECOMES significant. It also becomes significant for small- 1416 ! scale problems. 1417 ! 1418 ! In such cases we strongly recommend you to use faster solver, 1419 ! CMatrixLUSolveFast() function. 1420 1421 INPUT PARAMETERS 1422 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result 1423 P - array[0..N-1], pivots array, CMatrixLU result 1424 N - size of A 1425 B - array[0..N-1], right part 1426 1427 OUTPUT PARAMETERS 1428 Info - return code: 1429 * -3 matrix is very badly conditioned or exactly singular. 1430 * -1 N<=0 was passed 1431 * 1 task is solved (but matrix A may be ill-conditioned, 1432 check R1/RInf parameters for condition numbers). 1433 Rep - additional report, following fields are set: 1434 * rep.r1 condition number in 1-norm 1435 * rep.rinf condition number in inf-norm 1436 X - array[N], it contains: 1437 * info>0 => solution 1438 * info=-3 => filled by zeros 1439 1440 -- ALGLIB -- 1441 Copyright 27.01.2010 by Bochkanov Sergey 1442 *************************************************************************/ 1443 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, const xparams _xparams = alglib::xdefault); 1444 1445 1446 /************************************************************************* 1447 Complex dense linear solver for A*x=b with N*N complex A given by its LU 1448 decomposition and N*1 vectors x and b. This is fast lightweight version 1449 of solver, which is significantly faster than CMatrixLUSolve(), but does 1450 not provide additional information (like condition numbers). 1451 1452 Algorithm features: 1453 * O(N^2) complexity 1454 * no additional time-consuming features, just triangular solver 1455 1456 INPUT PARAMETERS 1457 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result 1458 P - array[0..N-1], pivots array, CMatrixLU result 1459 N - size of A 1460 B - array[0..N-1], right part 1461 1462 OUTPUT PARAMETERS 1463 Info - return code: 1464 * -3 matrix is exactly singular (ill conditioned matrices 1465 are not recognized). 1466 * -1 N<=0 was passed 1467 * 1 task is solved 1468 B - array[N]: 1469 * info>0 => overwritten by solution 1470 * info=-3 => filled by zeros 1471 1472 NOTE: unlike CMatrixLUSolve(), this function does NOT check for 1473 near-degeneracy of input matrix. It checks for EXACT degeneracy, 1474 because this check is easy to do. However, very badly conditioned 1475 matrices may went unnoticed. 1476 1477 1478 -- ALGLIB -- 1479 Copyright 27.01.2010 by Bochkanov Sergey 1480 *************************************************************************/ 1481 void cmatrixlusolvefast(const complex_2d_array &lua, const integer_1d_array &p, const ae_int_t n, const complex_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 1482 1483 1484 /************************************************************************* 1485 Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices. 1486 1487 Algorithm features: 1488 * automatic detection of degenerate cases 1489 * condition number estimation 1490 * iterative refinement 1491 * O(M*N^2) complexity 1492 1493 INPUT PARAMETERS 1494 A - array[0..N-1,0..N-1], system matrix 1495 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result 1496 P - array[0..N-1], pivots array, CMatrixLU result 1497 N - size of A 1498 B - array[0..N-1,0..M-1], right part 1499 M - right part size 1500 1501 OUTPUT PARAMETERS 1502 Info - return code: 1503 * -3 matrix is very badly conditioned or exactly singular. 1504 * -1 N<=0 was passed 1505 * 1 task is solved (but matrix A may be ill-conditioned, 1506 check R1/RInf parameters for condition numbers). 1507 Rep - additional report, following fields are set: 1508 * rep.r1 condition number in 1-norm 1509 * rep.rinf condition number in inf-norm 1510 X - array[N,M], it contains: 1511 * info>0 => solution 1512 * info=-3 => filled by zeros 1513 1514 -- ALGLIB -- 1515 Copyright 27.01.2010 by Bochkanov Sergey 1516 *************************************************************************/ 1517 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, const xparams _xparams = alglib::xdefault); 1518 1519 1520 /************************************************************************* 1521 Dense solver. Same as RMatrixMixedSolve(), but for complex matrices. 1522 1523 Algorithm features: 1524 * automatic detection of degenerate cases 1525 * condition number estimation 1526 * iterative refinement 1527 * O(N^2) complexity 1528 1529 INPUT PARAMETERS 1530 A - array[0..N-1,0..N-1], system matrix 1531 LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result 1532 P - array[0..N-1], pivots array, CMatrixLU result 1533 N - size of A 1534 B - array[0..N-1], right part 1535 1536 OUTPUT PARAMETERS 1537 Info - return code: 1538 * -3 matrix is very badly conditioned or exactly singular. 1539 * -1 N<=0 was passed 1540 * 1 task is solved (but matrix A may be ill-conditioned, 1541 check R1/RInf parameters for condition numbers). 1542 Rep - additional report, following fields are set: 1543 * rep.r1 condition number in 1-norm 1544 * rep.rinf condition number in inf-norm 1545 X - array[N], it contains: 1546 * info>0 => solution 1547 * info=-3 => filled by zeros 1548 1549 -- ALGLIB -- 1550 Copyright 27.01.2010 by Bochkanov Sergey 1551 *************************************************************************/ 1552 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, const xparams _xparams = alglib::xdefault); 1553 1554 1555 /************************************************************************* 1556 Dense solver for A*X=B with N*N symmetric positive definite matrix A, and 1557 N*M vectors X and B. It is "slow-but-feature-rich" version of the solver. 1558 1559 Algorithm features: 1560 * automatic detection of degenerate cases 1561 * condition number estimation 1562 * O(N^3+M*N^2) complexity 1563 * matrix is represented by its upper or lower triangle 1564 1565 No iterative refinement is provided because such partial representation of 1566 matrix does not allow efficient calculation of extra-precise matrix-vector 1567 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 1568 need iterative refinement. 1569 1570 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1571 ! by ALGLIB. It estimates condition number of linear system, 1572 ! which results in significant performance penalty when 1573 ! compared with "fast" version which just performs Cholesky 1574 ! decomposition and calls triangular solver. 1575 ! 1576 ! This performance penalty is especially visible in the 1577 ! multithreaded mode, because both condition number estimation 1578 ! and iterative refinement are inherently sequential 1579 ! calculations. 1580 ! 1581 ! Thus, if you need high performance and if you are pretty sure 1582 ! that your system is well conditioned, we strongly recommend 1583 ! you to use faster solver, SPDMatrixSolveMFast() function. 1584 1585 ! COMMERCIAL EDITION OF ALGLIB: 1586 ! 1587 ! Commercial Edition of ALGLIB includes following important improvements 1588 ! of this function: 1589 ! * high-performance native backend with same C# interface (C# version) 1590 ! * multithreading support (C++ and C# versions) 1591 ! * hardware vendor (Intel) implementations of linear algebra primitives 1592 ! (C++ and C# versions, x86/x64 platform) 1593 ! 1594 ! We recommend you to read 'Working with commercial version' section of 1595 ! ALGLIB Reference Manual in order to find out how to use performance- 1596 ! related features provided by commercial edition of ALGLIB. 1597 1598 INPUT PARAMETERS 1599 A - array[0..N-1,0..N-1], system matrix 1600 N - size of A 1601 IsUpper - what half of A is provided 1602 B - array[0..N-1,0..M-1], right part 1603 M - right part size 1604 1605 OUTPUT PARAMETERS 1606 Info - return code: 1607 * -3 matrix is very badly conditioned or non-SPD. 1608 * -1 N<=0 was passed 1609 * 1 task is solved (but matrix A may be ill-conditioned, 1610 check R1/RInf parameters for condition numbers). 1611 Rep - additional report, following fields are set: 1612 * rep.r1 condition number in 1-norm 1613 * rep.rinf condition number in inf-norm 1614 X - array[N,M], it contains: 1615 * info>0 => solution 1616 * info=-3 => filled by zeros 1617 1618 -- ALGLIB -- 1619 Copyright 27.01.2010 by Bochkanov Sergey 1620 *************************************************************************/ 1621 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, const xparams _xparams = alglib::xdefault); 1622 1623 1624 /************************************************************************* 1625 Dense solver for A*X=B with N*N symmetric positive definite matrix A, and 1626 N*M vectors X and B. It is "fast-but-lightweight" version of the solver. 1627 1628 Algorithm features: 1629 * O(N^3+M*N^2) complexity 1630 * matrix is represented by its upper or lower triangle 1631 * no additional time consuming features 1632 1633 ! COMMERCIAL EDITION OF ALGLIB: 1634 ! 1635 ! Commercial Edition of ALGLIB includes following important improvements 1636 ! of this function: 1637 ! * high-performance native backend with same C# interface (C# version) 1638 ! * multithreading support (C++ and C# versions) 1639 ! * hardware vendor (Intel) implementations of linear algebra primitives 1640 ! (C++ and C# versions, x86/x64 platform) 1641 ! 1642 ! We recommend you to read 'Working with commercial version' section of 1643 ! ALGLIB Reference Manual in order to find out how to use performance- 1644 ! related features provided by commercial edition of ALGLIB. 1645 1646 INPUT PARAMETERS 1647 A - array[0..N-1,0..N-1], system matrix 1648 N - size of A 1649 IsUpper - what half of A is provided 1650 B - array[0..N-1,0..M-1], right part 1651 M - right part size 1652 1653 OUTPUT PARAMETERS 1654 Info - return code: 1655 * -3 A is is exactly singular 1656 * -1 N<=0 was passed 1657 * 1 task was solved 1658 B - array[N,M], it contains: 1659 * info>0 => solution 1660 * info=-3 => filled by zeros 1661 1662 -- ALGLIB -- 1663 Copyright 17.03.2015 by Bochkanov Sergey 1664 *************************************************************************/ 1665 void spdmatrixsolvemfast(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, const xparams _xparams = alglib::xdefault); 1666 1667 1668 /************************************************************************* 1669 Dense linear solver for A*x=b with N*N real symmetric positive definite 1670 matrix A, N*1 vectors x and b. "Slow-but-feature-rich" version of the 1671 solver. 1672 1673 Algorithm features: 1674 * automatic detection of degenerate cases 1675 * condition number estimation 1676 * O(N^3) complexity 1677 * matrix is represented by its upper or lower triangle 1678 1679 No iterative refinement is provided because such partial representation of 1680 matrix does not allow efficient calculation of extra-precise matrix-vector 1681 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 1682 need iterative refinement. 1683 1684 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1685 ! by ALGLIB. It estimates condition number of linear system, 1686 ! which results in significant performance penalty when 1687 ! compared with "fast" version which just performs Cholesky 1688 ! decomposition and calls triangular solver. 1689 ! 1690 ! This performance penalty is especially visible in the 1691 ! multithreaded mode, because both condition number estimation 1692 ! and iterative refinement are inherently sequential 1693 ! calculations. 1694 ! 1695 ! Thus, if you need high performance and if you are pretty sure 1696 ! that your system is well conditioned, we strongly recommend 1697 ! you to use faster solver, SPDMatrixSolveFast() function. 1698 1699 ! COMMERCIAL EDITION OF ALGLIB: 1700 ! 1701 ! Commercial Edition of ALGLIB includes following important improvements 1702 ! of this function: 1703 ! * high-performance native backend with same C# interface (C# version) 1704 ! * multithreading support (C++ and C# versions) 1705 ! * hardware vendor (Intel) implementations of linear algebra primitives 1706 ! (C++ and C# versions, x86/x64 platform) 1707 ! 1708 ! We recommend you to read 'Working with commercial version' section of 1709 ! ALGLIB Reference Manual in order to find out how to use performance- 1710 ! related features provided by commercial edition of ALGLIB. 1711 1712 INPUT PARAMETERS 1713 A - array[0..N-1,0..N-1], system matrix 1714 N - size of A 1715 IsUpper - what half of A is provided 1716 B - array[0..N-1], right part 1717 1718 OUTPUT PARAMETERS 1719 Info - return code: 1720 * -3 matrix is very badly conditioned or non-SPD. 1721 * -1 N<=0 was passed 1722 * 1 task is solved (but matrix A may be ill-conditioned, 1723 check R1/RInf parameters for condition numbers). 1724 Rep - additional report, following fields are set: 1725 * rep.r1 condition number in 1-norm 1726 * rep.rinf condition number in inf-norm 1727 X - array[N], it contains: 1728 * info>0 => solution 1729 * info=-3 => filled by zeros 1730 1731 -- ALGLIB -- 1732 Copyright 27.01.2010 by Bochkanov Sergey 1733 *************************************************************************/ 1734 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, const xparams _xparams = alglib::xdefault); 1735 1736 1737 /************************************************************************* 1738 Dense linear solver for A*x=b with N*N real symmetric positive definite 1739 matrix A, N*1 vectors x and b. "Fast-but-lightweight" version of the 1740 solver. 1741 1742 Algorithm features: 1743 * O(N^3) complexity 1744 * matrix is represented by its upper or lower triangle 1745 * no additional time consuming features like condition number estimation 1746 1747 ! COMMERCIAL EDITION OF ALGLIB: 1748 ! 1749 ! Commercial Edition of ALGLIB includes following important improvements 1750 ! of this function: 1751 ! * high-performance native backend with same C# interface (C# version) 1752 ! * multithreading support (C++ and C# versions) 1753 ! * hardware vendor (Intel) implementations of linear algebra primitives 1754 ! (C++ and C# versions, x86/x64 platform) 1755 ! 1756 ! We recommend you to read 'Working with commercial version' section of 1757 ! ALGLIB Reference Manual in order to find out how to use performance- 1758 ! related features provided by commercial edition of ALGLIB. 1759 1760 INPUT PARAMETERS 1761 A - array[0..N-1,0..N-1], system matrix 1762 N - size of A 1763 IsUpper - what half of A is provided 1764 B - array[0..N-1], right part 1765 1766 OUTPUT PARAMETERS 1767 Info - return code: 1768 * -3 A is is exactly singular or non-SPD 1769 * -1 N<=0 was passed 1770 * 1 task was solved 1771 B - array[N], it contains: 1772 * info>0 => solution 1773 * info=-3 => filled by zeros 1774 1775 -- ALGLIB -- 1776 Copyright 17.03.2015 by Bochkanov Sergey 1777 *************************************************************************/ 1778 void spdmatrixsolvefast(const real_2d_array &a, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 1779 1780 1781 /************************************************************************* 1782 Dense solver for A*X=B with N*N symmetric positive definite matrix A given 1783 by its Cholesky decomposition, and N*M vectors X and B. It is "slow-but- 1784 feature-rich" version of the solver which estimates condition number of 1785 the system. 1786 1787 Algorithm features: 1788 * automatic detection of degenerate cases 1789 * O(M*N^2) complexity 1790 * condition number estimation 1791 * matrix is represented by its upper or lower triangle 1792 1793 No iterative refinement is provided because such partial representation of 1794 matrix does not allow efficient calculation of extra-precise matrix-vector 1795 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 1796 need iterative refinement. 1797 1798 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1799 ! by ALGLIB. It estimates condition number of linear system, 1800 ! which results in significant performance penalty when 1801 ! compared with "fast" version which just calls triangular 1802 ! solver. Amount of overhead introduced depends on M (the 1803 ! larger - the more efficient). 1804 ! 1805 ! This performance penalty is insignificant when compared with 1806 ! cost of large LU decomposition. However, if you call this 1807 ! function many times for the same left side, this overhead 1808 ! BECOMES significant. It also becomes significant for small- 1809 ! scale problems (N<50). 1810 ! 1811 ! In such cases we strongly recommend you to use faster solver, 1812 ! SPDMatrixCholeskySolveMFast() function. 1813 1814 INPUT PARAMETERS 1815 CHA - array[0..N-1,0..N-1], Cholesky decomposition, 1816 SPDMatrixCholesky result 1817 N - size of CHA 1818 IsUpper - what half of CHA is provided 1819 B - array[0..N-1,0..M-1], right part 1820 M - right part size 1821 1822 OUTPUT PARAMETERS 1823 Info - return code: 1824 * -3 A is is exactly singular or badly conditioned 1825 X is filled by zeros in such cases. 1826 * -1 N<=0 was passed 1827 * 1 task was solved 1828 Rep - additional report, following fields are set: 1829 * rep.r1 condition number in 1-norm 1830 * rep.rinf condition number in inf-norm 1831 X - array[N]: 1832 * for info>0 contains solution 1833 * for info=-3 filled by zeros 1834 1835 -- ALGLIB -- 1836 Copyright 27.01.2010 by Bochkanov Sergey 1837 *************************************************************************/ 1838 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, const xparams _xparams = alglib::xdefault); 1839 1840 1841 /************************************************************************* 1842 Dense solver for A*X=B with N*N symmetric positive definite matrix A given 1843 by its Cholesky decomposition, and N*M vectors X and B. It is "fast-but- 1844 lightweight" version of the solver which just solves linear system, 1845 without any additional functions. 1846 1847 Algorithm features: 1848 * O(M*N^2) complexity 1849 * matrix is represented by its upper or lower triangle 1850 * no additional functionality 1851 1852 INPUT PARAMETERS 1853 CHA - array[N,N], Cholesky decomposition, 1854 SPDMatrixCholesky result 1855 N - size of CHA 1856 IsUpper - what half of CHA is provided 1857 B - array[N,M], right part 1858 M - right part size 1859 1860 OUTPUT PARAMETERS 1861 Info - return code: 1862 * -3 A is is exactly singular or badly conditioned 1863 X is filled by zeros in such cases. 1864 * -1 N<=0 was passed 1865 * 1 task was solved 1866 B - array[N]: 1867 * for info>0 overwritten by solution 1868 * for info=-3 filled by zeros 1869 1870 -- ALGLIB -- 1871 Copyright 18.03.2015 by Bochkanov Sergey 1872 *************************************************************************/ 1873 void spdmatrixcholeskysolvemfast(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, const xparams _xparams = alglib::xdefault); 1874 1875 1876 /************************************************************************* 1877 Dense solver for A*x=b with N*N symmetric positive definite matrix A given 1878 by its Cholesky decomposition, and N*1 real vectors x and b. This is "slow- 1879 but-feature-rich" version of the solver which, in addition to the 1880 solution, performs condition number estimation. 1881 1882 Algorithm features: 1883 * automatic detection of degenerate cases 1884 * O(N^2) complexity 1885 * condition number estimation 1886 * matrix is represented by its upper or lower triangle 1887 1888 No iterative refinement is provided because such partial representation of 1889 matrix does not allow efficient calculation of extra-precise matrix-vector 1890 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 1891 need iterative refinement. 1892 1893 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1894 ! by ALGLIB. It estimates condition number of linear system, 1895 ! which results in 10-15x performance penalty when compared 1896 ! with "fast" version which just calls triangular solver. 1897 ! 1898 ! This performance penalty is insignificant when compared with 1899 ! cost of large LU decomposition. However, if you call this 1900 ! function many times for the same left side, this overhead 1901 ! BECOMES significant. It also becomes significant for small- 1902 ! scale problems (N<50). 1903 ! 1904 ! In such cases we strongly recommend you to use faster solver, 1905 ! SPDMatrixCholeskySolveFast() function. 1906 1907 INPUT PARAMETERS 1908 CHA - array[N,N], Cholesky decomposition, 1909 SPDMatrixCholesky result 1910 N - size of A 1911 IsUpper - what half of CHA is provided 1912 B - array[N], right part 1913 1914 OUTPUT PARAMETERS 1915 Info - return code: 1916 * -3 A is is exactly singular or ill conditioned 1917 X is filled by zeros in such cases. 1918 * -1 N<=0 was passed 1919 * 1 task is solved 1920 Rep - additional report, following fields are set: 1921 * rep.r1 condition number in 1-norm 1922 * rep.rinf condition number in inf-norm 1923 X - array[N]: 1924 * for info>0 - solution 1925 * for info=-3 - filled by zeros 1926 1927 -- ALGLIB -- 1928 Copyright 27.01.2010 by Bochkanov Sergey 1929 *************************************************************************/ 1930 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, const xparams _xparams = alglib::xdefault); 1931 1932 1933 /************************************************************************* 1934 Dense solver for A*x=b with N*N symmetric positive definite matrix A given 1935 by its Cholesky decomposition, and N*1 real vectors x and b. This is "fast- 1936 but-lightweight" version of the solver. 1937 1938 Algorithm features: 1939 * O(N^2) complexity 1940 * matrix is represented by its upper or lower triangle 1941 * no additional features 1942 1943 INPUT PARAMETERS 1944 CHA - array[N,N], Cholesky decomposition, 1945 SPDMatrixCholesky result 1946 N - size of A 1947 IsUpper - what half of CHA is provided 1948 B - array[N], right part 1949 1950 OUTPUT PARAMETERS 1951 Info - return code: 1952 * -3 A is is exactly singular or ill conditioned 1953 X is filled by zeros in such cases. 1954 * -1 N<=0 was passed 1955 * 1 task is solved 1956 B - array[N]: 1957 * for info>0 - overwritten by solution 1958 * for info=-3 - filled by zeros 1959 1960 -- ALGLIB -- 1961 Copyright 27.01.2010 by Bochkanov Sergey 1962 *************************************************************************/ 1963 void spdmatrixcholeskysolvefast(const real_2d_array &cha, const ae_int_t n, const bool isupper, const real_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 1964 1965 1966 /************************************************************************* 1967 Dense solver for A*X=B, with N*N Hermitian positive definite matrix A and 1968 N*M complex matrices X and B. "Slow-but-feature-rich" version of the 1969 solver. 1970 1971 Algorithm features: 1972 * automatic detection of degenerate cases 1973 * condition number estimation 1974 * O(N^3+M*N^2) complexity 1975 * matrix is represented by its upper or lower triangle 1976 1977 No iterative refinement is provided because such partial representation of 1978 matrix does not allow efficient calculation of extra-precise matrix-vector 1979 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 1980 need iterative refinement. 1981 1982 IMPORTANT: ! this function is NOT the most efficient linear solver provided 1983 ! by ALGLIB. It estimates condition number of linear system, 1984 ! which results in significant performance penalty when 1985 ! compared with "fast" version which just calls triangular 1986 ! solver. 1987 ! 1988 ! This performance penalty is especially apparent when you use 1989 ! ALGLIB parallel capabilities (condition number estimation is 1990 ! inherently sequential). It also becomes significant for 1991 ! small-scale problems (N<100). 1992 ! 1993 ! In such cases we strongly recommend you to use faster solver, 1994 ! HPDMatrixSolveMFast() function. 1995 1996 ! COMMERCIAL EDITION OF ALGLIB: 1997 ! 1998 ! Commercial Edition of ALGLIB includes following important improvements 1999 ! of this function: 2000 ! * high-performance native backend with same C# interface (C# version) 2001 ! * multithreading support (C++ and C# versions) 2002 ! * hardware vendor (Intel) implementations of linear algebra primitives 2003 ! (C++ and C# versions, x86/x64 platform) 2004 ! 2005 ! We recommend you to read 'Working with commercial version' section of 2006 ! ALGLIB Reference Manual in order to find out how to use performance- 2007 ! related features provided by commercial edition of ALGLIB. 2008 2009 INPUT PARAMETERS 2010 A - array[0..N-1,0..N-1], system matrix 2011 N - size of A 2012 IsUpper - what half of A is provided 2013 B - array[0..N-1,0..M-1], right part 2014 M - right part size 2015 2016 OUTPUT PARAMETERS 2017 Info - same as in RMatrixSolve. 2018 Returns -3 for non-HPD matrices. 2019 Rep - same as in RMatrixSolve 2020 X - same as in RMatrixSolve 2021 2022 -- ALGLIB -- 2023 Copyright 27.01.2010 by Bochkanov Sergey 2024 *************************************************************************/ 2025 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, const xparams _xparams = alglib::xdefault); 2026 2027 2028 /************************************************************************* 2029 Dense solver for A*X=B, with N*N Hermitian positive definite matrix A and 2030 N*M complex matrices X and B. "Fast-but-lightweight" version of the solver. 2031 2032 Algorithm features: 2033 * O(N^3+M*N^2) complexity 2034 * matrix is represented by its upper or lower triangle 2035 * no additional time consuming features like condition number estimation 2036 2037 ! COMMERCIAL EDITION OF ALGLIB: 2038 ! 2039 ! Commercial Edition of ALGLIB includes following important improvements 2040 ! of this function: 2041 ! * high-performance native backend with same C# interface (C# version) 2042 ! * multithreading support (C++ and C# versions) 2043 ! * hardware vendor (Intel) implementations of linear algebra primitives 2044 ! (C++ and C# versions, x86/x64 platform) 2045 ! 2046 ! We recommend you to read 'Working with commercial version' section of 2047 ! ALGLIB Reference Manual in order to find out how to use performance- 2048 ! related features provided by commercial edition of ALGLIB. 2049 2050 INPUT PARAMETERS 2051 A - array[0..N-1,0..N-1], system matrix 2052 N - size of A 2053 IsUpper - what half of A is provided 2054 B - array[0..N-1,0..M-1], right part 2055 M - right part size 2056 2057 OUTPUT PARAMETERS 2058 Info - return code: 2059 * -3 A is is exactly singular or is not positive definite. 2060 B is filled by zeros in such cases. 2061 * -1 N<=0 was passed 2062 * 1 task is solved 2063 B - array[0..N-1]: 2064 * overwritten by solution 2065 * zeros, if problem was not solved 2066 2067 -- ALGLIB -- 2068 Copyright 17.03.2015 by Bochkanov Sergey 2069 *************************************************************************/ 2070 void hpdmatrixsolvemfast(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, const xparams _xparams = alglib::xdefault); 2071 2072 2073 /************************************************************************* 2074 Dense solver for A*x=b, with N*N Hermitian positive definite matrix A, and 2075 N*1 complex vectors x and b. "Slow-but-feature-rich" version of the 2076 solver. 2077 2078 Algorithm features: 2079 * automatic detection of degenerate cases 2080 * condition number estimation 2081 * O(N^3) complexity 2082 * matrix is represented by its upper or lower triangle 2083 2084 No iterative refinement is provided because such partial representation of 2085 matrix does not allow efficient calculation of extra-precise matrix-vector 2086 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 2087 need iterative refinement. 2088 2089 IMPORTANT: ! this function is NOT the most efficient linear solver provided 2090 ! by ALGLIB. It estimates condition number of linear system, 2091 ! which results in significant performance penalty when 2092 ! compared with "fast" version which just performs Cholesky 2093 ! decomposition and calls triangular solver. 2094 ! 2095 ! This performance penalty is especially visible in the 2096 ! multithreaded mode, because both condition number estimation 2097 ! and iterative refinement are inherently sequential 2098 ! calculations. 2099 ! 2100 ! Thus, if you need high performance and if you are pretty sure 2101 ! that your system is well conditioned, we strongly recommend 2102 ! you to use faster solver, HPDMatrixSolveFast() function. 2103 2104 ! COMMERCIAL EDITION OF ALGLIB: 2105 ! 2106 ! Commercial Edition of ALGLIB includes following important improvements 2107 ! of this function: 2108 ! * high-performance native backend with same C# interface (C# version) 2109 ! * multithreading support (C++ and C# versions) 2110 ! * hardware vendor (Intel) implementations of linear algebra primitives 2111 ! (C++ and C# versions, x86/x64 platform) 2112 ! 2113 ! We recommend you to read 'Working with commercial version' section of 2114 ! ALGLIB Reference Manual in order to find out how to use performance- 2115 ! related features provided by commercial edition of ALGLIB. 2116 2117 INPUT PARAMETERS 2118 A - array[0..N-1,0..N-1], system matrix 2119 N - size of A 2120 IsUpper - what half of A is provided 2121 B - array[0..N-1], right part 2122 2123 OUTPUT PARAMETERS 2124 Info - same as in RMatrixSolve 2125 Returns -3 for non-HPD matrices. 2126 Rep - same as in RMatrixSolve 2127 X - same as in RMatrixSolve 2128 2129 -- ALGLIB -- 2130 Copyright 27.01.2010 by Bochkanov Sergey 2131 *************************************************************************/ 2132 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, const xparams _xparams = alglib::xdefault); 2133 2134 2135 /************************************************************************* 2136 Dense solver for A*x=b, with N*N Hermitian positive definite matrix A, and 2137 N*1 complex vectors x and b. "Fast-but-lightweight" version of the 2138 solver without additional functions. 2139 2140 Algorithm features: 2141 * O(N^3) complexity 2142 * matrix is represented by its upper or lower triangle 2143 * no additional time consuming functions 2144 2145 ! COMMERCIAL EDITION OF ALGLIB: 2146 ! 2147 ! Commercial Edition of ALGLIB includes following important improvements 2148 ! of this function: 2149 ! * high-performance native backend with same C# interface (C# version) 2150 ! * multithreading support (C++ and C# versions) 2151 ! * hardware vendor (Intel) implementations of linear algebra primitives 2152 ! (C++ and C# versions, x86/x64 platform) 2153 ! 2154 ! We recommend you to read 'Working with commercial version' section of 2155 ! ALGLIB Reference Manual in order to find out how to use performance- 2156 ! related features provided by commercial edition of ALGLIB. 2157 2158 INPUT PARAMETERS 2159 A - array[0..N-1,0..N-1], system matrix 2160 N - size of A 2161 IsUpper - what half of A is provided 2162 B - array[0..N-1], right part 2163 2164 OUTPUT PARAMETERS 2165 Info - return code: 2166 * -3 A is is exactly singular or not positive definite 2167 X is filled by zeros in such cases. 2168 * -1 N<=0 was passed 2169 * 1 task was solved 2170 B - array[0..N-1]: 2171 * overwritten by solution 2172 * zeros, if A is exactly singular (diagonal of its LU 2173 decomposition has exact zeros). 2174 2175 -- ALGLIB -- 2176 Copyright 17.03.2015 by Bochkanov Sergey 2177 *************************************************************************/ 2178 void hpdmatrixsolvefast(const complex_2d_array &a, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 2179 2180 2181 /************************************************************************* 2182 Dense solver for A*X=B with N*N Hermitian positive definite matrix A given 2183 by its Cholesky decomposition and N*M complex matrices X and B. This is 2184 "slow-but-feature-rich" version of the solver which, in addition to the 2185 solution, estimates condition number of the system. 2186 2187 Algorithm features: 2188 * automatic detection of degenerate cases 2189 * O(M*N^2) complexity 2190 * condition number estimation 2191 * matrix is represented by its upper or lower triangle 2192 2193 No iterative refinement is provided because such partial representation of 2194 matrix does not allow efficient calculation of extra-precise matrix-vector 2195 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 2196 need iterative refinement. 2197 2198 IMPORTANT: ! this function is NOT the most efficient linear solver provided 2199 ! by ALGLIB. It estimates condition number of linear system, 2200 ! which results in significant performance penalty when 2201 ! compared with "fast" version which just calls triangular 2202 ! solver. Amount of overhead introduced depends on M (the 2203 ! larger - the more efficient). 2204 ! 2205 ! This performance penalty is insignificant when compared with 2206 ! cost of large Cholesky decomposition. However, if you call 2207 ! this function many times for the same left side, this 2208 ! overhead BECOMES significant. It also becomes significant 2209 ! for small-scale problems (N<50). 2210 ! 2211 ! In such cases we strongly recommend you to use faster solver, 2212 ! HPDMatrixCholeskySolveMFast() function. 2213 2214 2215 INPUT PARAMETERS 2216 CHA - array[N,N], Cholesky decomposition, 2217 HPDMatrixCholesky result 2218 N - size of CHA 2219 IsUpper - what half of CHA is provided 2220 B - array[N,M], right part 2221 M - right part size 2222 2223 OUTPUT PARAMETERS: 2224 Info - return code: 2225 * -3 A is singular, or VERY close to singular. 2226 X is filled by zeros in such cases. 2227 * -1 N<=0 was passed 2228 * 1 task was solved 2229 Rep - additional report, following fields are set: 2230 * rep.r1 condition number in 1-norm 2231 * rep.rinf condition number in inf-norm 2232 X - array[N]: 2233 * for info>0 contains solution 2234 * for info=-3 filled by zeros 2235 2236 -- ALGLIB -- 2237 Copyright 27.01.2010 by Bochkanov Sergey 2238 *************************************************************************/ 2239 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, const xparams _xparams = alglib::xdefault); 2240 2241 2242 /************************************************************************* 2243 Dense solver for A*X=B with N*N Hermitian positive definite matrix A given 2244 by its Cholesky decomposition and N*M complex matrices X and B. This is 2245 "fast-but-lightweight" version of the solver. 2246 2247 Algorithm features: 2248 * O(M*N^2) complexity 2249 * matrix is represented by its upper or lower triangle 2250 * no additional time-consuming features 2251 2252 INPUT PARAMETERS 2253 CHA - array[N,N], Cholesky decomposition, 2254 HPDMatrixCholesky result 2255 N - size of CHA 2256 IsUpper - what half of CHA is provided 2257 B - array[N,M], right part 2258 M - right part size 2259 2260 OUTPUT PARAMETERS: 2261 Info - return code: 2262 * -3 A is singular, or VERY close to singular. 2263 X is filled by zeros in such cases. 2264 * -1 N<=0 was passed 2265 * 1 task was solved 2266 B - array[N]: 2267 * for info>0 overwritten by solution 2268 * for info=-3 filled by zeros 2269 2270 -- ALGLIB -- 2271 Copyright 18.03.2015 by Bochkanov Sergey 2272 *************************************************************************/ 2273 void hpdmatrixcholeskysolvemfast(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, const xparams _xparams = alglib::xdefault); 2274 2275 2276 /************************************************************************* 2277 Dense solver for A*x=b with N*N Hermitian positive definite matrix A given 2278 by its Cholesky decomposition, and N*1 complex vectors x and b. This is 2279 "slow-but-feature-rich" version of the solver which estimates condition 2280 number of the system. 2281 2282 Algorithm features: 2283 * automatic detection of degenerate cases 2284 * O(N^2) complexity 2285 * condition number estimation 2286 * matrix is represented by its upper or lower triangle 2287 2288 No iterative refinement is provided because such partial representation of 2289 matrix does not allow efficient calculation of extra-precise matrix-vector 2290 products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you 2291 need iterative refinement. 2292 2293 IMPORTANT: ! this function is NOT the most efficient linear solver provided 2294 ! by ALGLIB. It estimates condition number of linear system, 2295 ! which results in 10-15x performance penalty when compared 2296 ! with "fast" version which just calls triangular solver. 2297 ! 2298 ! This performance penalty is insignificant when compared with 2299 ! cost of large LU decomposition. However, if you call this 2300 ! function many times for the same left side, this overhead 2301 ! BECOMES significant. It also becomes significant for small- 2302 ! scale problems (N<50). 2303 ! 2304 ! In such cases we strongly recommend you to use faster solver, 2305 ! HPDMatrixCholeskySolveFast() function. 2306 2307 INPUT PARAMETERS 2308 CHA - array[0..N-1,0..N-1], Cholesky decomposition, 2309 SPDMatrixCholesky result 2310 N - size of A 2311 IsUpper - what half of CHA is provided 2312 B - array[0..N-1], right part 2313 2314 OUTPUT PARAMETERS 2315 Info - return code: 2316 * -3 A is is exactly singular or ill conditioned 2317 X is filled by zeros in such cases. 2318 * -1 N<=0 was passed 2319 * 1 task is solved 2320 Rep - additional report, following fields are set: 2321 * rep.r1 condition number in 1-norm 2322 * rep.rinf condition number in inf-norm 2323 X - array[N]: 2324 * for info>0 - solution 2325 * for info=-3 - filled by zeros 2326 2327 -- ALGLIB -- 2328 Copyright 27.01.2010 by Bochkanov Sergey 2329 *************************************************************************/ 2330 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, const xparams _xparams = alglib::xdefault); 2331 2332 2333 /************************************************************************* 2334 Dense solver for A*x=b with N*N Hermitian positive definite matrix A given 2335 by its Cholesky decomposition, and N*1 complex vectors x and b. This is 2336 "fast-but-lightweight" version of the solver. 2337 2338 Algorithm features: 2339 * O(N^2) complexity 2340 * matrix is represented by its upper or lower triangle 2341 * no additional time-consuming features 2342 2343 INPUT PARAMETERS 2344 CHA - array[0..N-1,0..N-1], Cholesky decomposition, 2345 SPDMatrixCholesky result 2346 N - size of A 2347 IsUpper - what half of CHA is provided 2348 B - array[0..N-1], right part 2349 2350 OUTPUT PARAMETERS 2351 Info - return code: 2352 * -3 A is is exactly singular or ill conditioned 2353 B is filled by zeros in such cases. 2354 * -1 N<=0 was passed 2355 * 1 task is solved 2356 B - array[N]: 2357 * for info>0 - overwritten by solution 2358 * for info=-3 - filled by zeros 2359 2360 -- ALGLIB -- 2361 Copyright 18.03.2015 by Bochkanov Sergey 2362 *************************************************************************/ 2363 void hpdmatrixcholeskysolvefast(const complex_2d_array &cha, const ae_int_t n, const bool isupper, const complex_1d_array &b, ae_int_t &info, const xparams _xparams = alglib::xdefault); 2364 2365 2366 /************************************************************************* 2367 Dense solver. 2368 2369 This subroutine finds solution of the linear system A*X=B with non-square, 2370 possibly degenerate A. System is solved in the least squares sense, and 2371 general least squares solution X = X0 + CX*y which minimizes |A*X-B| is 2372 returned. If A is non-degenerate, solution in the usual sense is returned. 2373 2374 Algorithm features: 2375 * automatic detection (and correct handling!) of degenerate cases 2376 * iterative refinement 2377 * O(N^3) complexity 2378 2379 ! COMMERCIAL EDITION OF ALGLIB: 2380 ! 2381 ! Commercial Edition of ALGLIB includes following important improvements 2382 ! of this function: 2383 ! * high-performance native backend with same C# interface (C# version) 2384 ! * multithreading support (C++ and C# versions) 2385 ! * hardware vendor (Intel) implementations of linear algebra primitives 2386 ! (C++ and C# versions, x86/x64 platform) 2387 ! 2388 ! We recommend you to read 'Working with commercial version' section of 2389 ! ALGLIB Reference Manual in order to find out how to use performance- 2390 ! related features provided by commercial edition of ALGLIB. 2391 2392 INPUT PARAMETERS 2393 A - array[0..NRows-1,0..NCols-1], system matrix 2394 NRows - vertical size of A 2395 NCols - horizontal size of A 2396 B - array[0..NCols-1], right part 2397 Threshold- a number in [0,1]. Singular values beyond Threshold are 2398 considered zero. Set it to 0.0, if you don't understand 2399 what it means, so the solver will choose good value on its 2400 own. 2401 2402 OUTPUT PARAMETERS 2403 Info - return code: 2404 * -4 SVD subroutine failed 2405 * -1 if NRows<=0 or NCols<=0 or Threshold<0 was passed 2406 * 1 if task is solved 2407 Rep - solver report, see below for more info 2408 X - array[0..N-1,0..M-1], it contains: 2409 * solution of A*X=B (even for singular A) 2410 * zeros, if SVD subroutine failed 2411 2412 SOLVER REPORT 2413 2414 Subroutine sets following fields of the Rep structure: 2415 * R2 reciprocal of condition number: 1/cond(A), 2-norm. 2416 * N = NCols 2417 * K dim(Null(A)) 2418 * CX array[0..N-1,0..K-1], kernel of A. 2419 Columns of CX store such vectors that A*CX[i]=0. 2420 2421 -- ALGLIB -- 2422 Copyright 24.08.2009 by Bochkanov Sergey 2423 *************************************************************************/ 2424 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, const xparams _xparams = alglib::xdefault); 2425 #endif 2426 2427 #if defined(AE_COMPILE_LINLSQR) || !defined(AE_PARTIAL_BUILD) 2428 /************************************************************************* 2429 This function initializes linear LSQR Solver. This solver is used to solve 2430 non-symmetric (and, possibly, non-square) problems. Least squares solution 2431 is returned for non-compatible systems. 2432 2433 USAGE: 2434 1. User initializes algorithm state with LinLSQRCreate() call 2435 2. User tunes solver parameters with LinLSQRSetCond() and other functions 2436 3. User calls LinLSQRSolveSparse() function which takes algorithm state 2437 and SparseMatrix object. 2438 4. User calls LinLSQRResults() to get solution 2439 5. Optionally, user may call LinLSQRSolveSparse() again to solve another 2440 problem with different matrix and/or right part without reinitializing 2441 LinLSQRState structure. 2442 2443 INPUT PARAMETERS: 2444 M - number of rows in A 2445 N - number of variables, N>0 2446 2447 OUTPUT PARAMETERS: 2448 State - structure which stores algorithm state 2449 2450 NOTE: see also linlsqrcreatebuf() for version which reuses previously 2451 allocated place as much as possible. 2452 2453 -- ALGLIB -- 2454 Copyright 30.11.2011 by Bochkanov Sergey 2455 *************************************************************************/ 2456 void linlsqrcreate(const ae_int_t m, const ae_int_t n, linlsqrstate &state, const xparams _xparams = alglib::xdefault); 2457 2458 2459 /************************************************************************* 2460 This function initializes linear LSQR Solver. It provides exactly same 2461 functionality as linlsqrcreate(), but reuses previously allocated space 2462 as much as possible. 2463 2464 INPUT PARAMETERS: 2465 M - number of rows in A 2466 N - number of variables, N>0 2467 2468 OUTPUT PARAMETERS: 2469 State - structure which stores algorithm state 2470 2471 -- ALGLIB -- 2472 Copyright 14.11.2018 by Bochkanov Sergey 2473 *************************************************************************/ 2474 void linlsqrcreatebuf(const ae_int_t m, const ae_int_t n, const linlsqrstate &state, const xparams _xparams = alglib::xdefault); 2475 2476 2477 /************************************************************************* 2478 This function changes preconditioning settings of LinLSQQSolveSparse() 2479 function. By default, SolveSparse() uses diagonal preconditioner, but if 2480 you want to use solver without preconditioning, you can call this function 2481 which forces solver to use unit matrix for preconditioning. 2482 2483 INPUT PARAMETERS: 2484 State - structure which stores algorithm state 2485 2486 -- ALGLIB -- 2487 Copyright 19.11.2012 by Bochkanov Sergey 2488 *************************************************************************/ 2489 void linlsqrsetprecunit(const linlsqrstate &state, const xparams _xparams = alglib::xdefault); 2490 2491 2492 /************************************************************************* 2493 This function changes preconditioning settings of LinCGSolveSparse() 2494 function. LinCGSolveSparse() will use diagonal of the system matrix as 2495 preconditioner. This preconditioning mode is active by default. 2496 2497 INPUT PARAMETERS: 2498 State - structure which stores algorithm state 2499 2500 -- ALGLIB -- 2501 Copyright 19.11.2012 by Bochkanov Sergey 2502 *************************************************************************/ 2503 void linlsqrsetprecdiag(const linlsqrstate &state, const xparams _xparams = alglib::xdefault); 2504 2505 2506 /************************************************************************* 2507 This function sets optional Tikhonov regularization coefficient. 2508 It is zero by default. 2509 2510 INPUT PARAMETERS: 2511 LambdaI - regularization factor, LambdaI>=0 2512 2513 OUTPUT PARAMETERS: 2514 State - structure which stores algorithm state 2515 2516 -- ALGLIB -- 2517 Copyright 30.11.2011 by Bochkanov Sergey 2518 *************************************************************************/ 2519 void linlsqrsetlambdai(const linlsqrstate &state, const double lambdai, const xparams _xparams = alglib::xdefault); 2520 2521 2522 /************************************************************************* 2523 Procedure for solution of A*x=b with sparse A. 2524 2525 INPUT PARAMETERS: 2526 State - algorithm state 2527 A - sparse M*N matrix in the CRS format (you MUST contvert it 2528 to CRS format by calling SparseConvertToCRS() function 2529 BEFORE you pass it to this function). 2530 B - right part, array[M] 2531 2532 RESULT: 2533 This function returns no result. 2534 You can get solution by calling LinCGResults() 2535 2536 NOTE: this function uses lightweight preconditioning - multiplication by 2537 inverse of diag(A). If you want, you can turn preconditioning off by 2538 calling LinLSQRSetPrecUnit(). However, preconditioning cost is low 2539 and preconditioner is very important for solution of badly scaled 2540 problems. 2541 2542 -- ALGLIB -- 2543 Copyright 30.11.2011 by Bochkanov Sergey 2544 *************************************************************************/ 2545 void linlsqrsolvesparse(const linlsqrstate &state, const sparsematrix &a, const real_1d_array &b, const xparams _xparams = alglib::xdefault); 2546 2547 2548 /************************************************************************* 2549 This function sets stopping criteria. 2550 2551 INPUT PARAMETERS: 2552 EpsA - algorithm will be stopped if ||A^T*Rk||/(||A||*||Rk||)<=EpsA. 2553 EpsB - algorithm will be stopped if ||Rk||<=EpsB*||B|| 2554 MaxIts - algorithm will be stopped if number of iterations 2555 more than MaxIts. 2556 2557 OUTPUT PARAMETERS: 2558 State - structure which stores algorithm state 2559 2560 NOTE: if EpsA,EpsB,EpsC and MaxIts are zero then these variables will 2561 be setted as default values. 2562 2563 -- ALGLIB -- 2564 Copyright 30.11.2011 by Bochkanov Sergey 2565 *************************************************************************/ 2566 void linlsqrsetcond(const linlsqrstate &state, const double epsa, const double epsb, const ae_int_t maxits, const xparams _xparams = alglib::xdefault); 2567 2568 2569 /************************************************************************* 2570 LSQR solver: results. 2571 2572 This function must be called after LinLSQRSolve 2573 2574 INPUT PARAMETERS: 2575 State - algorithm state 2576 2577 OUTPUT PARAMETERS: 2578 X - array[N], solution 2579 Rep - optimization report: 2580 * Rep.TerminationType completetion code: 2581 * 1 ||Rk||<=EpsB*||B|| 2582 * 4 ||A^T*Rk||/(||A||*||Rk||)<=EpsA 2583 * 5 MaxIts steps was taken 2584 * 7 rounding errors prevent further progress, 2585 X contains best point found so far. 2586 (sometimes returned on singular systems) 2587 * 8 user requested termination via calling 2588 linlsqrrequesttermination() 2589 * Rep.IterationsCount contains iterations count 2590 * NMV countains number of matrix-vector calculations 2591 2592 -- ALGLIB -- 2593 Copyright 30.11.2011 by Bochkanov Sergey 2594 *************************************************************************/ 2595 void linlsqrresults(const linlsqrstate &state, real_1d_array &x, linlsqrreport &rep, const xparams _xparams = alglib::xdefault); 2596 2597 2598 /************************************************************************* 2599 This function turns on/off reporting. 2600 2601 INPUT PARAMETERS: 2602 State - structure which stores algorithm state 2603 NeedXRep- whether iteration reports are needed or not 2604 2605 If NeedXRep is True, algorithm will call rep() callback function if it is 2606 provided to MinCGOptimize(). 2607 2608 -- ALGLIB -- 2609 Copyright 30.11.2011 by Bochkanov Sergey 2610 *************************************************************************/ 2611 void linlsqrsetxrep(const linlsqrstate &state, const bool needxrep, const xparams _xparams = alglib::xdefault); 2612 2613 2614 /************************************************************************* 2615 This function is used to peek into LSQR solver and get current iteration 2616 counter. You can safely "peek" into the solver from another thread. 2617 2618 INPUT PARAMETERS: 2619 S - solver object 2620 2621 RESULT: 2622 iteration counter, in [0,INF) 2623 2624 -- ALGLIB -- 2625 Copyright 21.05.2018 by Bochkanov Sergey 2626 *************************************************************************/ 2627 ae_int_t linlsqrpeekiterationscount(const linlsqrstate &s, const xparams _xparams = alglib::xdefault); 2628 2629 2630 /************************************************************************* 2631 This subroutine submits request for termination of the running solver. It 2632 can be called from some other thread which wants LSQR solver to terminate 2633 (obviously, the thread running LSQR solver can not request termination 2634 because it is already busy working on LSQR). 2635 2636 As result, solver stops at point which was "current accepted" when 2637 termination request was submitted and returns error code 8 (successful 2638 termination). Such termination is a smooth process which properly 2639 deallocates all temporaries. 2640 2641 INPUT PARAMETERS: 2642 State - solver structure 2643 2644 NOTE: calling this function on solver which is NOT running will have no 2645 effect. 2646 2647 NOTE: multiple calls to this function are possible. First call is counted, 2648 subsequent calls are silently ignored. 2649 2650 NOTE: solver clears termination flag on its start, it means that if some 2651 other thread will request termination too soon, its request will went 2652 unnoticed. 2653 2654 -- ALGLIB -- 2655 Copyright 08.10.2014 by Bochkanov Sergey 2656 *************************************************************************/ 2657 void linlsqrrequesttermination(const linlsqrstate &state, const xparams _xparams = alglib::xdefault); 2658 #endif 2659 2660 #if defined(AE_COMPILE_POLYNOMIALSOLVER) || !defined(AE_PARTIAL_BUILD) 2661 /************************************************************************* 2662 Polynomial root finding. 2663 2664 This function returns all roots of the polynomial 2665 P(x) = a0 + a1*x + a2*x^2 + ... + an*x^n 2666 Both real and complex roots are returned (see below). 2667 2668 INPUT PARAMETERS: 2669 A - array[N+1], polynomial coefficients: 2670 * A[0] is constant term 2671 * A[N] is a coefficient of X^N 2672 N - polynomial degree 2673 2674 OUTPUT PARAMETERS: 2675 X - array of complex roots: 2676 * for isolated real root, X[I] is strictly real: IMAGE(X[I])=0 2677 * complex roots are always returned in pairs - roots occupy 2678 positions I and I+1, with: 2679 * X[I+1]=Conj(X[I]) 2680 * IMAGE(X[I]) > 0 2681 * IMAGE(X[I+1]) = -IMAGE(X[I]) < 0 2682 * multiple real roots may have non-zero imaginary part due 2683 to roundoff errors. There is no reliable way to distinguish 2684 real root of multiplicity 2 from two complex roots in 2685 the presence of roundoff errors. 2686 Rep - report, additional information, following fields are set: 2687 * Rep.MaxErr - max( |P(xi)| ) for i=0..N-1. This field 2688 allows to quickly estimate "quality" of the roots being 2689 returned. 2690 2691 NOTE: this function uses companion matrix method to find roots. In case 2692 internal EVD solver fails do find eigenvalues, exception is 2693 generated. 2694 2695 NOTE: roots are not "polished" and no matrix balancing is performed 2696 for them. 2697 2698 -- ALGLIB -- 2699 Copyright 24.02.2014 by Bochkanov Sergey 2700 *************************************************************************/ 2701 void polynomialsolve(const real_1d_array &a, const ae_int_t n, complex_1d_array &x, polynomialsolverreport &rep, const xparams _xparams = alglib::xdefault); 2702 #endif 2703 2704 #if defined(AE_COMPILE_NLEQ) || !defined(AE_PARTIAL_BUILD) 2705 /************************************************************************* 2706 LEVENBERG-MARQUARDT-LIKE NONLINEAR SOLVER 2707 2708 DESCRIPTION: 2709 This algorithm solves system of nonlinear equations 2710 F[0](x[0], ..., x[n-1]) = 0 2711 F[1](x[0], ..., x[n-1]) = 0 2712 ... 2713 F[M-1](x[0], ..., x[n-1]) = 0 2714 with M/N do not necessarily coincide. Algorithm converges quadratically 2715 under following conditions: 2716 * the solution set XS is nonempty 2717 * for some xs in XS there exist such neighbourhood N(xs) that: 2718 * vector function F(x) and its Jacobian J(x) are continuously 2719 differentiable on N 2720 * ||F(x)|| provides local error bound on N, i.e. there exists such 2721 c1, that ||F(x)||>c1*distance(x,XS) 2722 Note that these conditions are much more weaker than usual non-singularity 2723 conditions. For example, algorithm will converge for any affine function 2724 F (whether its Jacobian singular or not). 2725 2726 2727 REQUIREMENTS: 2728 Algorithm will request following information during its operation: 2729 * function vector F[] and Jacobian matrix at given point X 2730 * value of merit function f(x)=F[0]^2(x)+...+F[M-1]^2(x) at given point X 2731 2732 2733 USAGE: 2734 1. User initializes algorithm state with NLEQCreateLM() call 2735 2. User tunes solver parameters with NLEQSetCond(), NLEQSetStpMax() and 2736 other functions 2737 3. User calls NLEQSolve() function which takes algorithm state and 2738 pointers (delegates, etc.) to callback functions which calculate merit 2739 function value and Jacobian. 2740 4. User calls NLEQResults() to get solution 2741 5. Optionally, user may call NLEQRestartFrom() to solve another problem 2742 with same parameters (N/M) but another starting point and/or another 2743 function vector. NLEQRestartFrom() allows to reuse already initialized 2744 structure. 2745 2746 2747 INPUT PARAMETERS: 2748 N - space dimension, N>1: 2749 * if provided, only leading N elements of X are used 2750 * if not provided, determined automatically from size of X 2751 M - system size 2752 X - starting point 2753 2754 2755 OUTPUT PARAMETERS: 2756 State - structure which stores algorithm state 2757 2758 2759 NOTES: 2760 1. you may tune stopping conditions with NLEQSetCond() function 2761 2. if target function contains exp() or other fast growing functions, and 2762 optimization algorithm makes too large steps which leads to overflow, 2763 use NLEQSetStpMax() function to bound algorithm's steps. 2764 3. this algorithm is a slightly modified implementation of the method 2765 described in 'Levenberg-Marquardt method for constrained nonlinear 2766 equations with strong local convergence properties' by Christian Kanzow 2767 Nobuo Yamashita and Masao Fukushima and further developed in 'On the 2768 convergence of a New Levenberg-Marquardt Method' by Jin-yan Fan and 2769 Ya-Xiang Yuan. 2770 2771 2772 -- ALGLIB -- 2773 Copyright 20.08.2009 by Bochkanov Sergey 2774 *************************************************************************/ 2775 void nleqcreatelm(const ae_int_t n, const ae_int_t m, const real_1d_array &x, nleqstate &state, const xparams _xparams = alglib::xdefault); 2776 void nleqcreatelm(const ae_int_t m, const real_1d_array &x, nleqstate &state, const xparams _xparams = alglib::xdefault); 2777 2778 2779 /************************************************************************* 2780 This function sets stopping conditions for the nonlinear solver 2781 2782 INPUT PARAMETERS: 2783 State - structure which stores algorithm state 2784 EpsF - >=0 2785 The subroutine finishes its work if on k+1-th iteration 2786 the condition ||F||<=EpsF is satisfied 2787 MaxIts - maximum number of iterations. If MaxIts=0, the number of 2788 iterations is unlimited. 2789 2790 Passing EpsF=0 and MaxIts=0 simultaneously will lead to automatic 2791 stopping criterion selection (small EpsF). 2792 2793 NOTES: 2794 2795 -- ALGLIB -- 2796 Copyright 20.08.2010 by Bochkanov Sergey 2797 *************************************************************************/ 2798 void nleqsetcond(const nleqstate &state, const double epsf, const ae_int_t maxits, const xparams _xparams = alglib::xdefault); 2799 2800 2801 /************************************************************************* 2802 This function turns on/off reporting. 2803 2804 INPUT PARAMETERS: 2805 State - structure which stores algorithm state 2806 NeedXRep- whether iteration reports are needed or not 2807 2808 If NeedXRep is True, algorithm will call rep() callback function if it is 2809 provided to NLEQSolve(). 2810 2811 -- ALGLIB -- 2812 Copyright 20.08.2010 by Bochkanov Sergey 2813 *************************************************************************/ 2814 void nleqsetxrep(const nleqstate &state, const bool needxrep, const xparams _xparams = alglib::xdefault); 2815 2816 2817 /************************************************************************* 2818 This function sets maximum step length 2819 2820 INPUT PARAMETERS: 2821 State - structure which stores algorithm state 2822 StpMax - maximum step length, >=0. Set StpMax to 0.0, if you don't 2823 want to limit step length. 2824 2825 Use this subroutine when target function contains exp() or other fast 2826 growing functions, and algorithm makes too large steps which lead to 2827 overflow. This function allows us to reject steps that are too large (and 2828 therefore expose us to the possible overflow) without actually calculating 2829 function value at the x+stp*d. 2830 2831 -- ALGLIB -- 2832 Copyright 20.08.2010 by Bochkanov Sergey 2833 *************************************************************************/ 2834 void nleqsetstpmax(const nleqstate &state, const double stpmax, const xparams _xparams = alglib::xdefault); 2835 2836 2837 /************************************************************************* 2838 This function provides reverse communication interface 2839 Reverse communication interface is not documented or recommended to use. 2840 See below for functions which provide better documented API 2841 *************************************************************************/ 2842 bool nleqiteration(const nleqstate &state, const xparams _xparams = alglib::xdefault); 2843 2844 2845 /************************************************************************* 2846 This family of functions is used to launcn iterations of nonlinear solver 2847 2848 These functions accept following parameters: 2849 state - algorithm state 2850 func - callback which calculates function (or merit function) 2851 value func at given point x 2852 jac - callback which calculates function vector fi[] 2853 and Jacobian jac at given point x 2854 rep - optional callback which is called after each iteration 2855 can be NULL 2856 ptr - optional pointer which is passed to func/grad/hess/jac/rep 2857 can be NULL 2858 2859 2860 -- ALGLIB -- 2861 Copyright 20.03.2009 by Bochkanov Sergey 2862 2863 *************************************************************************/ 2864 void nleqsolve(nleqstate &state, 2865 void (*func)(const real_1d_array &x, double &func, void *ptr), 2866 void (*jac)(const real_1d_array &x, real_1d_array &fi, real_2d_array &jac, void *ptr), 2867 void (*rep)(const real_1d_array &x, double func, void *ptr) = NULL, 2868 void *ptr = NULL, 2869 const xparams _xparams = alglib::xdefault); 2870 2871 2872 /************************************************************************* 2873 NLEQ solver results 2874 2875 INPUT PARAMETERS: 2876 State - algorithm state. 2877 2878 OUTPUT PARAMETERS: 2879 X - array[0..N-1], solution 2880 Rep - optimization report: 2881 * Rep.TerminationType completetion code: 2882 * -4 ERROR: algorithm has converged to the 2883 stationary point Xf which is local minimum of 2884 f=F[0]^2+...+F[m-1]^2, but is not solution of 2885 nonlinear system. 2886 * 1 sqrt(f)<=EpsF. 2887 * 5 MaxIts steps was taken 2888 * 7 stopping conditions are too stringent, 2889 further improvement is impossible 2890 * Rep.IterationsCount contains iterations count 2891 * NFEV countains number of function calculations 2892 * ActiveConstraints contains number of active constraints 2893 2894 -- ALGLIB -- 2895 Copyright 20.08.2009 by Bochkanov Sergey 2896 *************************************************************************/ 2897 void nleqresults(const nleqstate &state, real_1d_array &x, nleqreport &rep, const xparams _xparams = alglib::xdefault); 2898 2899 2900 /************************************************************************* 2901 NLEQ solver results 2902 2903 Buffered implementation of NLEQResults(), which uses pre-allocated buffer 2904 to store X[]. If buffer size is too small, it resizes buffer. It is 2905 intended to be used in the inner cycles of performance critical algorithms 2906 where array reallocation penalty is too large to be ignored. 2907 2908 -- ALGLIB -- 2909 Copyright 20.08.2009 by Bochkanov Sergey 2910 *************************************************************************/ 2911 void nleqresultsbuf(const nleqstate &state, real_1d_array &x, nleqreport &rep, const xparams _xparams = alglib::xdefault); 2912 2913 2914 /************************************************************************* 2915 This subroutine restarts CG algorithm from new point. All optimization 2916 parameters are left unchanged. 2917 2918 This function allows to solve multiple optimization problems (which 2919 must have same number of dimensions) without object reallocation penalty. 2920 2921 INPUT PARAMETERS: 2922 State - structure used for reverse communication previously 2923 allocated with MinCGCreate call. 2924 X - new starting point. 2925 BndL - new lower bounds 2926 BndU - new upper bounds 2927 2928 -- ALGLIB -- 2929 Copyright 30.07.2010 by Bochkanov Sergey 2930 *************************************************************************/ 2931 void nleqrestartfrom(const nleqstate &state, const real_1d_array &x, const xparams _xparams = alglib::xdefault); 2932 #endif 2933 2934 #if defined(AE_COMPILE_DIRECTSPARSESOLVERS) || !defined(AE_PARTIAL_BUILD) 2935 /************************************************************************* 2936 Sparse linear solver for A*x=b with N*N sparse real symmetric positive 2937 definite matrix A, N*1 vectors x and b. 2938 2939 This solver converts input matrix to SKS format, performs Cholesky 2940 factorization using SKS Cholesky subroutine (works well for limited 2941 bandwidth matrices) and uses sparse triangular solvers to get solution of 2942 the original system. 2943 2944 INPUT PARAMETERS 2945 A - sparse matrix, must be NxN exactly 2946 N - size of A, N>0 2947 IsUpper - which half of A is provided (another half is ignored) 2948 B - array[0..N-1], right part 2949 2950 OUTPUT PARAMETERS 2951 Rep - solver report, following fields are set: 2952 * rep.terminationtype - solver status; >0 for success, 2953 set to -3 on failure (degenerate or non-SPD system). 2954 X - array[N], it contains: 2955 * rep.terminationtype>0 => solution 2956 * rep.terminationtype=-3 => filled by zeros 2957 2958 -- ALGLIB -- 2959 Copyright 26.12.2017 by Bochkanov Sergey 2960 *************************************************************************/ 2961 void sparsesolvesks(const sparsematrix &a, const ae_int_t n, const bool isupper, const real_1d_array &b, sparsesolverreport &rep, real_1d_array &x, const xparams _xparams = alglib::xdefault); 2962 2963 2964 /************************************************************************* 2965 Sparse linear solver for A*x=b with N*N real symmetric positive definite 2966 matrix A given by its Cholesky decomposition, and N*1 vectors x and b. 2967 2968 IMPORTANT: this solver requires input matrix to be in the SKS (Skyline) 2969 sparse storage format. An exception will be generated if you 2970 pass matrix in some other format (HASH or CRS). 2971 2972 INPUT PARAMETERS 2973 A - sparse NxN matrix stored in SKS format, must be NxN exactly 2974 N - size of A, N>0 2975 IsUpper - which half of A is provided (another half is ignored) 2976 B - array[N], right part 2977 2978 OUTPUT PARAMETERS 2979 Rep - solver report, following fields are set: 2980 * rep.terminationtype - solver status; >0 for success, 2981 set to -3 on failure (degenerate or non-SPD system). 2982 X - array[N], it contains: 2983 * rep.terminationtype>0 => solution 2984 * rep.terminationtype=-3 => filled by zeros 2985 2986 -- ALGLIB -- 2987 Copyright 26.12.2017 by Bochkanov Sergey 2988 *************************************************************************/ 2989 void sparsecholeskysolvesks(const sparsematrix &a, const ae_int_t n, const bool isupper, const real_1d_array &b, sparsesolverreport &rep, real_1d_array &x, const xparams _xparams = alglib::xdefault); 2990 #endif 2991 2992 #if defined(AE_COMPILE_LINCG) || !defined(AE_PARTIAL_BUILD) 2993 /************************************************************************* 2994 This function initializes linear CG Solver. This solver is used to solve 2995 symmetric positive definite problems. If you want to solve nonsymmetric 2996 (or non-positive definite) problem you may use LinLSQR solver provided by 2997 ALGLIB. 2998 2999 USAGE: 3000 1. User initializes algorithm state with LinCGCreate() call 3001 2. User tunes solver parameters with LinCGSetCond() and other functions 3002 3. Optionally, user sets starting point with LinCGSetStartingPoint() 3003 4. User calls LinCGSolveSparse() function which takes algorithm state and 3004 SparseMatrix object. 3005 5. User calls LinCGResults() to get solution 3006 6. Optionally, user may call LinCGSolveSparse() again to solve another 3007 problem with different matrix and/or right part without reinitializing 3008 LinCGState structure. 3009 3010 INPUT PARAMETERS: 3011 N - problem dimension, N>0 3012 3013 OUTPUT PARAMETERS: 3014 State - structure which stores algorithm state 3015 3016 -- ALGLIB -- 3017 Copyright 14.11.2011 by Bochkanov Sergey 3018 *************************************************************************/ 3019 void lincgcreate(const ae_int_t n, lincgstate &state, const xparams _xparams = alglib::xdefault); 3020 3021 3022 /************************************************************************* 3023 This function sets starting point. 3024 By default, zero starting point is used. 3025 3026 INPUT PARAMETERS: 3027 X - starting point, array[N] 3028 3029 OUTPUT PARAMETERS: 3030 State - structure which stores algorithm state 3031 3032 -- ALGLIB -- 3033 Copyright 14.11.2011 by Bochkanov Sergey 3034 *************************************************************************/ 3035 void lincgsetstartingpoint(const lincgstate &state, const real_1d_array &x, const xparams _xparams = alglib::xdefault); 3036 3037 3038 /************************************************************************* 3039 This function changes preconditioning settings of LinCGSolveSparse() 3040 function. By default, SolveSparse() uses diagonal preconditioner, but if 3041 you want to use solver without preconditioning, you can call this function 3042 which forces solver to use unit matrix for preconditioning. 3043 3044 INPUT PARAMETERS: 3045 State - structure which stores algorithm state 3046 3047 -- ALGLIB -- 3048 Copyright 19.11.2012 by Bochkanov Sergey 3049 *************************************************************************/ 3050 void lincgsetprecunit(const lincgstate &state, const xparams _xparams = alglib::xdefault); 3051 3052 3053 /************************************************************************* 3054 This function changes preconditioning settings of LinCGSolveSparse() 3055 function. LinCGSolveSparse() will use diagonal of the system matrix as 3056 preconditioner. This preconditioning mode is active by default. 3057 3058 INPUT PARAMETERS: 3059 State - structure which stores algorithm state 3060 3061 -- ALGLIB -- 3062 Copyright 19.11.2012 by Bochkanov Sergey 3063 *************************************************************************/ 3064 void lincgsetprecdiag(const lincgstate &state, const xparams _xparams = alglib::xdefault); 3065 3066 3067 /************************************************************************* 3068 This function sets stopping criteria. 3069 3070 INPUT PARAMETERS: 3071 EpsF - algorithm will be stopped if norm of residual is less than 3072 EpsF*||b||. 3073 MaxIts - algorithm will be stopped if number of iterations is more 3074 than MaxIts. 3075 3076 OUTPUT PARAMETERS: 3077 State - structure which stores algorithm state 3078 3079 NOTES: 3080 If both EpsF and MaxIts are zero then small EpsF will be set to small 3081 value. 3082 3083 -- ALGLIB -- 3084 Copyright 14.11.2011 by Bochkanov Sergey 3085 *************************************************************************/ 3086 void lincgsetcond(const lincgstate &state, const double epsf, const ae_int_t maxits, const xparams _xparams = alglib::xdefault); 3087 3088 3089 /************************************************************************* 3090 Procedure for solution of A*x=b with sparse A. 3091 3092 INPUT PARAMETERS: 3093 State - algorithm state 3094 A - sparse matrix in the CRS format (you MUST contvert it to 3095 CRS format by calling SparseConvertToCRS() function). 3096 IsUpper - whether upper or lower triangle of A is used: 3097 * IsUpper=True => only upper triangle is used and lower 3098 triangle is not referenced at all 3099 * IsUpper=False => only lower triangle is used and upper 3100 triangle is not referenced at all 3101 B - right part, array[N] 3102 3103 RESULT: 3104 This function returns no result. 3105 You can get solution by calling LinCGResults() 3106 3107 NOTE: this function uses lightweight preconditioning - multiplication by 3108 inverse of diag(A). If you want, you can turn preconditioning off by 3109 calling LinCGSetPrecUnit(). However, preconditioning cost is low and 3110 preconditioner is very important for solution of badly scaled 3111 problems. 3112 3113 -- ALGLIB -- 3114 Copyright 14.11.2011 by Bochkanov Sergey 3115 *************************************************************************/ 3116 void lincgsolvesparse(const lincgstate &state, const sparsematrix &a, const bool isupper, const real_1d_array &b, const xparams _xparams = alglib::xdefault); 3117 3118 3119 /************************************************************************* 3120 CG-solver: results. 3121 3122 This function must be called after LinCGSolve 3123 3124 INPUT PARAMETERS: 3125 State - algorithm state 3126 3127 OUTPUT PARAMETERS: 3128 X - array[N], solution 3129 Rep - optimization report: 3130 * Rep.TerminationType completetion code: 3131 * -5 input matrix is either not positive definite, 3132 too large or too small 3133 * -4 overflow/underflow during solution 3134 (ill conditioned problem) 3135 * 1 ||residual||<=EpsF*||b|| 3136 * 5 MaxIts steps was taken 3137 * 7 rounding errors prevent further progress, 3138 best point found is returned 3139 * Rep.IterationsCount contains iterations count 3140 * NMV countains number of matrix-vector calculations 3141 3142 -- ALGLIB -- 3143 Copyright 14.11.2011 by Bochkanov Sergey 3144 *************************************************************************/ 3145 void lincgresults(const lincgstate &state, real_1d_array &x, lincgreport &rep, const xparams _xparams = alglib::xdefault); 3146 3147 3148 /************************************************************************* 3149 This function sets restart frequency. By default, algorithm is restarted 3150 after N subsequent iterations. 3151 3152 -- ALGLIB -- 3153 Copyright 14.11.2011 by Bochkanov Sergey 3154 *************************************************************************/ 3155 void lincgsetrestartfreq(const lincgstate &state, const ae_int_t srf, const xparams _xparams = alglib::xdefault); 3156 3157 3158 /************************************************************************* 3159 This function sets frequency of residual recalculations. 3160 3161 Algorithm updates residual r_k using iterative formula, but recalculates 3162 it from scratch after each 10 iterations. It is done to avoid accumulation 3163 of numerical errors and to stop algorithm when r_k starts to grow. 3164 3165 Such low update frequence (1/10) gives very little overhead, but makes 3166 algorithm a bit more robust against numerical errors. However, you may 3167 change it 3168 3169 INPUT PARAMETERS: 3170 Freq - desired update frequency, Freq>=0. 3171 Zero value means that no updates will be done. 3172 3173 -- ALGLIB -- 3174 Copyright 14.11.2011 by Bochkanov Sergey 3175 *************************************************************************/ 3176 void lincgsetrupdatefreq(const lincgstate &state, const ae_int_t freq, const xparams _xparams = alglib::xdefault); 3177 3178 3179 /************************************************************************* 3180 This function turns on/off reporting. 3181 3182 INPUT PARAMETERS: 3183 State - structure which stores algorithm state 3184 NeedXRep- whether iteration reports are needed or not 3185 3186 If NeedXRep is True, algorithm will call rep() callback function if it is 3187 provided to MinCGOptimize(). 3188 3189 -- ALGLIB -- 3190 Copyright 14.11.2011 by Bochkanov Sergey 3191 *************************************************************************/ 3192 void lincgsetxrep(const lincgstate &state, const bool needxrep, const xparams _xparams = alglib::xdefault); 3193 #endif 3194 } 3195 3196 ///////////////////////////////////////////////////////////////////////// 3197 // 3198 // THIS SECTION CONTAINS COMPUTATIONAL CORE DECLARATIONS (FUNCTIONS) 3199 // 3200 ///////////////////////////////////////////////////////////////////////// 3201 namespace alglib_impl 3202 { 3203 #if defined(AE_COMPILE_DIRECTDENSESOLVERS) || !defined(AE_PARTIAL_BUILD) 3204 void rmatrixsolve(/* Real */ ae_matrix* a, 3205 ae_int_t n, 3206 /* Real */ ae_vector* b, 3207 ae_int_t* info, 3208 densesolverreport* rep, 3209 /* Real */ ae_vector* x, 3210 ae_state *_state); 3211 void rmatrixsolvefast(/* Real */ ae_matrix* a, 3212 ae_int_t n, 3213 /* Real */ ae_vector* b, 3214 ae_int_t* info, 3215 ae_state *_state); 3216 void rmatrixsolvem(/* Real */ ae_matrix* a, 3217 ae_int_t n, 3218 /* Real */ ae_matrix* b, 3219 ae_int_t m, 3220 ae_bool rfs, 3221 ae_int_t* info, 3222 densesolverreport* rep, 3223 /* Real */ ae_matrix* x, 3224 ae_state *_state); 3225 void rmatrixsolvemfast(/* Real */ ae_matrix* a, 3226 ae_int_t n, 3227 /* Real */ ae_matrix* b, 3228 ae_int_t m, 3229 ae_int_t* info, 3230 ae_state *_state); 3231 void rmatrixlusolve(/* Real */ ae_matrix* lua, 3232 /* Integer */ ae_vector* p, 3233 ae_int_t n, 3234 /* Real */ ae_vector* b, 3235 ae_int_t* info, 3236 densesolverreport* rep, 3237 /* Real */ ae_vector* x, 3238 ae_state *_state); 3239 void rmatrixlusolvefast(/* Real */ ae_matrix* lua, 3240 /* Integer */ ae_vector* p, 3241 ae_int_t n, 3242 /* Real */ ae_vector* b, 3243 ae_int_t* info, 3244 ae_state *_state); 3245 void rmatrixlusolvem(/* Real */ ae_matrix* lua, 3246 /* Integer */ ae_vector* p, 3247 ae_int_t n, 3248 /* Real */ ae_matrix* b, 3249 ae_int_t m, 3250 ae_int_t* info, 3251 densesolverreport* rep, 3252 /* Real */ ae_matrix* x, 3253 ae_state *_state); 3254 void rmatrixlusolvemfast(/* Real */ ae_matrix* lua, 3255 /* Integer */ ae_vector* p, 3256 ae_int_t n, 3257 /* Real */ ae_matrix* b, 3258 ae_int_t m, 3259 ae_int_t* info, 3260 ae_state *_state); 3261 void rmatrixmixedsolve(/* Real */ ae_matrix* a, 3262 /* Real */ ae_matrix* lua, 3263 /* Integer */ ae_vector* p, 3264 ae_int_t n, 3265 /* Real */ ae_vector* b, 3266 ae_int_t* info, 3267 densesolverreport* rep, 3268 /* Real */ ae_vector* x, 3269 ae_state *_state); 3270 void rmatrixmixedsolvem(/* Real */ ae_matrix* a, 3271 /* Real */ ae_matrix* lua, 3272 /* Integer */ ae_vector* p, 3273 ae_int_t n, 3274 /* Real */ ae_matrix* b, 3275 ae_int_t m, 3276 ae_int_t* info, 3277 densesolverreport* rep, 3278 /* Real */ ae_matrix* x, 3279 ae_state *_state); 3280 void cmatrixsolvem(/* Complex */ ae_matrix* a, 3281 ae_int_t n, 3282 /* Complex */ ae_matrix* b, 3283 ae_int_t m, 3284 ae_bool rfs, 3285 ae_int_t* info, 3286 densesolverreport* rep, 3287 /* Complex */ ae_matrix* x, 3288 ae_state *_state); 3289 void cmatrixsolvemfast(/* Complex */ ae_matrix* a, 3290 ae_int_t n, 3291 /* Complex */ ae_matrix* b, 3292 ae_int_t m, 3293 ae_int_t* info, 3294 ae_state *_state); 3295 void cmatrixsolve(/* Complex */ ae_matrix* a, 3296 ae_int_t n, 3297 /* Complex */ ae_vector* b, 3298 ae_int_t* info, 3299 densesolverreport* rep, 3300 /* Complex */ ae_vector* x, 3301 ae_state *_state); 3302 void cmatrixsolvefast(/* Complex */ ae_matrix* a, 3303 ae_int_t n, 3304 /* Complex */ ae_vector* b, 3305 ae_int_t* info, 3306 ae_state *_state); 3307 void cmatrixlusolvem(/* Complex */ ae_matrix* lua, 3308 /* Integer */ ae_vector* p, 3309 ae_int_t n, 3310 /* Complex */ ae_matrix* b, 3311 ae_int_t m, 3312 ae_int_t* info, 3313 densesolverreport* rep, 3314 /* Complex */ ae_matrix* x, 3315 ae_state *_state); 3316 void cmatrixlusolvemfast(/* Complex */ ae_matrix* lua, 3317 /* Integer */ ae_vector* p, 3318 ae_int_t n, 3319 /* Complex */ ae_matrix* b, 3320 ae_int_t m, 3321 ae_int_t* info, 3322 ae_state *_state); 3323 void cmatrixlusolve(/* Complex */ ae_matrix* lua, 3324 /* Integer */ ae_vector* p, 3325 ae_int_t n, 3326 /* Complex */ ae_vector* b, 3327 ae_int_t* info, 3328 densesolverreport* rep, 3329 /* Complex */ ae_vector* x, 3330 ae_state *_state); 3331 void cmatrixlusolvefast(/* Complex */ ae_matrix* lua, 3332 /* Integer */ ae_vector* p, 3333 ae_int_t n, 3334 /* Complex */ ae_vector* b, 3335 ae_int_t* info, 3336 ae_state *_state); 3337 void cmatrixmixedsolvem(/* Complex */ ae_matrix* a, 3338 /* Complex */ ae_matrix* lua, 3339 /* Integer */ ae_vector* p, 3340 ae_int_t n, 3341 /* Complex */ ae_matrix* b, 3342 ae_int_t m, 3343 ae_int_t* info, 3344 densesolverreport* rep, 3345 /* Complex */ ae_matrix* x, 3346 ae_state *_state); 3347 void cmatrixmixedsolve(/* Complex */ ae_matrix* a, 3348 /* Complex */ ae_matrix* lua, 3349 /* Integer */ ae_vector* p, 3350 ae_int_t n, 3351 /* Complex */ ae_vector* b, 3352 ae_int_t* info, 3353 densesolverreport* rep, 3354 /* Complex */ ae_vector* x, 3355 ae_state *_state); 3356 void spdmatrixsolvem(/* Real */ ae_matrix* a, 3357 ae_int_t n, 3358 ae_bool isupper, 3359 /* Real */ ae_matrix* b, 3360 ae_int_t m, 3361 ae_int_t* info, 3362 densesolverreport* rep, 3363 /* Real */ ae_matrix* x, 3364 ae_state *_state); 3365 void spdmatrixsolvemfast(/* Real */ ae_matrix* a, 3366 ae_int_t n, 3367 ae_bool isupper, 3368 /* Real */ ae_matrix* b, 3369 ae_int_t m, 3370 ae_int_t* info, 3371 ae_state *_state); 3372 void spdmatrixsolve(/* Real */ ae_matrix* a, 3373 ae_int_t n, 3374 ae_bool isupper, 3375 /* Real */ ae_vector* b, 3376 ae_int_t* info, 3377 densesolverreport* rep, 3378 /* Real */ ae_vector* x, 3379 ae_state *_state); 3380 void spdmatrixsolvefast(/* Real */ ae_matrix* a, 3381 ae_int_t n, 3382 ae_bool isupper, 3383 /* Real */ ae_vector* b, 3384 ae_int_t* info, 3385 ae_state *_state); 3386 void spdmatrixcholeskysolvem(/* Real */ ae_matrix* cha, 3387 ae_int_t n, 3388 ae_bool isupper, 3389 /* Real */ ae_matrix* b, 3390 ae_int_t m, 3391 ae_int_t* info, 3392 densesolverreport* rep, 3393 /* Real */ ae_matrix* x, 3394 ae_state *_state); 3395 void spdmatrixcholeskysolvemfast(/* Real */ ae_matrix* cha, 3396 ae_int_t n, 3397 ae_bool isupper, 3398 /* Real */ ae_matrix* b, 3399 ae_int_t m, 3400 ae_int_t* info, 3401 ae_state *_state); 3402 void spdmatrixcholeskysolve(/* Real */ ae_matrix* cha, 3403 ae_int_t n, 3404 ae_bool isupper, 3405 /* Real */ ae_vector* b, 3406 ae_int_t* info, 3407 densesolverreport* rep, 3408 /* Real */ ae_vector* x, 3409 ae_state *_state); 3410 void spdmatrixcholeskysolvefast(/* Real */ ae_matrix* cha, 3411 ae_int_t n, 3412 ae_bool isupper, 3413 /* Real */ ae_vector* b, 3414 ae_int_t* info, 3415 ae_state *_state); 3416 void hpdmatrixsolvem(/* Complex */ ae_matrix* a, 3417 ae_int_t n, 3418 ae_bool isupper, 3419 /* Complex */ ae_matrix* b, 3420 ae_int_t m, 3421 ae_int_t* info, 3422 densesolverreport* rep, 3423 /* Complex */ ae_matrix* x, 3424 ae_state *_state); 3425 void hpdmatrixsolvemfast(/* Complex */ ae_matrix* a, 3426 ae_int_t n, 3427 ae_bool isupper, 3428 /* Complex */ ae_matrix* b, 3429 ae_int_t m, 3430 ae_int_t* info, 3431 ae_state *_state); 3432 void hpdmatrixsolve(/* Complex */ ae_matrix* a, 3433 ae_int_t n, 3434 ae_bool isupper, 3435 /* Complex */ ae_vector* b, 3436 ae_int_t* info, 3437 densesolverreport* rep, 3438 /* Complex */ ae_vector* x, 3439 ae_state *_state); 3440 void hpdmatrixsolvefast(/* Complex */ ae_matrix* a, 3441 ae_int_t n, 3442 ae_bool isupper, 3443 /* Complex */ ae_vector* b, 3444 ae_int_t* info, 3445 ae_state *_state); 3446 void hpdmatrixcholeskysolvem(/* Complex */ ae_matrix* cha, 3447 ae_int_t n, 3448 ae_bool isupper, 3449 /* Complex */ ae_matrix* b, 3450 ae_int_t m, 3451 ae_int_t* info, 3452 densesolverreport* rep, 3453 /* Complex */ ae_matrix* x, 3454 ae_state *_state); 3455 void hpdmatrixcholeskysolvemfast(/* Complex */ ae_matrix* cha, 3456 ae_int_t n, 3457 ae_bool isupper, 3458 /* Complex */ ae_matrix* b, 3459 ae_int_t m, 3460 ae_int_t* info, 3461 ae_state *_state); 3462 void hpdmatrixcholeskysolve(/* Complex */ ae_matrix* cha, 3463 ae_int_t n, 3464 ae_bool isupper, 3465 /* Complex */ ae_vector* b, 3466 ae_int_t* info, 3467 densesolverreport* rep, 3468 /* Complex */ ae_vector* x, 3469 ae_state *_state); 3470 void hpdmatrixcholeskysolvefast(/* Complex */ ae_matrix* cha, 3471 ae_int_t n, 3472 ae_bool isupper, 3473 /* Complex */ ae_vector* b, 3474 ae_int_t* info, 3475 ae_state *_state); 3476 void rmatrixsolvels(/* Real */ ae_matrix* a, 3477 ae_int_t nrows, 3478 ae_int_t ncols, 3479 /* Real */ ae_vector* b, 3480 double threshold, 3481 ae_int_t* info, 3482 densesolverlsreport* rep, 3483 /* Real */ ae_vector* x, 3484 ae_state *_state); 3485 void _densesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3486 void _densesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3487 void _densesolverreport_clear(void* _p); 3488 void _densesolverreport_destroy(void* _p); 3489 void _densesolverlsreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3490 void _densesolverlsreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3491 void _densesolverlsreport_clear(void* _p); 3492 void _densesolverlsreport_destroy(void* _p); 3493 #endif 3494 #if defined(AE_COMPILE_LINLSQR) || !defined(AE_PARTIAL_BUILD) 3495 void linlsqrcreate(ae_int_t m, 3496 ae_int_t n, 3497 linlsqrstate* state, 3498 ae_state *_state); 3499 void linlsqrcreatebuf(ae_int_t m, 3500 ae_int_t n, 3501 linlsqrstate* state, 3502 ae_state *_state); 3503 void linlsqrsetb(linlsqrstate* state, 3504 /* Real */ ae_vector* b, 3505 ae_state *_state); 3506 void linlsqrsetprecunit(linlsqrstate* state, ae_state *_state); 3507 void linlsqrsetprecdiag(linlsqrstate* state, ae_state *_state); 3508 void linlsqrsetlambdai(linlsqrstate* state, 3509 double lambdai, 3510 ae_state *_state); 3511 ae_bool linlsqriteration(linlsqrstate* state, ae_state *_state); 3512 void linlsqrsolvesparse(linlsqrstate* state, 3513 sparsematrix* a, 3514 /* Real */ ae_vector* b, 3515 ae_state *_state); 3516 void linlsqrsetcond(linlsqrstate* state, 3517 double epsa, 3518 double epsb, 3519 ae_int_t maxits, 3520 ae_state *_state); 3521 void linlsqrresults(linlsqrstate* state, 3522 /* Real */ ae_vector* x, 3523 linlsqrreport* rep, 3524 ae_state *_state); 3525 void linlsqrsetxrep(linlsqrstate* state, 3526 ae_bool needxrep, 3527 ae_state *_state); 3528 void linlsqrrestart(linlsqrstate* state, ae_state *_state); 3529 ae_int_t linlsqrpeekiterationscount(linlsqrstate* s, ae_state *_state); 3530 void linlsqrrequesttermination(linlsqrstate* state, ae_state *_state); 3531 void _linlsqrstate_init(void* _p, ae_state *_state, ae_bool make_automatic); 3532 void _linlsqrstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3533 void _linlsqrstate_clear(void* _p); 3534 void _linlsqrstate_destroy(void* _p); 3535 void _linlsqrreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3536 void _linlsqrreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3537 void _linlsqrreport_clear(void* _p); 3538 void _linlsqrreport_destroy(void* _p); 3539 #endif 3540 #if defined(AE_COMPILE_POLYNOMIALSOLVER) || !defined(AE_PARTIAL_BUILD) 3541 void polynomialsolve(/* Real */ ae_vector* a, 3542 ae_int_t n, 3543 /* Complex */ ae_vector* x, 3544 polynomialsolverreport* rep, 3545 ae_state *_state); 3546 void _polynomialsolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3547 void _polynomialsolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3548 void _polynomialsolverreport_clear(void* _p); 3549 void _polynomialsolverreport_destroy(void* _p); 3550 #endif 3551 #if defined(AE_COMPILE_NLEQ) || !defined(AE_PARTIAL_BUILD) 3552 void nleqcreatelm(ae_int_t n, 3553 ae_int_t m, 3554 /* Real */ ae_vector* x, 3555 nleqstate* state, 3556 ae_state *_state); 3557 void nleqsetcond(nleqstate* state, 3558 double epsf, 3559 ae_int_t maxits, 3560 ae_state *_state); 3561 void nleqsetxrep(nleqstate* state, ae_bool needxrep, ae_state *_state); 3562 void nleqsetstpmax(nleqstate* state, double stpmax, ae_state *_state); 3563 ae_bool nleqiteration(nleqstate* state, ae_state *_state); 3564 void nleqresults(nleqstate* state, 3565 /* Real */ ae_vector* x, 3566 nleqreport* rep, 3567 ae_state *_state); 3568 void nleqresultsbuf(nleqstate* state, 3569 /* Real */ ae_vector* x, 3570 nleqreport* rep, 3571 ae_state *_state); 3572 void nleqrestartfrom(nleqstate* state, 3573 /* Real */ ae_vector* x, 3574 ae_state *_state); 3575 void _nleqstate_init(void* _p, ae_state *_state, ae_bool make_automatic); 3576 void _nleqstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3577 void _nleqstate_clear(void* _p); 3578 void _nleqstate_destroy(void* _p); 3579 void _nleqreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3580 void _nleqreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3581 void _nleqreport_clear(void* _p); 3582 void _nleqreport_destroy(void* _p); 3583 #endif 3584 #if defined(AE_COMPILE_DIRECTSPARSESOLVERS) || !defined(AE_PARTIAL_BUILD) 3585 void sparsesolvesks(sparsematrix* a, 3586 ae_int_t n, 3587 ae_bool isupper, 3588 /* Real */ ae_vector* b, 3589 sparsesolverreport* rep, 3590 /* Real */ ae_vector* x, 3591 ae_state *_state); 3592 void sparsecholeskysolvesks(sparsematrix* a, 3593 ae_int_t n, 3594 ae_bool isupper, 3595 /* Real */ ae_vector* b, 3596 sparsesolverreport* rep, 3597 /* Real */ ae_vector* x, 3598 ae_state *_state); 3599 void _sparsesolverreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3600 void _sparsesolverreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3601 void _sparsesolverreport_clear(void* _p); 3602 void _sparsesolverreport_destroy(void* _p); 3603 #endif 3604 #if defined(AE_COMPILE_LINCG) || !defined(AE_PARTIAL_BUILD) 3605 void lincgcreate(ae_int_t n, lincgstate* state, ae_state *_state); 3606 void lincgsetstartingpoint(lincgstate* state, 3607 /* Real */ ae_vector* x, 3608 ae_state *_state); 3609 void lincgsetb(lincgstate* state, 3610 /* Real */ ae_vector* b, 3611 ae_state *_state); 3612 void lincgsetprecunit(lincgstate* state, ae_state *_state); 3613 void lincgsetprecdiag(lincgstate* state, ae_state *_state); 3614 void lincgsetcond(lincgstate* state, 3615 double epsf, 3616 ae_int_t maxits, 3617 ae_state *_state); 3618 ae_bool lincgiteration(lincgstate* state, ae_state *_state); 3619 void lincgsolvesparse(lincgstate* state, 3620 sparsematrix* a, 3621 ae_bool isupper, 3622 /* Real */ ae_vector* b, 3623 ae_state *_state); 3624 void lincgresults(lincgstate* state, 3625 /* Real */ ae_vector* x, 3626 lincgreport* rep, 3627 ae_state *_state); 3628 void lincgsetrestartfreq(lincgstate* state, 3629 ae_int_t srf, 3630 ae_state *_state); 3631 void lincgsetrupdatefreq(lincgstate* state, 3632 ae_int_t freq, 3633 ae_state *_state); 3634 void lincgsetxrep(lincgstate* state, ae_bool needxrep, ae_state *_state); 3635 void lincgrestart(lincgstate* state, ae_state *_state); 3636 void _lincgstate_init(void* _p, ae_state *_state, ae_bool make_automatic); 3637 void _lincgstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3638 void _lincgstate_clear(void* _p); 3639 void _lincgstate_destroy(void* _p); 3640 void _lincgreport_init(void* _p, ae_state *_state, ae_bool make_automatic); 3641 void _lincgreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic); 3642 void _lincgreport_clear(void* _p); 3643 void _lincgreport_destroy(void* _p); 3644 #endif 3645 3646 } 3647 #endif 3648 3649