Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/densesolver.cs @ 3839

Visit:

Last change on this file since 3839 was 3839, checked in by mkommend, 14 years ago
implemented first version of LR (ticket #1012)
File size: 94.8 KB

Line
1	/*************************************************************************
2	Copyright (c) 2007-2008, 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
18	>>> END OF LICENSE >>>
19	*************************************************************************/
20
21	using System;
22
23	namespace alglib
24	{
25	public class densesolver
26	{
27	public struct densesolverreport
28	{
29	public double r1;
30	public double rinf;
31	};
32
33
34	public struct densesolverlsreport
35	{
36	public double r2;
37	public double[,] cx;
38	public int n;
39	public int k;
40	};
41
42
43
44
45	/*************************************************************************
46	Dense solver.
47
48	This subroutine solves a system A*x=b, where A is NxN non-denegerate
49	real matrix, x and b are vectors.
50
51	Algorithm features:
52	* automatic detection of degenerate cases
53	* condition number estimation
54	* iterative refinement
55	* O(N^3) complexity
56
57	INPUT PARAMETERS
58	A - array[0..N-1,0..N-1], system matrix
59	N - size of A
60	B - array[0..N-1], right part
61
62	OUTPUT PARAMETERS
63	Info - return code:
64	* -3 A is singular, or VERY close to singular.
65	X is filled by zeros in such cases.
66	* -1 N<=0 was passed
67	* 1 task is solved (but matrix A may be ill-conditioned,
68	check R1/RInf parameters for condition numbers).
69	Rep - solver report, see below for more info
70	X - array[0..N-1], it contains:
71	* solution of A*x=b if A is non-singular (well-conditioned
72	or ill-conditioned, but not very close to singular)
73	* zeros, if A is singular or VERY close to singular
74	(in this case Info=-3).
75
76	SOLVER REPORT
77
78	Subroutine sets following fields of the Rep structure:
79	* R1 reciprocal of condition number: 1/cond(A), 1-norm.
80	* RInf reciprocal of condition number: 1/cond(A), inf-norm.
81
82	-- ALGLIB --
83	Copyright 27.01.2010 by Bochkanov Sergey
84	*************************************************************************/
85	public static void rmatrixsolve(ref double[,] a,
86	int n,
87	ref double[] b,
88	ref int info,
89	ref densesolverreport rep,
90	ref double[] x)
91	{
92	double[,] bm = new double[0,0];
93	double[,] xm = new double[0,0];
94	int i_ = 0;
95
96	if( n<=0 )
97	{
98	info = -1;
99	return;
100	}
101	bm = new double[n, 1];
102	for(i_=0; i_<=n-1;i_++)
103	{
104	bm[i_,0] = b[i_];
105	}
106	rmatrixsolvem(ref a, n, ref bm, 1, true, ref info, ref rep, ref xm);
107	x = new double[n];
108	for(i_=0; i_<=n-1;i_++)
109	{
110	x[i_] = xm[i_,0];
111	}
112	}
113
114
115	/*************************************************************************
116	Dense solver.
117
118	Similar to RMatrixSolve() but solves task with multiple right parts (where
119	b and x are NxM matrices).
120
121	Algorithm features:
122	* automatic detection of degenerate cases
123	* condition number estimation
124	* optional iterative refinement
125	* O(N^3+M*N^2) complexity
126
127	INPUT PARAMETERS
128	A - array[0..N-1,0..N-1], system matrix
129	N - size of A
130	B - array[0..N-1,0..M-1], right part
131	M - right part size
132	RFS - iterative refinement switch:
133	* True - refinement is used.
134	Less performance, more precision.
135	* False - refinement is not used.
136	More performance, less precision.
137
138	OUTPUT PARAMETERS
139	Info - same as in RMatrixSolve
140	Rep - same as in RMatrixSolve
141	X - same as in RMatrixSolve
142
143	-- ALGLIB --
144	Copyright 27.01.2010 by Bochkanov Sergey
145	*************************************************************************/
146	public static void rmatrixsolvem(ref double[,] a,
147	int n,
148	ref double[,] b,
149	int m,
150	bool rfs,
151	ref int info,
152	ref densesolverreport rep,
153	ref double[,] x)
154	{
155	double[,] da = new double[0,0];
156	double[,] emptya = new double[0,0];
157	int[] p = new int[0];
158	double scalea = 0;
159	int i = 0;
160	int j = 0;
161	int i_ = 0;
162
163
164	//
165	// prepare: check inputs, allocate space...
166	//
167	if( n<=0 \| m<=0 )
168	{
169	info = -1;
170	return;
171	}
172	da = new double[n, n];
173
174	//
175	// 1. scale matrix, max(\|A[i,j]\|)
176	// 2. factorize scaled matrix
177	// 3. solve
178	//
179	scalea = 0;
180	for(i=0; i<=n-1; i++)
181	{
182	for(j=0; j<=n-1; j++)
183	{
184	scalea = Math.Max(scalea, Math.Abs(a[i,j]));
185	}
186	}
187	if( (double)(scalea)==(double)(0) )
188	{
189	scalea = 1;
190	}
191	scalea = 1/scalea;
192	for(i=0; i<=n-1; i++)
193	{
194	for(i_=0; i_<=n-1;i_++)
195	{
196	da[i,i_] = a[i,i_];
197	}
198	}
199	trfac.rmatrixlu(ref da, n, n, ref p);
200	if( rfs )
201	{
202	rmatrixlusolveinternal(ref da, ref p, scalea, n, ref a, true, ref b, m, ref info, ref rep, ref x);
203	}
204	else
205	{
206	rmatrixlusolveinternal(ref da, ref p, scalea, n, ref emptya, false, ref b, m, ref info, ref rep, ref x);
207	}
208	}
209
210
211	/*************************************************************************
212	Dense solver.
213
214	This subroutine solves a system A*X=B, where A is NxN non-denegerate
215	real matrix given by its LU decomposition, X and B are NxM real matrices.
216
217	Algorithm features:
218	* automatic detection of degenerate cases
219	* O(N^2) complexity
220	* condition number estimation
221
222	No iterative refinement is provided because exact form of original matrix
223	is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
224	need iterative refinement.
225
226	INPUT PARAMETERS
227	LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
228	P - array[0..N-1], pivots array, RMatrixLU result
229	N - size of A
230	B - array[0..N-1], right part
231
232	OUTPUT PARAMETERS
233	Info - same as in RMatrixSolve
234	Rep - same as in RMatrixSolve
235	X - same as in RMatrixSolve
236
237	-- ALGLIB --
238	Copyright 27.01.2010 by Bochkanov Sergey
239	*************************************************************************/
240	public static void rmatrixlusolve(ref double[,] lua,
241	ref int[] p,
242	int n,
243	ref double[] b,
244	ref int info,
245	ref densesolverreport rep,
246	ref double[] x)
247	{
248	double[,] bm = new double[0,0];
249	double[,] xm = new double[0,0];
250	int i_ = 0;
251
252	if( n<=0 )
253	{
254	info = -1;
255	return;
256	}
257	bm = new double[n, 1];
258	for(i_=0; i_<=n-1;i_++)
259	{
260	bm[i_,0] = b[i_];
261	}
262	rmatrixlusolvem(ref lua, ref p, n, ref bm, 1, ref info, ref rep, ref xm);
263	x = new double[n];
264	for(i_=0; i_<=n-1;i_++)
265	{
266	x[i_] = xm[i_,0];
267	}
268	}
269
270
271	/*************************************************************************
272	Dense solver.
273
274	Similar to RMatrixLUSolve() but solves task with multiple right parts
275	(where b and x are NxM matrices).
276
277	Algorithm features:
278	* automatic detection of degenerate cases
279	* O(M*N^2) complexity
280	* condition number estimation
281
282	No iterative refinement is provided because exact form of original matrix
283	is not known to subroutine. Use RMatrixSolve or RMatrixMixedSolve if you
284	need iterative refinement.
285
286	INPUT PARAMETERS
287	LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
288	P - array[0..N-1], pivots array, RMatrixLU result
289	N - size of A
290	B - array[0..N-1,0..M-1], right part
291	M - right part size
292
293	OUTPUT PARAMETERS
294	Info - same as in RMatrixSolve
295	Rep - same as in RMatrixSolve
296	X - same as in RMatrixSolve
297
298	-- ALGLIB --
299	Copyright 27.01.2010 by Bochkanov Sergey
300	*************************************************************************/
301	public static void rmatrixlusolvem(ref double[,] lua,
302	ref int[] p,
303	int n,
304	ref double[,] b,
305	int m,
306	ref int info,
307	ref densesolverreport rep,
308	ref double[,] x)
309	{
310	double[,] emptya = new double[0,0];
311	int i = 0;
312	int j = 0;
313	double scalea = 0;
314
315
316	//
317	// prepare: check inputs, allocate space...
318	//
319	if( n<=0 \| m<=0 )
320	{
321	info = -1;
322	return;
323	}
324
325	//
326	// 1. scale matrix, max(\|U[i,j]\|)
327	// we assume that LU is in its normal form, i.e. \|L[i,j]\|<=1
328	// 2. solve
329	//
330	scalea = 0;
331	for(i=0; i<=n-1; i++)
332	{
333	for(j=i; j<=n-1; j++)
334	{
335	scalea = Math.Max(scalea, Math.Abs(lua[i,j]));
336	}
337	}
338	if( (double)(scalea)==(double)(0) )
339	{
340	scalea = 1;
341	}
342	scalea = 1/scalea;
343	rmatrixlusolveinternal(ref lua, ref p, scalea, n, ref emptya, false, ref b, m, ref info, ref rep, ref x);
344	}
345
346
347	/*************************************************************************
348	Dense solver.
349
350	This subroutine solves a system A*x=b, where BOTH ORIGINAL A AND ITS
351	LU DECOMPOSITION ARE KNOWN. You can use it if for some reasons you have
352	both A and its LU decomposition.
353
354	Algorithm features:
355	* automatic detection of degenerate cases
356	* condition number estimation
357	* iterative refinement
358	* O(N^2) complexity
359
360	INPUT PARAMETERS
361	A - array[0..N-1,0..N-1], system matrix
362	LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
363	P - array[0..N-1], pivots array, RMatrixLU result
364	N - size of A
365	B - array[0..N-1], right part
366
367	OUTPUT PARAMETERS
368	Info - same as in RMatrixSolveM
369	Rep - same as in RMatrixSolveM
370	X - same as in RMatrixSolveM
371
372	-- ALGLIB --
373	Copyright 27.01.2010 by Bochkanov Sergey
374	*************************************************************************/
375	public static void rmatrixmixedsolve(ref double[,] a,
376	ref double[,] lua,
377	ref int[] p,
378	int n,
379	ref double[] b,
380	ref int info,
381	ref densesolverreport rep,
382	ref double[] x)
383	{
384	double[,] bm = new double[0,0];
385	double[,] xm = new double[0,0];
386	int i_ = 0;
387
388	if( n<=0 )
389	{
390	info = -1;
391	return;
392	}
393	bm = new double[n, 1];
394	for(i_=0; i_<=n-1;i_++)
395	{
396	bm[i_,0] = b[i_];
397	}
398	rmatrixmixedsolvem(ref a, ref lua, ref p, n, ref bm, 1, ref info, ref rep, ref xm);
399	x = new double[n];
400	for(i_=0; i_<=n-1;i_++)
401	{
402	x[i_] = xm[i_,0];
403	}
404	}
405
406
407	/*************************************************************************
408	Dense solver.
409
410	Similar to RMatrixMixedSolve() but solves task with multiple right parts
411	(where b and x are NxM matrices).
412
413	Algorithm features:
414	* automatic detection of degenerate cases
415	* condition number estimation
416	* iterative refinement
417	* O(M*N^2) complexity
418
419	INPUT PARAMETERS
420	A - array[0..N-1,0..N-1], system matrix
421	LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
422	P - array[0..N-1], pivots array, RMatrixLU result
423	N - size of A
424	B - array[0..N-1,0..M-1], right part
425	M - right part size
426
427	OUTPUT PARAMETERS
428	Info - same as in RMatrixSolveM
429	Rep - same as in RMatrixSolveM
430	X - same as in RMatrixSolveM
431
432	-- ALGLIB --
433	Copyright 27.01.2010 by Bochkanov Sergey
434	*************************************************************************/
435	public static void rmatrixmixedsolvem(ref double[,] a,
436	ref double[,] lua,
437	ref int[] p,
438	int n,
439	ref double[,] b,
440	int m,
441	ref int info,
442	ref densesolverreport rep,
443	ref double[,] x)
444	{
445	double scalea = 0;
446	int i = 0;
447	int j = 0;
448
449
450	//
451	// prepare: check inputs, allocate space...
452	//
453	if( n<=0 \| m<=0 )
454	{
455	info = -1;
456	return;
457	}
458
459	//
460	// 1. scale matrix, max(\|A[i,j]\|)
461	// 2. factorize scaled matrix
462	// 3. solve
463	//
464	scalea = 0;
465	for(i=0; i<=n-1; i++)
466	{
467	for(j=0; j<=n-1; j++)
468	{
469	scalea = Math.Max(scalea, Math.Abs(a[i,j]));
470	}
471	}
472	if( (double)(scalea)==(double)(0) )
473	{
474	scalea = 1;
475	}
476	scalea = 1/scalea;
477	rmatrixlusolveinternal(ref lua, ref p, scalea, n, ref a, true, ref b, m, ref info, ref rep, ref x);
478	}
479
480
481	/*************************************************************************
482	Dense solver. Same as RMatrixSolveM(), but for complex matrices.
483
484	Algorithm features:
485	* automatic detection of degenerate cases
486	* condition number estimation
487	* iterative refinement
488	* O(N^3+M*N^2) complexity
489
490	INPUT PARAMETERS
491	A - array[0..N-1,0..N-1], system matrix
492	N - size of A
493	B - array[0..N-1,0..M-1], right part
494	M - right part size
495	RFS - iterative refinement switch:
496	* True - refinement is used.
497	Less performance, more precision.
498	* False - refinement is not used.
499	More performance, less precision.
500
501	OUTPUT PARAMETERS
502	Info - same as in RMatrixSolve
503	Rep - same as in RMatrixSolve
504	X - same as in RMatrixSolve
505
506	-- ALGLIB --
507	Copyright 27.01.2010 by Bochkanov Sergey
508	*************************************************************************/
509	public static void cmatrixsolvem(ref AP.Complex[,] a,
510	int n,
511	ref AP.Complex[,] b,
512	int m,
513	bool rfs,
514	ref int info,
515	ref densesolverreport rep,
516	ref AP.Complex[,] x)
517	{
518	AP.Complex[,] da = new AP.Complex[0,0];
519	AP.Complex[,] emptya = new AP.Complex[0,0];
520	int[] p = new int[0];
521	double scalea = 0;
522	int i = 0;
523	int j = 0;
524	int i_ = 0;
525
526
527	//
528	// prepare: check inputs, allocate space...
529	//
530	if( n<=0 \| m<=0 )
531	{
532	info = -1;
533	return;
534	}
535	da = new AP.Complex[n, n];
536
537	//
538	// 1. scale matrix, max(\|A[i,j]\|)
539	// 2. factorize scaled matrix
540	// 3. solve
541	//
542	scalea = 0;
543	for(i=0; i<=n-1; i++)
544	{
545	for(j=0; j<=n-1; j++)
546	{
547	scalea = Math.Max(scalea, AP.Math.AbsComplex(a[i,j]));
548	}
549	}
550	if( (double)(scalea)==(double)(0) )
551	{
552	scalea = 1;
553	}
554	scalea = 1/scalea;
555	for(i=0; i<=n-1; i++)
556	{
557	for(i_=0; i_<=n-1;i_++)
558	{
559	da[i,i_] = a[i,i_];
560	}
561	}
562	trfac.cmatrixlu(ref da, n, n, ref p);
563	if( rfs )
564	{
565	cmatrixlusolveinternal(ref da, ref p, scalea, n, ref a, true, ref b, m, ref info, ref rep, ref x);
566	}
567	else
568	{
569	cmatrixlusolveinternal(ref da, ref p, scalea, n, ref emptya, false, ref b, m, ref info, ref rep, ref x);
570	}
571	}
572
573
574	/*************************************************************************
575	Dense solver. Same as RMatrixSolve(), but for complex matrices.
576
577	Algorithm features:
578	* automatic detection of degenerate cases
579	* condition number estimation
580	* iterative refinement
581	* O(N^3) complexity
582
583	INPUT PARAMETERS
584	A - array[0..N-1,0..N-1], system matrix
585	N - size of A
586	B - array[0..N-1], right part
587
588	OUTPUT PARAMETERS
589	Info - same as in RMatrixSolve
590	Rep - same as in RMatrixSolve
591	X - same as in RMatrixSolve
592
593	-- ALGLIB --
594	Copyright 27.01.2010 by Bochkanov Sergey
595	*************************************************************************/
596	public static void cmatrixsolve(ref AP.Complex[,] a,
597	int n,
598	ref AP.Complex[] b,
599	ref int info,
600	ref densesolverreport rep,
601	ref AP.Complex[] x)
602	{
603	AP.Complex[,] bm = new AP.Complex[0,0];
604	AP.Complex[,] xm = new AP.Complex[0,0];
605	int i_ = 0;
606
607	if( n<=0 )
608	{
609	info = -1;
610	return;
611	}
612	bm = new AP.Complex[n, 1];
613	for(i_=0; i_<=n-1;i_++)
614	{
615	bm[i_,0] = b[i_];
616	}
617	cmatrixsolvem(ref a, n, ref bm, 1, true, ref info, ref rep, ref xm);
618	x = new AP.Complex[n];
619	for(i_=0; i_<=n-1;i_++)
620	{
621	x[i_] = xm[i_,0];
622	}
623	}
624
625
626	/*************************************************************************
627	Dense solver. Same as RMatrixLUSolveM(), but for complex matrices.
628
629	Algorithm features:
630	* automatic detection of degenerate cases
631	* O(M*N^2) complexity
632	* condition number estimation
633
634	No iterative refinement is provided because exact form of original matrix
635	is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
636	need iterative refinement.
637
638	INPUT PARAMETERS
639	LUA - array[0..N-1,0..N-1], LU decomposition, RMatrixLU result
640	P - array[0..N-1], pivots array, RMatrixLU result
641	N - size of A
642	B - array[0..N-1,0..M-1], right part
643	M - right part size
644
645	OUTPUT PARAMETERS
646	Info - same as in RMatrixSolve
647	Rep - same as in RMatrixSolve
648	X - same as in RMatrixSolve
649
650	-- ALGLIB --
651	Copyright 27.01.2010 by Bochkanov Sergey
652	*************************************************************************/
653	public static void cmatrixlusolvem(ref AP.Complex[,] lua,
654	ref int[] p,
655	int n,
656	ref AP.Complex[,] b,
657	int m,
658	ref int info,
659	ref densesolverreport rep,
660	ref AP.Complex[,] x)
661	{
662	AP.Complex[,] emptya = new AP.Complex[0,0];
663	int i = 0;
664	int j = 0;
665	double scalea = 0;
666
667
668	//
669	// prepare: check inputs, allocate space...
670	//
671	if( n<=0 \| m<=0 )
672	{
673	info = -1;
674	return;
675	}
676
677	//
678	// 1. scale matrix, max(\|U[i,j]\|)
679	// we assume that LU is in its normal form, i.e. \|L[i,j]\|<=1
680	// 2. solve
681	//
682	scalea = 0;
683	for(i=0; i<=n-1; i++)
684	{
685	for(j=i; j<=n-1; j++)
686	{
687	scalea = Math.Max(scalea, AP.Math.AbsComplex(lua[i,j]));
688	}
689	}
690	if( (double)(scalea)==(double)(0) )
691	{
692	scalea = 1;
693	}
694	scalea = 1/scalea;
695	cmatrixlusolveinternal(ref lua, ref p, scalea, n, ref emptya, false, ref b, m, ref info, ref rep, ref x);
696	}
697
698
699	/*************************************************************************
700	Dense solver. Same as RMatrixLUSolve(), but for complex matrices.
701
702	Algorithm features:
703	* automatic detection of degenerate cases
704	* O(N^2) complexity
705	* condition number estimation
706
707	No iterative refinement is provided because exact form of original matrix
708	is not known to subroutine. Use CMatrixSolve or CMatrixMixedSolve if you
709	need iterative refinement.
710
711	INPUT PARAMETERS
712	LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
713	P - array[0..N-1], pivots array, CMatrixLU result
714	N - size of A
715	B - array[0..N-1], right part
716
717	OUTPUT PARAMETERS
718	Info - same as in RMatrixSolve
719	Rep - same as in RMatrixSolve
720	X - same as in RMatrixSolve
721
722	-- ALGLIB --
723	Copyright 27.01.2010 by Bochkanov Sergey
724	*************************************************************************/
725	public static void cmatrixlusolve(ref AP.Complex[,] lua,
726	ref int[] p,
727	int n,
728	ref AP.Complex[] b,
729	ref int info,
730	ref densesolverreport rep,
731	ref AP.Complex[] x)
732	{
733	AP.Complex[,] bm = new AP.Complex[0,0];
734	AP.Complex[,] xm = new AP.Complex[0,0];
735	int i_ = 0;
736
737	if( n<=0 )
738	{
739	info = -1;
740	return;
741	}
742	bm = new AP.Complex[n, 1];
743	for(i_=0; i_<=n-1;i_++)
744	{
745	bm[i_,0] = b[i_];
746	}
747	cmatrixlusolvem(ref lua, ref p, n, ref bm, 1, ref info, ref rep, ref xm);
748	x = new AP.Complex[n];
749	for(i_=0; i_<=n-1;i_++)
750	{
751	x[i_] = xm[i_,0];
752	}
753	}
754
755
756	/*************************************************************************
757	Dense solver. Same as RMatrixMixedSolveM(), but for complex matrices.
758
759	Algorithm features:
760	* automatic detection of degenerate cases
761	* condition number estimation
762	* iterative refinement
763	* O(M*N^2) complexity
764
765	INPUT PARAMETERS
766	A - array[0..N-1,0..N-1], system matrix
767	LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
768	P - array[0..N-1], pivots array, CMatrixLU result
769	N - size of A
770	B - array[0..N-1,0..M-1], right part
771	M - right part size
772
773	OUTPUT PARAMETERS
774	Info - same as in RMatrixSolveM
775	Rep - same as in RMatrixSolveM
776	X - same as in RMatrixSolveM
777
778	-- ALGLIB --
779	Copyright 27.01.2010 by Bochkanov Sergey
780	*************************************************************************/
781	public static void cmatrixmixedsolvem(ref AP.Complex[,] a,
782	ref AP.Complex[,] lua,
783	ref int[] p,
784	int n,
785	ref AP.Complex[,] b,
786	int m,
787	ref int info,
788	ref densesolverreport rep,
789	ref AP.Complex[,] x)
790	{
791	double scalea = 0;
792	int i = 0;
793	int j = 0;
794
795
796	//
797	// prepare: check inputs, allocate space...
798	//
799	if( n<=0 \| m<=0 )
800	{
801	info = -1;
802	return;
803	}
804
805	//
806	// 1. scale matrix, max(\|A[i,j]\|)
807	// 2. factorize scaled matrix
808	// 3. solve
809	//
810	scalea = 0;
811	for(i=0; i<=n-1; i++)
812	{
813	for(j=0; j<=n-1; j++)
814	{
815	scalea = Math.Max(scalea, AP.Math.AbsComplex(a[i,j]));
816	}
817	}
818	if( (double)(scalea)==(double)(0) )
819	{
820	scalea = 1;
821	}
822	scalea = 1/scalea;
823	cmatrixlusolveinternal(ref lua, ref p, scalea, n, ref a, true, ref b, m, ref info, ref rep, ref x);
824	}
825
826
827	/*************************************************************************
828	Dense solver. Same as RMatrixMixedSolve(), but for complex matrices.
829
830	Algorithm features:
831	* automatic detection of degenerate cases
832	* condition number estimation
833	* iterative refinement
834	* O(N^2) complexity
835
836	INPUT PARAMETERS
837	A - array[0..N-1,0..N-1], system matrix
838	LUA - array[0..N-1,0..N-1], LU decomposition, CMatrixLU result
839	P - array[0..N-1], pivots array, CMatrixLU result
840	N - size of A
841	B - array[0..N-1], right part
842
843	OUTPUT PARAMETERS
844	Info - same as in RMatrixSolveM
845	Rep - same as in RMatrixSolveM
846	X - same as in RMatrixSolveM
847
848	-- ALGLIB --
849	Copyright 27.01.2010 by Bochkanov Sergey
850	*************************************************************************/
851	public static void cmatrixmixedsolve(ref AP.Complex[,] a,
852	ref AP.Complex[,] lua,
853	ref int[] p,
854	int n,
855	ref AP.Complex[] b,
856	ref int info,
857	ref densesolverreport rep,
858	ref AP.Complex[] x)
859	{
860	AP.Complex[,] bm = new AP.Complex[0,0];
861	AP.Complex[,] xm = new AP.Complex[0,0];
862	int i_ = 0;
863
864	if( n<=0 )
865	{
866	info = -1;
867	return;
868	}
869	bm = new AP.Complex[n, 1];
870	for(i_=0; i_<=n-1;i_++)
871	{
872	bm[i_,0] = b[i_];
873	}
874	cmatrixmixedsolvem(ref a, ref lua, ref p, n, ref bm, 1, ref info, ref rep, ref xm);
875	x = new AP.Complex[n];
876	for(i_=0; i_<=n-1;i_++)
877	{
878	x[i_] = xm[i_,0];
879	}
880	}
881
882
883	/*************************************************************************
884	Dense solver. Same as RMatrixSolveM(), but for symmetric positive definite
885	matrices.
886
887	Algorithm features:
888	* automatic detection of degenerate cases
889	* condition number estimation
890	* O(N^3+M*N^2) complexity
891	* matrix is represented by its upper or lower triangle
892
893	No iterative refinement is provided because such partial representation of
894	matrix does not allow efficient calculation of extra-precise matrix-vector
895	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
896	need iterative refinement.
897
898	INPUT PARAMETERS
899	A - array[0..N-1,0..N-1], system matrix
900	N - size of A
901	IsUpper - what half of A is provided
902	B - array[0..N-1,0..M-1], right part
903	M - right part size
904
905	OUTPUT PARAMETERS
906	Info - same as in RMatrixSolve.
907	Returns -3 for non-SPD matrices.
908	Rep - same as in RMatrixSolve
909	X - same as in RMatrixSolve
910
911	-- ALGLIB --
912	Copyright 27.01.2010 by Bochkanov Sergey
913	*************************************************************************/
914	public static void spdmatrixsolvem(ref double[,] a,
915	int n,
916	bool isupper,
917	ref double[,] b,
918	int m,
919	ref int info,
920	ref densesolverreport rep,
921	ref double[,] x)
922	{
923	double[,] da = new double[0,0];
924	double sqrtscalea = 0;
925	int i = 0;
926	int j = 0;
927	int j1 = 0;
928	int j2 = 0;
929	int i_ = 0;
930
931
932	//
933	// prepare: check inputs, allocate space...
934	//
935	if( n<=0 \| m<=0 )
936	{
937	info = -1;
938	return;
939	}
940	da = new double[n, n];
941
942	//
943	// 1. scale matrix, max(\|A[i,j]\|)
944	// 2. factorize scaled matrix
945	// 3. solve
946	//
947	sqrtscalea = 0;
948	for(i=0; i<=n-1; i++)
949	{
950	if( isupper )
951	{
952	j1 = i;
953	j2 = n-1;
954	}
955	else
956	{
957	j1 = 0;
958	j2 = i;
959	}
960	for(j=j1; j<=j2; j++)
961	{
962	sqrtscalea = Math.Max(sqrtscalea, Math.Abs(a[i,j]));
963	}
964	}
965	if( (double)(sqrtscalea)==(double)(0) )
966	{
967	sqrtscalea = 1;
968	}
969	sqrtscalea = 1/sqrtscalea;
970	sqrtscalea = Math.Sqrt(sqrtscalea);
971	for(i=0; i<=n-1; i++)
972	{
973	if( isupper )
974	{
975	j1 = i;
976	j2 = n-1;
977	}
978	else
979	{
980	j1 = 0;
981	j2 = i;
982	}
983	for(i_=j1; i_<=j2;i_++)
984	{
985	da[i,i_] = a[i,i_];
986	}
987	}
988	if( !trfac.spdmatrixcholesky(ref da, n, isupper) )
989	{
990	x = new double[n, m];
991	for(i=0; i<=n-1; i++)
992	{
993	for(j=0; j<=m-1; j++)
994	{
995	x[i,j] = 0;
996	}
997	}
998	rep.r1 = 0;
999	rep.rinf = 0;
1000	info = -3;
1001	return;
1002	}
1003	info = 1;
1004	spdmatrixcholeskysolveinternal(ref da, sqrtscalea, n, isupper, ref a, true, ref b, m, ref info, ref rep, ref x);
1005	}
1006
1007
1008	/*************************************************************************
1009	Dense solver. Same as RMatrixSolve(), but for SPD matrices.
1010
1011	Algorithm features:
1012	* automatic detection of degenerate cases
1013	* condition number estimation
1014	* O(N^3) complexity
1015	* matrix is represented by its upper or lower triangle
1016
1017	No iterative refinement is provided because such partial representation of
1018	matrix does not allow efficient calculation of extra-precise matrix-vector
1019	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1020	need iterative refinement.
1021
1022	INPUT PARAMETERS
1023	A - array[0..N-1,0..N-1], system matrix
1024	N - size of A
1025	IsUpper - what half of A is provided
1026	B - array[0..N-1], right part
1027
1028	OUTPUT PARAMETERS
1029	Info - same as in RMatrixSolve
1030	Returns -3 for non-SPD matrices.
1031	Rep - same as in RMatrixSolve
1032	X - same as in RMatrixSolve
1033
1034	-- ALGLIB --
1035	Copyright 27.01.2010 by Bochkanov Sergey
1036	*************************************************************************/
1037	public static void spdmatrixsolve(ref double[,] a,
1038	int n,
1039	bool isupper,
1040	ref double[] b,
1041	ref int info,
1042	ref densesolverreport rep,
1043	ref double[] x)
1044	{
1045	double[,] bm = new double[0,0];
1046	double[,] xm = new double[0,0];
1047	int i_ = 0;
1048
1049	if( n<=0 )
1050	{
1051	info = -1;
1052	return;
1053	}
1054	bm = new double[n, 1];
1055	for(i_=0; i_<=n-1;i_++)
1056	{
1057	bm[i_,0] = b[i_];
1058	}
1059	spdmatrixsolvem(ref a, n, isupper, ref bm, 1, ref info, ref rep, ref xm);
1060	x = new double[n];
1061	for(i_=0; i_<=n-1;i_++)
1062	{
1063	x[i_] = xm[i_,0];
1064	}
1065	}
1066
1067
1068	/*************************************************************************
1069	Dense solver. Same as RMatrixLUSolveM(), but for SPD matrices represented
1070	by their Cholesky decomposition.
1071
1072	Algorithm features:
1073	* automatic detection of degenerate cases
1074	* O(M*N^2) complexity
1075	* condition number estimation
1076	* matrix is represented by its upper or lower triangle
1077
1078	No iterative refinement is provided because such partial representation of
1079	matrix does not allow efficient calculation of extra-precise matrix-vector
1080	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1081	need iterative refinement.
1082
1083	INPUT PARAMETERS
1084	CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1085	SPDMatrixCholesky result
1086	N - size of CHA
1087	IsUpper - what half of CHA is provided
1088	B - array[0..N-1,0..M-1], right part
1089	M - right part size
1090
1091	OUTPUT PARAMETERS
1092	Info - same as in RMatrixSolve
1093	Rep - same as in RMatrixSolve
1094	X - same as in RMatrixSolve
1095
1096	-- ALGLIB --
1097	Copyright 27.01.2010 by Bochkanov Sergey
1098	*************************************************************************/
1099	public static void spdmatrixcholeskysolvem(ref double[,] cha,
1100	int n,
1101	bool isupper,
1102	ref double[,] b,
1103	int m,
1104	ref int info,
1105	ref densesolverreport rep,
1106	ref double[,] x)
1107	{
1108	double[,] emptya = new double[0,0];
1109	double sqrtscalea = 0;
1110	int i = 0;
1111	int j = 0;
1112	int j1 = 0;
1113	int j2 = 0;
1114
1115
1116	//
1117	// prepare: check inputs, allocate space...
1118	//
1119	if( n<=0 \| m<=0 )
1120	{
1121	info = -1;
1122	return;
1123	}
1124
1125	//
1126	// 1. scale matrix, max(\|U[i,j]\|)
1127	// 2. factorize scaled matrix
1128	// 3. solve
1129	//
1130	sqrtscalea = 0;
1131	for(i=0; i<=n-1; i++)
1132	{
1133	if( isupper )
1134	{
1135	j1 = i;
1136	j2 = n-1;
1137	}
1138	else
1139	{
1140	j1 = 0;
1141	j2 = i;
1142	}
1143	for(j=j1; j<=j2; j++)
1144	{
1145	sqrtscalea = Math.Max(sqrtscalea, Math.Abs(cha[i,j]));
1146	}
1147	}
1148	if( (double)(sqrtscalea)==(double)(0) )
1149	{
1150	sqrtscalea = 1;
1151	}
1152	sqrtscalea = 1/sqrtscalea;
1153	spdmatrixcholeskysolveinternal(ref cha, sqrtscalea, n, isupper, ref emptya, false, ref b, m, ref info, ref rep, ref x);
1154	}
1155
1156
1157	/*************************************************************************
1158	Dense solver. Same as RMatrixLUSolve(), but for SPD matrices represented
1159	by their Cholesky decomposition.
1160
1161	Algorithm features:
1162	* automatic detection of degenerate cases
1163	* O(N^2) complexity
1164	* condition number estimation
1165	* matrix is represented by its upper or lower triangle
1166
1167	No iterative refinement is provided because such partial representation of
1168	matrix does not allow efficient calculation of extra-precise matrix-vector
1169	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1170	need iterative refinement.
1171
1172	INPUT PARAMETERS
1173	CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1174	SPDMatrixCholesky result
1175	N - size of A
1176	IsUpper - what half of CHA is provided
1177	B - array[0..N-1], right part
1178
1179	OUTPUT PARAMETERS
1180	Info - same as in RMatrixSolve
1181	Rep - same as in RMatrixSolve
1182	X - same as in RMatrixSolve
1183
1184	-- ALGLIB --
1185	Copyright 27.01.2010 by Bochkanov Sergey
1186	*************************************************************************/
1187	public static void spdmatrixcholeskysolve(ref double[,] cha,
1188	int n,
1189	bool isupper,
1190	ref double[] b,
1191	ref int info,
1192	ref densesolverreport rep,
1193	ref double[] x)
1194	{
1195	double[,] bm = new double[0,0];
1196	double[,] xm = new double[0,0];
1197	int i_ = 0;
1198
1199	if( n<=0 )
1200	{
1201	info = -1;
1202	return;
1203	}
1204	bm = new double[n, 1];
1205	for(i_=0; i_<=n-1;i_++)
1206	{
1207	bm[i_,0] = b[i_];
1208	}
1209	spdmatrixcholeskysolvem(ref cha, n, isupper, ref bm, 1, ref info, ref rep, ref xm);
1210	x = new double[n];
1211	for(i_=0; i_<=n-1;i_++)
1212	{
1213	x[i_] = xm[i_,0];
1214	}
1215	}
1216
1217
1218	/*************************************************************************
1219	Dense solver. Same as RMatrixSolveM(), but for Hermitian positive definite
1220	matrices.
1221
1222	Algorithm features:
1223	* automatic detection of degenerate cases
1224	* condition number estimation
1225	* O(N^3+M*N^2) complexity
1226	* matrix is represented by its upper or lower triangle
1227
1228	No iterative refinement is provided because such partial representation of
1229	matrix does not allow efficient calculation of extra-precise matrix-vector
1230	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1231	need iterative refinement.
1232
1233	INPUT PARAMETERS
1234	A - array[0..N-1,0..N-1], system matrix
1235	N - size of A
1236	IsUpper - what half of A is provided
1237	B - array[0..N-1,0..M-1], right part
1238	M - right part size
1239
1240	OUTPUT PARAMETERS
1241	Info - same as in RMatrixSolve.
1242	Returns -3 for non-HPD matrices.
1243	Rep - same as in RMatrixSolve
1244	X - same as in RMatrixSolve
1245
1246	-- ALGLIB --
1247	Copyright 27.01.2010 by Bochkanov Sergey
1248	*************************************************************************/
1249	public static void hpdmatrixsolvem(ref AP.Complex[,] a,
1250	int n,
1251	bool isupper,
1252	ref AP.Complex[,] b,
1253	int m,
1254	ref int info,
1255	ref densesolverreport rep,
1256	ref AP.Complex[,] x)
1257	{
1258	AP.Complex[,] da = new AP.Complex[0,0];
1259	double sqrtscalea = 0;
1260	int i = 0;
1261	int j = 0;
1262	int j1 = 0;
1263	int j2 = 0;
1264	int i_ = 0;
1265
1266
1267	//
1268	// prepare: check inputs, allocate space...
1269	//
1270	if( n<=0 \| m<=0 )
1271	{
1272	info = -1;
1273	return;
1274	}
1275	da = new AP.Complex[n, n];
1276
1277	//
1278	// 1. scale matrix, max(\|A[i,j]\|)
1279	// 2. factorize scaled matrix
1280	// 3. solve
1281	//
1282	sqrtscalea = 0;
1283	for(i=0; i<=n-1; i++)
1284	{
1285	if( isupper )
1286	{
1287	j1 = i;
1288	j2 = n-1;
1289	}
1290	else
1291	{
1292	j1 = 0;
1293	j2 = i;
1294	}
1295	for(j=j1; j<=j2; j++)
1296	{
1297	sqrtscalea = Math.Max(sqrtscalea, AP.Math.AbsComplex(a[i,j]));
1298	}
1299	}
1300	if( (double)(sqrtscalea)==(double)(0) )
1301	{
1302	sqrtscalea = 1;
1303	}
1304	sqrtscalea = 1/sqrtscalea;
1305	sqrtscalea = Math.Sqrt(sqrtscalea);
1306	for(i=0; i<=n-1; i++)
1307	{
1308	if( isupper )
1309	{
1310	j1 = i;
1311	j2 = n-1;
1312	}
1313	else
1314	{
1315	j1 = 0;
1316	j2 = i;
1317	}
1318	for(i_=j1; i_<=j2;i_++)
1319	{
1320	da[i,i_] = a[i,i_];
1321	}
1322	}
1323	if( !trfac.hpdmatrixcholesky(ref da, n, isupper) )
1324	{
1325	x = new AP.Complex[n, m];
1326	for(i=0; i<=n-1; i++)
1327	{
1328	for(j=0; j<=m-1; j++)
1329	{
1330	x[i,j] = 0;
1331	}
1332	}
1333	rep.r1 = 0;
1334	rep.rinf = 0;
1335	info = -3;
1336	return;
1337	}
1338	info = 1;
1339	hpdmatrixcholeskysolveinternal(ref da, sqrtscalea, n, isupper, ref a, true, ref b, m, ref info, ref rep, ref x);
1340	}
1341
1342
1343	/*************************************************************************
1344	Dense solver. Same as RMatrixSolve(), but for Hermitian positive definite
1345	matrices.
1346
1347	Algorithm features:
1348	* automatic detection of degenerate cases
1349	* condition number estimation
1350	* O(N^3) complexity
1351	* matrix is represented by its upper or lower triangle
1352
1353	No iterative refinement is provided because such partial representation of
1354	matrix does not allow efficient calculation of extra-precise matrix-vector
1355	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1356	need iterative refinement.
1357
1358	INPUT PARAMETERS
1359	A - array[0..N-1,0..N-1], system matrix
1360	N - size of A
1361	IsUpper - what half of A is provided
1362	B - array[0..N-1], right part
1363
1364	OUTPUT PARAMETERS
1365	Info - same as in RMatrixSolve
1366	Returns -3 for non-HPD matrices.
1367	Rep - same as in RMatrixSolve
1368	X - same as in RMatrixSolve
1369
1370	-- ALGLIB --
1371	Copyright 27.01.2010 by Bochkanov Sergey
1372	*************************************************************************/
1373	public static void hpdmatrixsolve(ref AP.Complex[,] a,
1374	int n,
1375	bool isupper,
1376	ref AP.Complex[] b,
1377	ref int info,
1378	ref densesolverreport rep,
1379	ref AP.Complex[] x)
1380	{
1381	AP.Complex[,] bm = new AP.Complex[0,0];
1382	AP.Complex[,] xm = new AP.Complex[0,0];
1383	int i_ = 0;
1384
1385	if( n<=0 )
1386	{
1387	info = -1;
1388	return;
1389	}
1390	bm = new AP.Complex[n, 1];
1391	for(i_=0; i_<=n-1;i_++)
1392	{
1393	bm[i_,0] = b[i_];
1394	}
1395	hpdmatrixsolvem(ref a, n, isupper, ref bm, 1, ref info, ref rep, ref xm);
1396	x = new AP.Complex[n];
1397	for(i_=0; i_<=n-1;i_++)
1398	{
1399	x[i_] = xm[i_,0];
1400	}
1401	}
1402
1403
1404	/*************************************************************************
1405	Dense solver. Same as RMatrixLUSolveM(), but for HPD matrices represented
1406	by their Cholesky decomposition.
1407
1408	Algorithm features:
1409	* automatic detection of degenerate cases
1410	* O(M*N^2) complexity
1411	* condition number estimation
1412	* matrix is represented by its upper or lower triangle
1413
1414	No iterative refinement is provided because such partial representation of
1415	matrix does not allow efficient calculation of extra-precise matrix-vector
1416	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1417	need iterative refinement.
1418
1419	INPUT PARAMETERS
1420	CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1421	HPDMatrixCholesky result
1422	N - size of CHA
1423	IsUpper - what half of CHA is provided
1424	B - array[0..N-1,0..M-1], right part
1425	M - right part size
1426
1427	OUTPUT PARAMETERS
1428	Info - same as in RMatrixSolve
1429	Rep - same as in RMatrixSolve
1430	X - same as in RMatrixSolve
1431
1432	-- ALGLIB --
1433	Copyright 27.01.2010 by Bochkanov Sergey
1434	*************************************************************************/
1435	public static void hpdmatrixcholeskysolvem(ref AP.Complex[,] cha,
1436	int n,
1437	bool isupper,
1438	ref AP.Complex[,] b,
1439	int m,
1440	ref int info,
1441	ref densesolverreport rep,
1442	ref AP.Complex[,] x)
1443	{
1444	AP.Complex[,] emptya = new AP.Complex[0,0];
1445	double sqrtscalea = 0;
1446	int i = 0;
1447	int j = 0;
1448	int j1 = 0;
1449	int j2 = 0;
1450
1451
1452	//
1453	// prepare: check inputs, allocate space...
1454	//
1455	if( n<=0 \| m<=0 )
1456	{
1457	info = -1;
1458	return;
1459	}
1460
1461	//
1462	// 1. scale matrix, max(\|U[i,j]\|)
1463	// 2. factorize scaled matrix
1464	// 3. solve
1465	//
1466	sqrtscalea = 0;
1467	for(i=0; i<=n-1; i++)
1468	{
1469	if( isupper )
1470	{
1471	j1 = i;
1472	j2 = n-1;
1473	}
1474	else
1475	{
1476	j1 = 0;
1477	j2 = i;
1478	}
1479	for(j=j1; j<=j2; j++)
1480	{
1481	sqrtscalea = Math.Max(sqrtscalea, AP.Math.AbsComplex(cha[i,j]));
1482	}
1483	}
1484	if( (double)(sqrtscalea)==(double)(0) )
1485	{
1486	sqrtscalea = 1;
1487	}
1488	sqrtscalea = 1/sqrtscalea;
1489	hpdmatrixcholeskysolveinternal(ref cha, sqrtscalea, n, isupper, ref emptya, false, ref b, m, ref info, ref rep, ref x);
1490	}
1491
1492
1493	/*************************************************************************
1494	Dense solver. Same as RMatrixLUSolve(), but for HPD matrices represented
1495	by their Cholesky decomposition.
1496
1497	Algorithm features:
1498	* automatic detection of degenerate cases
1499	* O(N^2) complexity
1500	* condition number estimation
1501	* matrix is represented by its upper or lower triangle
1502
1503	No iterative refinement is provided because such partial representation of
1504	matrix does not allow efficient calculation of extra-precise matrix-vector
1505	products for large matrices. Use RMatrixSolve or RMatrixMixedSolve if you
1506	need iterative refinement.
1507
1508	INPUT PARAMETERS
1509	CHA - array[0..N-1,0..N-1], Cholesky decomposition,
1510	SPDMatrixCholesky result
1511	N - size of A
1512	IsUpper - what half of CHA is provided
1513	B - array[0..N-1], right part
1514
1515	OUTPUT PARAMETERS
1516	Info - same as in RMatrixSolve
1517	Rep - same as in RMatrixSolve
1518	X - same as in RMatrixSolve
1519
1520	-- ALGLIB --
1521	Copyright 27.01.2010 by Bochkanov Sergey
1522	*************************************************************************/
1523	public static void hpdmatrixcholeskysolve(ref AP.Complex[,] cha,
1524	int n,
1525	bool isupper,
1526	ref AP.Complex[] b,
1527	ref int info,
1528	ref densesolverreport rep,
1529	ref AP.Complex[] x)
1530	{
1531	AP.Complex[,] bm = new AP.Complex[0,0];
1532	AP.Complex[,] xm = new AP.Complex[0,0];
1533	int i_ = 0;
1534
1535	if( n<=0 )
1536	{
1537	info = -1;
1538	return;
1539	}
1540	bm = new AP.Complex[n, 1];
1541	for(i_=0; i_<=n-1;i_++)
1542	{
1543	bm[i_,0] = b[i_];
1544	}
1545	hpdmatrixcholeskysolvem(ref cha, n, isupper, ref bm, 1, ref info, ref rep, ref xm);
1546	x = new AP.Complex[n];
1547	for(i_=0; i_<=n-1;i_++)
1548	{
1549	x[i_] = xm[i_,0];
1550	}
1551	}
1552
1553
1554	/*************************************************************************
1555	Dense solver.
1556
1557	This subroutine finds solution of the linear system A*X=B with non-square,
1558	possibly degenerate A. System is solved in the least squares sense, and
1559	general least squares solution X = X0 + CXy which minimizes \|AX-B\| is
1560	returned. If A is non-degenerate, solution in the usual sense is returned
1561
1562	Algorithm features:
1563	* automatic detection of degenerate cases
1564	* iterative refinement
1565	* O(N^3) complexity
1566
1567	INPUT PARAMETERS
1568	A - array[0..NRows-1,0..NCols-1], system matrix
1569	NRows - vertical size of A
1570	NCols - horizontal size of A
1571	B - array[0..NCols-1], right part
1572	Threshold- a number in [0,1]. Singular values beyond Threshold are
1573	considered zero. Set it to 0.0, if you don't understand
1574	what it means, so the solver will choose good value on its
1575	own.
1576
1577	OUTPUT PARAMETERS
1578	Info - return code:
1579	* -4 SVD subroutine failed
1580	* -1 if NRows<=0 or NCols<=0 or Threshold<0 was passed
1581	* 1 if task is solved
1582	Rep - solver report, see below for more info
1583	X - array[0..N-1,0..M-1], it contains:
1584	* solution of A*X=B if A is non-singular (well-conditioned
1585	or ill-conditioned, but not very close to singular)
1586	* zeros, if A is singular or VERY close to singular
1587	(in this case Info=-3).
1588
1589	SOLVER REPORT
1590
1591	Subroutine sets following fields of the Rep structure:
1592	* R2 reciprocal of condition number: 1/cond(A), 2-norm.
1593	* N = NCols
1594	* K dim(Null(A))
1595	* CX array[0..N-1,0..K-1], kernel of A.
1596	Columns of CX store such vectors that A*CX[i]=0.
1597
1598	-- ALGLIB --
1599	Copyright 24.08.2009 by Bochkanov Sergey
1600	*************************************************************************/
1601	public static void rmatrixsolvels(ref double[,] a,
1602	int nrows,
1603	int ncols,
1604	ref double[] b,
1605	double threshold,
1606	ref int info,
1607	ref densesolverlsreport rep,
1608	ref double[] x)
1609	{
1610	double[] sv = new double[0];
1611	double[,] u = new double[0,0];
1612	double[,] vt = new double[0,0];
1613	double[] rp = new double[0];
1614	double[] utb = new double[0];
1615	double[] sutb = new double[0];
1616	double[] tmp = new double[0];
1617	double[] ta = new double[0];
1618	double[] tx = new double[0];
1619	double[] buf = new double[0];
1620	double[] w = new double[0];
1621	int i = 0;
1622	int j = 0;
1623	int nsv = 0;
1624	int kernelidx = 0;
1625	double v = 0;
1626	double verr = 0;
1627	bool svdfailed = new bool();
1628	bool zeroa = new bool();
1629	int rfs = 0;
1630	int nrfs = 0;
1631	bool terminatenexttime = new bool();
1632	bool smallerr = new bool();
1633	int i_ = 0;
1634
1635	if( nrows<=0 \| ncols<=0 \| (double)(threshold)<(double)(0) )
1636	{
1637	info = -1;
1638	return;
1639	}
1640	if( (double)(threshold)==(double)(0) )
1641	{
1642	threshold = 1000*AP.Math.MachineEpsilon;
1643	}
1644
1645	//
1646	// Factorize A first
1647	//
1648	svdfailed = !svd.rmatrixsvd(a, nrows, ncols, 1, 2, 2, ref sv, ref u, ref vt);
1649	zeroa = (double)(sv[0])==(double)(0);
1650	if( svdfailed \| zeroa )
1651	{
1652	if( svdfailed )
1653	{
1654	info = -4;
1655	}
1656	else
1657	{
1658	info = 1;
1659	}
1660	x = new double[ncols];
1661	for(i=0; i<=ncols-1; i++)
1662	{
1663	x[i] = 0;
1664	}
1665	rep.n = ncols;
1666	rep.k = ncols;
1667	rep.cx = new double[ncols, ncols];
1668	for(i=0; i<=ncols-1; i++)
1669	{
1670	for(j=0; j<=ncols-1; j++)
1671	{
1672	if( i==j )
1673	{
1674	rep.cx[i,j] = 1;
1675	}
1676	else
1677	{
1678	rep.cx[i,j] = 0;
1679	}
1680	}
1681	}
1682	rep.r2 = 0;
1683	return;
1684	}
1685	nsv = Math.Min(ncols, nrows);
1686	if( nsv==ncols )
1687	{
1688	rep.r2 = sv[nsv-1]/sv[0];
1689	}
1690	else
1691	{
1692	rep.r2 = 0;
1693	}
1694	rep.n = ncols;
1695	info = 1;
1696
1697	//
1698	// Iterative refinement of xc combined with solution:
1699	// 1. xc = 0
1700	// 2. calculate r = bc-A*xc using extra-precise dot product
1701	// 3. solve A*y = r
1702	// 4. update x:=x+r
1703	// 5. goto 2
1704	//
1705	// This cycle is executed until one of two things happens:
1706	// 1. maximum number of iterations reached
1707	// 2. last iteration decreased error to the lower limit
1708	//
1709	utb = new double[nsv];
1710	sutb = new double[nsv];
1711	x = new double[ncols];
1712	tmp = new double[ncols];
1713	ta = new double[ncols+1];
1714	tx = new double[ncols+1];
1715	buf = new double[ncols+1];
1716	for(i=0; i<=ncols-1; i++)
1717	{
1718	x[i] = 0;
1719	}
1720	kernelidx = nsv;
1721	for(i=0; i<=nsv-1; i++)
1722	{
1723	if( (double)(sv[i])<=(double)(threshold*sv[0]) )
1724	{
1725	kernelidx = i;
1726	break;
1727	}
1728	}
1729	rep.k = ncols-kernelidx;
1730	nrfs = densesolverrfsmaxv2(ncols, rep.r2);
1731	terminatenexttime = false;
1732	rp = new double[nrows];
1733	for(rfs=0; rfs<=nrfs; rfs++)
1734	{
1735	if( terminatenexttime )
1736	{
1737	break;
1738	}
1739
1740	//
1741	// calculate right part
1742	//
1743	if( rfs==0 )
1744	{
1745	for(i_=0; i_<=nrows-1;i_++)
1746	{
1747	rp[i_] = b[i_];
1748	}
1749	}
1750	else
1751	{
1752	smallerr = true;
1753	for(i=0; i<=nrows-1; i++)
1754	{
1755	for(i_=0; i_<=ncols-1;i_++)
1756	{
1757	ta[i_] = a[i,i_];
1758	}
1759	ta[ncols] = -1;
1760	for(i_=0; i_<=ncols-1;i_++)
1761	{
1762	tx[i_] = x[i_];
1763	}
1764	tx[ncols] = b[i];
1765	xblas.xdot(ref ta, ref tx, ncols+1, ref buf, ref v, ref verr);
1766	rp[i] = -v;
1767	smallerr = smallerr & (double)(Math.Abs(v))<(double)(4*verr);
1768	}
1769	if( smallerr )
1770	{
1771	terminatenexttime = true;
1772	}
1773	}
1774
1775	//
1776	// solve A*dx = rp
1777	//
1778	for(i=0; i<=ncols-1; i++)
1779	{
1780	tmp[i] = 0;
1781	}
1782	for(i=0; i<=nsv-1; i++)
1783	{
1784	utb[i] = 0;
1785	}
1786	for(i=0; i<=nrows-1; i++)
1787	{
1788	v = rp[i];
1789	for(i_=0; i_<=nsv-1;i_++)
1790	{
1791	utb[i_] = utb[i_] + v*u[i,i_];
1792	}
1793	}
1794	for(i=0; i<=nsv-1; i++)
1795	{
1796	if( i<kernelidx )
1797	{
1798	sutb[i] = utb[i]/sv[i];
1799	}
1800	else
1801	{
1802	sutb[i] = 0;
1803	}
1804	}
1805	for(i=0; i<=nsv-1; i++)
1806	{
1807	v = sutb[i];
1808	for(i_=0; i_<=ncols-1;i_++)
1809	{
1810	tmp[i_] = tmp[i_] + v*vt[i,i_];
1811	}
1812	}
1813
1814	//
1815	// update x: x:=x+dx
1816	//
1817	for(i_=0; i_<=ncols-1;i_++)
1818	{
1819	x[i_] = x[i_] + tmp[i_];
1820	}
1821	}
1822
1823	//
1824	// fill CX
1825	//
1826	if( rep.k>0 )
1827	{
1828	rep.cx = new double[ncols, rep.k];
1829	for(i=0; i<=rep.k-1; i++)
1830	{
1831	for(i_=0; i_<=ncols-1;i_++)
1832	{
1833	rep.cx[i_,i] = vt[kernelidx+i,i_];
1834	}
1835	}
1836	}
1837	}
1838
1839
1840	/*************************************************************************
1841	Internal LU solver
1842
1843	-- ALGLIB --
1844	Copyright 27.01.2010 by Bochkanov Sergey
1845	*************************************************************************/
1846	private static void rmatrixlusolveinternal(ref double[,] lua,
1847	ref int[] p,
1848	double scalea,
1849	int n,
1850	ref double[,] a,
1851	bool havea,
1852	ref double[,] b,
1853	int m,
1854	ref int info,
1855	ref densesolverreport rep,
1856	ref double[,] x)
1857	{
1858	int i = 0;
1859	int j = 0;
1860	int k = 0;
1861	int rfs = 0;
1862	int nrfs = 0;
1863	double[] xc = new double[0];
1864	double[] y = new double[0];
1865	double[] bc = new double[0];
1866	double[] xa = new double[0];
1867	double[] xb = new double[0];
1868	double[] tx = new double[0];
1869	double v = 0;
1870	double verr = 0;
1871	double mxb = 0;
1872	double scaleright = 0;
1873	bool smallerr = new bool();
1874	bool terminatenexttime = new bool();
1875	int i_ = 0;
1876
1877	System.Diagnostics.Debug.Assert((double)(scalea)>(double)(0));
1878
1879	//
1880	// prepare: check inputs, allocate space...
1881	//
1882	if( n<=0 \| m<=0 )
1883	{
1884	info = -1;
1885	return;
1886	}
1887	for(i=0; i<=n-1; i++)
1888	{
1889	if( p[i]>n-1 \| p[i]<i )
1890	{
1891	info = -1;
1892	return;
1893	}
1894	}
1895	x = new double[n, m];
1896	y = new double[n];
1897	xc = new double[n];
1898	bc = new double[n];
1899	tx = new double[n+1];
1900	xa = new double[n+1];
1901	xb = new double[n+1];
1902
1903	//
1904	// estimate condition number, test for near singularity
1905	//
1906	rep.r1 = rcond.rmatrixlurcond1(ref lua, n);
1907	rep.rinf = rcond.rmatrixlurcondinf(ref lua, n);
1908	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
1909	{
1910	for(i=0; i<=n-1; i++)
1911	{
1912	for(j=0; j<=m-1; j++)
1913	{
1914	x[i,j] = 0;
1915	}
1916	}
1917	rep.r1 = 0;
1918	rep.rinf = 0;
1919	info = -3;
1920	return;
1921	}
1922	info = 1;
1923
1924	//
1925	// solve
1926	//
1927	for(k=0; k<=m-1; k++)
1928	{
1929
1930	//
1931	// copy B to contiguous storage
1932	//
1933	for(i_=0; i_<=n-1;i_++)
1934	{
1935	bc[i_] = b[i_,k];
1936	}
1937
1938	//
1939	// Scale right part:
1940	// * MX stores max(\|Bi\|)
1941	// * ScaleRight stores actual scaling applied to B when solving systems
1942	// it is chosen to make \|scaleRight*b\| close to 1.
1943	//
1944	mxb = 0;
1945	for(i=0; i<=n-1; i++)
1946	{
1947	mxb = Math.Max(mxb, Math.Abs(bc[i]));
1948	}
1949	if( (double)(mxb)==(double)(0) )
1950	{
1951	mxb = 1;
1952	}
1953	scaleright = 1/mxb;
1954
1955	//
1956	// First, non-iterative part of solution process.
1957	// We use separate code for this task because
1958	// XDot is quite slow and we want to save time.
1959	//
1960	for(i_=0; i_<=n-1;i_++)
1961	{
1962	xc[i_] = scaleright*bc[i_];
1963	}
1964	rbasiclusolve(ref lua, ref p, scalea, n, ref xc, ref tx);
1965
1966	//
1967	// Iterative refinement of xc:
1968	// * calculate r = bc-A*xc using extra-precise dot product
1969	// * solve A*y = r
1970	// * update x:=x+r
1971	//
1972	// This cycle is executed until one of two things happens:
1973	// 1. maximum number of iterations reached
1974	// 2. last iteration decreased error to the lower limit
1975	//
1976	if( havea )
1977	{
1978	nrfs = densesolverrfsmax(n, rep.r1, rep.rinf);
1979	terminatenexttime = false;
1980	for(rfs=0; rfs<=nrfs-1; rfs++)
1981	{
1982	if( terminatenexttime )
1983	{
1984	break;
1985	}
1986
1987	//
1988	// generate right part
1989	//
1990	smallerr = true;
1991	for(i_=0; i_<=n-1;i_++)
1992	{
1993	xb[i_] = xc[i_];
1994	}
1995	for(i=0; i<=n-1; i++)
1996	{
1997	for(i_=0; i_<=n-1;i_++)
1998	{
1999	xa[i_] = scalea*a[i,i_];
2000	}
2001	xa[n] = -1;
2002	xb[n] = scaleright*bc[i];
2003	xblas.xdot(ref xa, ref xb, n+1, ref tx, ref v, ref verr);
2004	y[i] = -v;
2005	smallerr = smallerr & (double)(Math.Abs(v))<(double)(4*verr);
2006	}
2007	if( smallerr )
2008	{
2009	terminatenexttime = true;
2010	}
2011
2012	//
2013	// solve and update
2014	//
2015	rbasiclusolve(ref lua, ref p, scalea, n, ref y, ref tx);
2016	for(i_=0; i_<=n-1;i_++)
2017	{
2018	xc[i_] = xc[i_] + y[i_];
2019	}
2020	}
2021	}
2022
2023	//
2024	// Store xc.
2025	// Post-scale result.
2026	//
2027	v = scalea*mxb;
2028	for(i_=0; i_<=n-1;i_++)
2029	{
2030	x[i_,k] = v*xc[i_];
2031	}
2032	}
2033	}
2034
2035
2036	/*************************************************************************
2037	Internal Cholesky solver
2038
2039	-- ALGLIB --
2040	Copyright 27.01.2010 by Bochkanov Sergey
2041	*************************************************************************/
2042	private static void spdmatrixcholeskysolveinternal(ref double[,] cha,
2043	double sqrtscalea,
2044	int n,
2045	bool isupper,
2046	ref double[,] a,
2047	bool havea,
2048	ref double[,] b,
2049	int m,
2050	ref int info,
2051	ref densesolverreport rep,
2052	ref double[,] x)
2053	{
2054	int i = 0;
2055	int j = 0;
2056	int k = 0;
2057	int rfs = 0;
2058	int nrfs = 0;
2059	double[] xc = new double[0];
2060	double[] y = new double[0];
2061	double[] bc = new double[0];
2062	double[] xa = new double[0];
2063	double[] xb = new double[0];
2064	double[] tx = new double[0];
2065	double v = 0;
2066	double verr = 0;
2067	double mxb = 0;
2068	double scaleright = 0;
2069	bool smallerr = new bool();
2070	bool terminatenexttime = new bool();
2071	int i_ = 0;
2072
2073	System.Diagnostics.Debug.Assert((double)(sqrtscalea)>(double)(0));
2074
2075	//
2076	// prepare: check inputs, allocate space...
2077	//
2078	if( n<=0 \| m<=0 )
2079	{
2080	info = -1;
2081	return;
2082	}
2083	x = new double[n, m];
2084	y = new double[n];
2085	xc = new double[n];
2086	bc = new double[n];
2087	tx = new double[n+1];
2088	xa = new double[n+1];
2089	xb = new double[n+1];
2090
2091	//
2092	// estimate condition number, test for near singularity
2093	//
2094	rep.r1 = rcond.spdmatrixcholeskyrcond(ref cha, n, isupper);
2095	rep.rinf = rep.r1;
2096	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) )
2097	{
2098	for(i=0; i<=n-1; i++)
2099	{
2100	for(j=0; j<=m-1; j++)
2101	{
2102	x[i,j] = 0;
2103	}
2104	}
2105	rep.r1 = 0;
2106	rep.rinf = 0;
2107	info = -3;
2108	return;
2109	}
2110	info = 1;
2111
2112	//
2113	// solve
2114	//
2115	for(k=0; k<=m-1; k++)
2116	{
2117
2118	//
2119	// copy B to contiguous storage
2120	//
2121	for(i_=0; i_<=n-1;i_++)
2122	{
2123	bc[i_] = b[i_,k];
2124	}
2125
2126	//
2127	// Scale right part:
2128	// * MX stores max(\|Bi\|)
2129	// * ScaleRight stores actual scaling applied to B when solving systems
2130	// it is chosen to make \|scaleRight*b\| close to 1.
2131	//
2132	mxb = 0;
2133	for(i=0; i<=n-1; i++)
2134	{
2135	mxb = Math.Max(mxb, Math.Abs(bc[i]));
2136	}
2137	if( (double)(mxb)==(double)(0) )
2138	{
2139	mxb = 1;
2140	}
2141	scaleright = 1/mxb;
2142
2143	//
2144	// First, non-iterative part of solution process.
2145	// We use separate code for this task because
2146	// XDot is quite slow and we want to save time.
2147	//
2148	for(i_=0; i_<=n-1;i_++)
2149	{
2150	xc[i_] = scaleright*bc[i_];
2151	}
2152	spdbasiccholeskysolve(ref cha, sqrtscalea, n, isupper, ref xc, ref tx);
2153
2154	//
2155	// Store xc.
2156	// Post-scale result.
2157	//
2158	v = AP.Math.Sqr(sqrtscalea)*mxb;
2159	for(i_=0; i_<=n-1;i_++)
2160	{
2161	x[i_,k] = v*xc[i_];
2162	}
2163	}
2164	}
2165
2166
2167	/*************************************************************************
2168	Internal LU solver
2169
2170	-- ALGLIB --
2171	Copyright 27.01.2010 by Bochkanov Sergey
2172	*************************************************************************/
2173	private static void cmatrixlusolveinternal(ref AP.Complex[,] lua,
2174	ref int[] p,
2175	double scalea,
2176	int n,
2177	ref AP.Complex[,] a,
2178	bool havea,
2179	ref AP.Complex[,] b,
2180	int m,
2181	ref int info,
2182	ref densesolverreport rep,
2183	ref AP.Complex[,] x)
2184	{
2185	int i = 0;
2186	int j = 0;
2187	int k = 0;
2188	int rfs = 0;
2189	int nrfs = 0;
2190	AP.Complex[] xc = new AP.Complex[0];
2191	AP.Complex[] y = new AP.Complex[0];
2192	AP.Complex[] bc = new AP.Complex[0];
2193	AP.Complex[] xa = new AP.Complex[0];
2194	AP.Complex[] xb = new AP.Complex[0];
2195	AP.Complex[] tx = new AP.Complex[0];
2196	double[] tmpbuf = new double[0];
2197	AP.Complex v = 0;
2198	double verr = 0;
2199	double mxb = 0;
2200	double scaleright = 0;
2201	bool smallerr = new bool();
2202	bool terminatenexttime = new bool();
2203	int i_ = 0;
2204
2205	System.Diagnostics.Debug.Assert((double)(scalea)>(double)(0));
2206
2207	//
2208	// prepare: check inputs, allocate space...
2209	//
2210	if( n<=0 \| m<=0 )
2211	{
2212	info = -1;
2213	return;
2214	}
2215	for(i=0; i<=n-1; i++)
2216	{
2217	if( p[i]>n-1 \| p[i]<i )
2218	{
2219	info = -1;
2220	return;
2221	}
2222	}
2223	x = new AP.Complex[n, m];
2224	y = new AP.Complex[n];
2225	xc = new AP.Complex[n];
2226	bc = new AP.Complex[n];
2227	tx = new AP.Complex[n];
2228	xa = new AP.Complex[n+1];
2229	xb = new AP.Complex[n+1];
2230	tmpbuf = new double[2*n+2];
2231
2232	//
2233	// estimate condition number, test for near singularity
2234	//
2235	rep.r1 = rcond.cmatrixlurcond1(ref lua, n);
2236	rep.rinf = rcond.cmatrixlurcondinf(ref lua, n);
2237	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
2238	{
2239	for(i=0; i<=n-1; i++)
2240	{
2241	for(j=0; j<=m-1; j++)
2242	{
2243	x[i,j] = 0;
2244	}
2245	}
2246	rep.r1 = 0;
2247	rep.rinf = 0;
2248	info = -3;
2249	return;
2250	}
2251	info = 1;
2252
2253	//
2254	// solve
2255	//
2256	for(k=0; k<=m-1; k++)
2257	{
2258
2259	//
2260	// copy B to contiguous storage
2261	//
2262	for(i_=0; i_<=n-1;i_++)
2263	{
2264	bc[i_] = b[i_,k];
2265	}
2266
2267	//
2268	// Scale right part:
2269	// * MX stores max(\|Bi\|)
2270	// * ScaleRight stores actual scaling applied to B when solving systems
2271	// it is chosen to make \|scaleRight*b\| close to 1.
2272	//
2273	mxb = 0;
2274	for(i=0; i<=n-1; i++)
2275	{
2276	mxb = Math.Max(mxb, AP.Math.AbsComplex(bc[i]));
2277	}
2278	if( (double)(mxb)==(double)(0) )
2279	{
2280	mxb = 1;
2281	}
2282	scaleright = 1/mxb;
2283
2284	//
2285	// First, non-iterative part of solution process.
2286	// We use separate code for this task because
2287	// XDot is quite slow and we want to save time.
2288	//
2289	for(i_=0; i_<=n-1;i_++)
2290	{
2291	xc[i_] = scaleright*bc[i_];
2292	}
2293	cbasiclusolve(ref lua, ref p, scalea, n, ref xc, ref tx);
2294
2295	//
2296	// Iterative refinement of xc:
2297	// * calculate r = bc-A*xc using extra-precise dot product
2298	// * solve A*y = r
2299	// * update x:=x+r
2300	//
2301	// This cycle is executed until one of two things happens:
2302	// 1. maximum number of iterations reached
2303	// 2. last iteration decreased error to the lower limit
2304	//
2305	if( havea )
2306	{
2307	nrfs = densesolverrfsmax(n, rep.r1, rep.rinf);
2308	terminatenexttime = false;
2309	for(rfs=0; rfs<=nrfs-1; rfs++)
2310	{
2311	if( terminatenexttime )
2312	{
2313	break;
2314	}
2315
2316	//
2317	// generate right part
2318	//
2319	smallerr = true;
2320	for(i_=0; i_<=n-1;i_++)
2321	{
2322	xb[i_] = xc[i_];
2323	}
2324	for(i=0; i<=n-1; i++)
2325	{
2326	for(i_=0; i_<=n-1;i_++)
2327	{
2328	xa[i_] = scalea*a[i,i_];
2329	}
2330	xa[n] = -1;
2331	xb[n] = scaleright*bc[i];
2332	xblas.xcdot(ref xa, ref xb, n+1, ref tmpbuf, ref v, ref verr);
2333	y[i] = -v;
2334	smallerr = smallerr & (double)(AP.Math.AbsComplex(v))<(double)(4*verr);
2335	}
2336	if( smallerr )
2337	{
2338	terminatenexttime = true;
2339	}
2340
2341	//
2342	// solve and update
2343	//
2344	cbasiclusolve(ref lua, ref p, scalea, n, ref y, ref tx);
2345	for(i_=0; i_<=n-1;i_++)
2346	{
2347	xc[i_] = xc[i_] + y[i_];
2348	}
2349	}
2350	}
2351
2352	//
2353	// Store xc.
2354	// Post-scale result.
2355	//
2356	v = scalea*mxb;
2357	for(i_=0; i_<=n-1;i_++)
2358	{
2359	x[i_,k] = v*xc[i_];
2360	}
2361	}
2362	}
2363
2364
2365	/*************************************************************************
2366	Internal Cholesky solver
2367
2368	-- ALGLIB --
2369	Copyright 27.01.2010 by Bochkanov Sergey
2370	*************************************************************************/
2371	private static void hpdmatrixcholeskysolveinternal(ref AP.Complex[,] cha,
2372	double sqrtscalea,
2373	int n,
2374	bool isupper,
2375	ref AP.Complex[,] a,
2376	bool havea,
2377	ref AP.Complex[,] b,
2378	int m,
2379	ref int info,
2380	ref densesolverreport rep,
2381	ref AP.Complex[,] x)
2382	{
2383	int i = 0;
2384	int j = 0;
2385	int k = 0;
2386	int rfs = 0;
2387	int nrfs = 0;
2388	AP.Complex[] xc = new AP.Complex[0];
2389	AP.Complex[] y = new AP.Complex[0];
2390	AP.Complex[] bc = new AP.Complex[0];
2391	AP.Complex[] xa = new AP.Complex[0];
2392	AP.Complex[] xb = new AP.Complex[0];
2393	AP.Complex[] tx = new AP.Complex[0];
2394	double v = 0;
2395	double verr = 0;
2396	double mxb = 0;
2397	double scaleright = 0;
2398	bool smallerr = new bool();
2399	bool terminatenexttime = new bool();
2400	int i_ = 0;
2401
2402	System.Diagnostics.Debug.Assert((double)(sqrtscalea)>(double)(0));
2403
2404	//
2405	// prepare: check inputs, allocate space...
2406	//
2407	if( n<=0 \| m<=0 )
2408	{
2409	info = -1;
2410	return;
2411	}
2412	x = new AP.Complex[n, m];
2413	y = new AP.Complex[n];
2414	xc = new AP.Complex[n];
2415	bc = new AP.Complex[n];
2416	tx = new AP.Complex[n+1];
2417	xa = new AP.Complex[n+1];
2418	xb = new AP.Complex[n+1];
2419
2420	//
2421	// estimate condition number, test for near singularity
2422	//
2423	rep.r1 = rcond.hpdmatrixcholeskyrcond(ref cha, n, isupper);
2424	rep.rinf = rep.r1;
2425	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) )
2426	{
2427	for(i=0; i<=n-1; i++)
2428	{
2429	for(j=0; j<=m-1; j++)
2430	{
2431	x[i,j] = 0;
2432	}
2433	}
2434	rep.r1 = 0;
2435	rep.rinf = 0;
2436	info = -3;
2437	return;
2438	}
2439	info = 1;
2440
2441	//
2442	// solve
2443	//
2444	for(k=0; k<=m-1; k++)
2445	{
2446
2447	//
2448	// copy B to contiguous storage
2449	//
2450	for(i_=0; i_<=n-1;i_++)
2451	{
2452	bc[i_] = b[i_,k];
2453	}
2454
2455	//
2456	// Scale right part:
2457	// * MX stores max(\|Bi\|)
2458	// * ScaleRight stores actual scaling applied to B when solving systems
2459	// it is chosen to make \|scaleRight*b\| close to 1.
2460	//
2461	mxb = 0;
2462	for(i=0; i<=n-1; i++)
2463	{
2464	mxb = Math.Max(mxb, AP.Math.AbsComplex(bc[i]));
2465	}
2466	if( (double)(mxb)==(double)(0) )
2467	{
2468	mxb = 1;
2469	}
2470	scaleright = 1/mxb;
2471
2472	//
2473	// First, non-iterative part of solution process.
2474	// We use separate code for this task because
2475	// XDot is quite slow and we want to save time.
2476	//
2477	for(i_=0; i_<=n-1;i_++)
2478	{
2479	xc[i_] = scaleright*bc[i_];
2480	}
2481	hpdbasiccholeskysolve(ref cha, sqrtscalea, n, isupper, ref xc, ref tx);
2482
2483	//
2484	// Store xc.
2485	// Post-scale result.
2486	//
2487	v = AP.Math.Sqr(sqrtscalea)*mxb;
2488	for(i_=0; i_<=n-1;i_++)
2489	{
2490	x[i_,k] = v*xc[i_];
2491	}
2492	}
2493	}
2494
2495
2496	/*************************************************************************
2497	Internal subroutine.
2498	Returns maximum count of RFS iterations as function of:
2499	1. machine epsilon
2500	2. task size.
2501	3. condition number
2502
2503	-- ALGLIB --
2504	Copyright 27.01.2010 by Bochkanov Sergey
2505	*************************************************************************/
2506	private static int densesolverrfsmax(int n,
2507	double r1,
2508	double rinf)
2509	{
2510	int result = 0;
2511
2512	result = 5;
2513	return result;
2514	}
2515
2516
2517	/*************************************************************************
2518	Internal subroutine.
2519	Returns maximum count of RFS iterations as function of:
2520	1. machine epsilon
2521	2. task size.
2522	3. norm-2 condition number
2523
2524	-- ALGLIB --
2525	Copyright 27.01.2010 by Bochkanov Sergey
2526	*************************************************************************/
2527	private static int densesolverrfsmaxv2(int n,
2528	double r2)
2529	{
2530	int result = 0;
2531
2532	result = densesolverrfsmax(n, 0, 0);
2533	return result;
2534	}
2535
2536
2537	/*************************************************************************
2538	Basic LU solver for ScaleAPLUx = y.
2539
2540	This subroutine assumes that:
2541	* L is well-scaled, and it is U which needs scaling by ScaleA.
2542	* A=PLU is well-conditioned, so no zero divisions or overflow may occur
2543
2544	-- ALGLIB --
2545	Copyright 27.01.2010 by Bochkanov Sergey
2546	*************************************************************************/
2547	private static void rbasiclusolve(ref double[,] lua,
2548	ref int[] p,
2549	double scalea,
2550	int n,
2551	ref double[] xb,
2552	ref double[] tmp)
2553	{
2554	int i = 0;
2555	double v = 0;
2556	int i_ = 0;
2557
2558	for(i=0; i<=n-1; i++)
2559	{
2560	if( p[i]!=i )
2561	{
2562	v = xb[i];
2563	xb[i] = xb[p[i]];
2564	xb[p[i]] = v;
2565	}
2566	}
2567	for(i=1; i<=n-1; i++)
2568	{
2569	v = 0.0;
2570	for(i_=0; i_<=i-1;i_++)
2571	{
2572	v += lua[i,i_]*xb[i_];
2573	}
2574	xb[i] = xb[i]-v;
2575	}
2576	xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
2577	for(i=n-2; i>=0; i--)
2578	{
2579	for(i_=i+1; i_<=n-1;i_++)
2580	{
2581	tmp[i_] = scalea*lua[i,i_];
2582	}
2583	v = 0.0;
2584	for(i_=i+1; i_<=n-1;i_++)
2585	{
2586	v += tmp[i_]*xb[i_];
2587	}
2588	xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
2589	}
2590	}
2591
2592
2593	/*************************************************************************
2594	Basic Cholesky solver for ScaleACholesky(A)'x = y.
2595
2596	This subroutine assumes that:
2597	* A*ScaleA is well scaled
2598	* A is well-conditioned, so no zero divisions or overflow may occur
2599
2600	-- ALGLIB --
2601	Copyright 27.01.2010 by Bochkanov Sergey
2602	*************************************************************************/
2603	private static void spdbasiccholeskysolve(ref double[,] cha,
2604	double sqrtscalea,
2605	int n,
2606	bool isupper,
2607	ref double[] xb,
2608	ref double[] tmp)
2609	{
2610	int i = 0;
2611	double v = 0;
2612	int i_ = 0;
2613
2614
2615	//
2616	// A = LL' or A=U'U
2617	//
2618	if( isupper )
2619	{
2620
2621	//
2622	// Solve U'*y=b first.
2623	//
2624	for(i=0; i<=n-1; i++)
2625	{
2626	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2627	if( i<n-1 )
2628	{
2629	v = xb[i];
2630	for(i_=i+1; i_<=n-1;i_++)
2631	{
2632	tmp[i_] = sqrtscalea*cha[i,i_];
2633	}
2634	for(i_=i+1; i_<=n-1;i_++)
2635	{
2636	xb[i_] = xb[i_] - v*tmp[i_];
2637	}
2638	}
2639	}
2640
2641	//
2642	// Solve U*x=y then.
2643	//
2644	for(i=n-1; i>=0; i--)
2645	{
2646	if( i<n-1 )
2647	{
2648	for(i_=i+1; i_<=n-1;i_++)
2649	{
2650	tmp[i_] = sqrtscalea*cha[i,i_];
2651	}
2652	v = 0.0;
2653	for(i_=i+1; i_<=n-1;i_++)
2654	{
2655	v += tmp[i_]*xb[i_];
2656	}
2657	xb[i] = xb[i]-v;
2658	}
2659	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2660	}
2661	}
2662	else
2663	{
2664
2665	//
2666	// Solve L*y=b first
2667	//
2668	for(i=0; i<=n-1; i++)
2669	{
2670	if( i>0 )
2671	{
2672	for(i_=0; i_<=i-1;i_++)
2673	{
2674	tmp[i_] = sqrtscalea*cha[i,i_];
2675	}
2676	v = 0.0;
2677	for(i_=0; i_<=i-1;i_++)
2678	{
2679	v += tmp[i_]*xb[i_];
2680	}
2681	xb[i] = xb[i]-v;
2682	}
2683	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2684	}
2685
2686	//
2687	// Solve L'*x=y then.
2688	//
2689	for(i=n-1; i>=0; i--)
2690	{
2691	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2692	if( i>0 )
2693	{
2694	v = xb[i];
2695	for(i_=0; i_<=i-1;i_++)
2696	{
2697	tmp[i_] = sqrtscalea*cha[i,i_];
2698	}
2699	for(i_=0; i_<=i-1;i_++)
2700	{
2701	xb[i_] = xb[i_] - v*tmp[i_];
2702	}
2703	}
2704	}
2705	}
2706	}
2707
2708
2709	/*************************************************************************
2710	Basic LU solver for ScaleAPLUx = y.
2711
2712	This subroutine assumes that:
2713	* L is well-scaled, and it is U which needs scaling by ScaleA.
2714	* A=PLU is well-conditioned, so no zero divisions or overflow may occur
2715
2716	-- ALGLIB --
2717	Copyright 27.01.2010 by Bochkanov Sergey
2718	*************************************************************************/
2719	private static void cbasiclusolve(ref AP.Complex[,] lua,
2720	ref int[] p,
2721	double scalea,
2722	int n,
2723	ref AP.Complex[] xb,
2724	ref AP.Complex[] tmp)
2725	{
2726	int i = 0;
2727	AP.Complex v = 0;
2728	int i_ = 0;
2729
2730	for(i=0; i<=n-1; i++)
2731	{
2732	if( p[i]!=i )
2733	{
2734	v = xb[i];
2735	xb[i] = xb[p[i]];
2736	xb[p[i]] = v;
2737	}
2738	}
2739	for(i=1; i<=n-1; i++)
2740	{
2741	v = 0.0;
2742	for(i_=0; i_<=i-1;i_++)
2743	{
2744	v += lua[i,i_]*xb[i_];
2745	}
2746	xb[i] = xb[i]-v;
2747	}
2748	xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
2749	for(i=n-2; i>=0; i--)
2750	{
2751	for(i_=i+1; i_<=n-1;i_++)
2752	{
2753	tmp[i_] = scalea*lua[i,i_];
2754	}
2755	v = 0.0;
2756	for(i_=i+1; i_<=n-1;i_++)
2757	{
2758	v += tmp[i_]*xb[i_];
2759	}
2760	xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
2761	}
2762	}
2763
2764
2765	/*************************************************************************
2766	Basic Cholesky solver for ScaleACholesky(A)'x = y.
2767
2768	This subroutine assumes that:
2769	* A*ScaleA is well scaled
2770	* A is well-conditioned, so no zero divisions or overflow may occur
2771
2772	-- ALGLIB --
2773	Copyright 27.01.2010 by Bochkanov Sergey
2774	*************************************************************************/
2775	private static void hpdbasiccholeskysolve(ref AP.Complex[,] cha,
2776	double sqrtscalea,
2777	int n,
2778	bool isupper,
2779	ref AP.Complex[] xb,
2780	ref AP.Complex[] tmp)
2781	{
2782	int i = 0;
2783	AP.Complex v = 0;
2784	int i_ = 0;
2785
2786
2787	//
2788	// A = LL' or A=U'U
2789	//
2790	if( isupper )
2791	{
2792
2793	//
2794	// Solve U'*y=b first.
2795	//
2796	for(i=0; i<=n-1; i++)
2797	{
2798	xb[i] = xb[i]/(sqrtscalea*AP.Math.Conj(cha[i,i]));
2799	if( i<n-1 )
2800	{
2801	v = xb[i];
2802	for(i_=i+1; i_<=n-1;i_++)
2803	{
2804	tmp[i_] = sqrtscalea*AP.Math.Conj(cha[i,i_]);
2805	}
2806	for(i_=i+1; i_<=n-1;i_++)
2807	{
2808	xb[i_] = xb[i_] - v*tmp[i_];
2809	}
2810	}
2811	}
2812
2813	//
2814	// Solve U*x=y then.
2815	//
2816	for(i=n-1; i>=0; i--)
2817	{
2818	if( i<n-1 )
2819	{
2820	for(i_=i+1; i_<=n-1;i_++)
2821	{
2822	tmp[i_] = sqrtscalea*cha[i,i_];
2823	}
2824	v = 0.0;
2825	for(i_=i+1; i_<=n-1;i_++)
2826	{
2827	v += tmp[i_]*xb[i_];
2828	}
2829	xb[i] = xb[i]-v;
2830	}
2831	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2832	}
2833	}
2834	else
2835	{
2836
2837	//
2838	// Solve L*y=b first
2839	//
2840	for(i=0; i<=n-1; i++)
2841	{
2842	if( i>0 )
2843	{
2844	for(i_=0; i_<=i-1;i_++)
2845	{
2846	tmp[i_] = sqrtscalea*cha[i,i_];
2847	}
2848	v = 0.0;
2849	for(i_=0; i_<=i-1;i_++)
2850	{
2851	v += tmp[i_]*xb[i_];
2852	}
2853	xb[i] = xb[i]-v;
2854	}
2855	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2856	}
2857
2858	//
2859	// Solve L'*x=y then.
2860	//
2861	for(i=n-1; i>=0; i--)
2862	{
2863	xb[i] = xb[i]/(sqrtscalea*AP.Math.Conj(cha[i,i]));
2864	if( i>0 )
2865	{
2866	v = xb[i];
2867	for(i_=0; i_<=i-1;i_++)
2868	{
2869	tmp[i_] = sqrtscalea*AP.Math.Conj(cha[i,i_]);
2870	}
2871	for(i_=0; i_<=i-1;i_++)
2872	{
2873	xb[i_] = xb[i_] - v*tmp[i_];
2874	}
2875	}
2876	}
2877	}
2878	}
2879	}
2880	}

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

Download in other formats:

Update cookies preferences