Context Navigation

source: branches/GeneralizedQAP/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.1.0/ALGLIB-3.1.0/solvers.cs @ 6843

Visit:

Last change on this file since 6843 was 4977, checked in by gkronber, 14 years ago
Added new alglib classes (version 3.1.0) #875
File size: 169.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences