Context Navigation

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/densesolver.cs @ 17712

Visit:

Last change on this file since 17712 was 2806, checked in by gkronber, 15 years ago
Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 94.4 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	x = new double[n, m];
1888	y = new double[n];
1889	xc = new double[n];
1890	bc = new double[n];
1891	tx = new double[n+1];
1892	xa = new double[n+1];
1893	xb = new double[n+1];
1894
1895	//
1896	// estimate condition number, test for near singularity
1897	//
1898	rep.r1 = rcond.rmatrixlurcond1(ref lua, n);
1899	rep.rinf = rcond.rmatrixlurcondinf(ref lua, n);
1900	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
1901	{
1902	for(i=0; i<=n-1; i++)
1903	{
1904	for(j=0; j<=m-1; j++)
1905	{
1906	x[i,j] = 0;
1907	}
1908	}
1909	rep.r1 = 0;
1910	rep.rinf = 0;
1911	info = -3;
1912	return;
1913	}
1914	info = 1;
1915
1916	//
1917	// solve
1918	//
1919	for(k=0; k<=m-1; k++)
1920	{
1921
1922	//
1923	// copy B to contiguous storage
1924	//
1925	for(i_=0; i_<=n-1;i_++)
1926	{
1927	bc[i_] = b[i_,k];
1928	}
1929
1930	//
1931	// Scale right part:
1932	// * MX stores max(\|Bi\|)
1933	// * ScaleRight stores actual scaling applied to B when solving systems
1934	// it is chosen to make \|scaleRight*b\| close to 1.
1935	//
1936	mxb = 0;
1937	for(i=0; i<=n-1; i++)
1938	{
1939	mxb = Math.Max(mxb, Math.Abs(bc[i]));
1940	}
1941	if( (double)(mxb)==(double)(0) )
1942	{
1943	mxb = 1;
1944	}
1945	scaleright = 1/mxb;
1946
1947	//
1948	// First, non-iterative part of solution process.
1949	// We use separate code for this task because
1950	// XDot is quite slow and we want to save time.
1951	//
1952	for(i_=0; i_<=n-1;i_++)
1953	{
1954	xc[i_] = scaleright*bc[i_];
1955	}
1956	rbasiclusolve(ref lua, ref p, scalea, n, ref xc, ref tx);
1957
1958	//
1959	// Iterative refinement of xc:
1960	// * calculate r = bc-A*xc using extra-precise dot product
1961	// * solve A*y = r
1962	// * update x:=x+r
1963	//
1964	// This cycle is executed until one of two things happens:
1965	// 1. maximum number of iterations reached
1966	// 2. last iteration decreased error to the lower limit
1967	//
1968	if( havea )
1969	{
1970	nrfs = densesolverrfsmax(n, rep.r1, rep.rinf);
1971	terminatenexttime = false;
1972	for(rfs=0; rfs<=nrfs-1; rfs++)
1973	{
1974	if( terminatenexttime )
1975	{
1976	break;
1977	}
1978
1979	//
1980	// generate right part
1981	//
1982	smallerr = true;
1983	for(i_=0; i_<=n-1;i_++)
1984	{
1985	xb[i_] = xc[i_];
1986	}
1987	for(i=0; i<=n-1; i++)
1988	{
1989	for(i_=0; i_<=n-1;i_++)
1990	{
1991	xa[i_] = scalea*a[i,i_];
1992	}
1993	xa[n] = -1;
1994	xb[n] = scaleright*bc[i];
1995	xblas.xdot(ref xa, ref xb, n+1, ref tx, ref v, ref verr);
1996	y[i] = -v;
1997	smallerr = smallerr & (double)(Math.Abs(v))<(double)(4*verr);
1998	}
1999	if( smallerr )
2000	{
2001	terminatenexttime = true;
2002	}
2003
2004	//
2005	// solve and update
2006	//
2007	rbasiclusolve(ref lua, ref p, scalea, n, ref y, ref tx);
2008	for(i_=0; i_<=n-1;i_++)
2009	{
2010	xc[i_] = xc[i_] + y[i_];
2011	}
2012	}
2013	}
2014
2015	//
2016	// Store xc.
2017	// Post-scale result.
2018	//
2019	v = scalea*mxb;
2020	for(i_=0; i_<=n-1;i_++)
2021	{
2022	x[i_,k] = v*xc[i_];
2023	}
2024	}
2025	}
2026
2027
2028	/*************************************************************************
2029	Internal Cholesky solver
2030
2031	-- ALGLIB --
2032	Copyright 27.01.2010 by Bochkanov Sergey
2033	*************************************************************************/
2034	private static void spdmatrixcholeskysolveinternal(ref double[,] cha,
2035	double sqrtscalea,
2036	int n,
2037	bool isupper,
2038	ref double[,] a,
2039	bool havea,
2040	ref double[,] b,
2041	int m,
2042	ref int info,
2043	ref densesolverreport rep,
2044	ref double[,] x)
2045	{
2046	int i = 0;
2047	int j = 0;
2048	int k = 0;
2049	int rfs = 0;
2050	int nrfs = 0;
2051	double[] xc = new double[0];
2052	double[] y = new double[0];
2053	double[] bc = new double[0];
2054	double[] xa = new double[0];
2055	double[] xb = new double[0];
2056	double[] tx = new double[0];
2057	double v = 0;
2058	double verr = 0;
2059	double mxb = 0;
2060	double scaleright = 0;
2061	bool smallerr = new bool();
2062	bool terminatenexttime = new bool();
2063	int i_ = 0;
2064
2065	System.Diagnostics.Debug.Assert((double)(sqrtscalea)>(double)(0));
2066
2067	//
2068	// prepare: check inputs, allocate space...
2069	//
2070	if( n<=0 \| m<=0 )
2071	{
2072	info = -1;
2073	return;
2074	}
2075	x = new double[n, m];
2076	y = new double[n];
2077	xc = new double[n];
2078	bc = new double[n];
2079	tx = new double[n+1];
2080	xa = new double[n+1];
2081	xb = new double[n+1];
2082
2083	//
2084	// estimate condition number, test for near singularity
2085	//
2086	rep.r1 = rcond.spdmatrixcholeskyrcond(ref cha, n, isupper);
2087	rep.rinf = rep.r1;
2088	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) )
2089	{
2090	for(i=0; i<=n-1; i++)
2091	{
2092	for(j=0; j<=m-1; j++)
2093	{
2094	x[i,j] = 0;
2095	}
2096	}
2097	rep.r1 = 0;
2098	rep.rinf = 0;
2099	info = -3;
2100	return;
2101	}
2102	info = 1;
2103
2104	//
2105	// solve
2106	//
2107	for(k=0; k<=m-1; k++)
2108	{
2109
2110	//
2111	// copy B to contiguous storage
2112	//
2113	for(i_=0; i_<=n-1;i_++)
2114	{
2115	bc[i_] = b[i_,k];
2116	}
2117
2118	//
2119	// Scale right part:
2120	// * MX stores max(\|Bi\|)
2121	// * ScaleRight stores actual scaling applied to B when solving systems
2122	// it is chosen to make \|scaleRight*b\| close to 1.
2123	//
2124	mxb = 0;
2125	for(i=0; i<=n-1; i++)
2126	{
2127	mxb = Math.Max(mxb, Math.Abs(bc[i]));
2128	}
2129	if( (double)(mxb)==(double)(0) )
2130	{
2131	mxb = 1;
2132	}
2133	scaleright = 1/mxb;
2134
2135	//
2136	// First, non-iterative part of solution process.
2137	// We use separate code for this task because
2138	// XDot is quite slow and we want to save time.
2139	//
2140	for(i_=0; i_<=n-1;i_++)
2141	{
2142	xc[i_] = scaleright*bc[i_];
2143	}
2144	spdbasiccholeskysolve(ref cha, sqrtscalea, n, isupper, ref xc, ref tx);
2145
2146	//
2147	// Store xc.
2148	// Post-scale result.
2149	//
2150	v = AP.Math.Sqr(sqrtscalea)*mxb;
2151	for(i_=0; i_<=n-1;i_++)
2152	{
2153	x[i_,k] = v*xc[i_];
2154	}
2155	}
2156	}
2157
2158
2159	/*************************************************************************
2160	Internal LU solver
2161
2162	-- ALGLIB --
2163	Copyright 27.01.2010 by Bochkanov Sergey
2164	*************************************************************************/
2165	private static void cmatrixlusolveinternal(ref AP.Complex[,] lua,
2166	ref int[] p,
2167	double scalea,
2168	int n,
2169	ref AP.Complex[,] a,
2170	bool havea,
2171	ref AP.Complex[,] b,
2172	int m,
2173	ref int info,
2174	ref densesolverreport rep,
2175	ref AP.Complex[,] x)
2176	{
2177	int i = 0;
2178	int j = 0;
2179	int k = 0;
2180	int rfs = 0;
2181	int nrfs = 0;
2182	AP.Complex[] xc = new AP.Complex[0];
2183	AP.Complex[] y = new AP.Complex[0];
2184	AP.Complex[] bc = new AP.Complex[0];
2185	AP.Complex[] xa = new AP.Complex[0];
2186	AP.Complex[] xb = new AP.Complex[0];
2187	AP.Complex[] tx = new AP.Complex[0];
2188	double[] tmpbuf = new double[0];
2189	AP.Complex v = 0;
2190	double verr = 0;
2191	double mxb = 0;
2192	double scaleright = 0;
2193	bool smallerr = new bool();
2194	bool terminatenexttime = new bool();
2195	int i_ = 0;
2196
2197	System.Diagnostics.Debug.Assert((double)(scalea)>(double)(0));
2198
2199	//
2200	// prepare: check inputs, allocate space...
2201	//
2202	if( n<=0 \| m<=0 )
2203	{
2204	info = -1;
2205	return;
2206	}
2207	x = new AP.Complex[n, m];
2208	y = new AP.Complex[n];
2209	xc = new AP.Complex[n];
2210	bc = new AP.Complex[n];
2211	tx = new AP.Complex[n];
2212	xa = new AP.Complex[n+1];
2213	xb = new AP.Complex[n+1];
2214	tmpbuf = new double[2*n+2];
2215
2216	//
2217	// estimate condition number, test for near singularity
2218	//
2219	rep.r1 = rcond.cmatrixlurcond1(ref lua, n);
2220	rep.rinf = rcond.cmatrixlurcondinf(ref lua, n);
2221	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) \| (double)(rep.rinf)<(double)(rcond.rcondthreshold()) )
2222	{
2223	for(i=0; i<=n-1; i++)
2224	{
2225	for(j=0; j<=m-1; j++)
2226	{
2227	x[i,j] = 0;
2228	}
2229	}
2230	rep.r1 = 0;
2231	rep.rinf = 0;
2232	info = -3;
2233	return;
2234	}
2235	info = 1;
2236
2237	//
2238	// solve
2239	//
2240	for(k=0; k<=m-1; k++)
2241	{
2242
2243	//
2244	// copy B to contiguous storage
2245	//
2246	for(i_=0; i_<=n-1;i_++)
2247	{
2248	bc[i_] = b[i_,k];
2249	}
2250
2251	//
2252	// Scale right part:
2253	// * MX stores max(\|Bi\|)
2254	// * ScaleRight stores actual scaling applied to B when solving systems
2255	// it is chosen to make \|scaleRight*b\| close to 1.
2256	//
2257	mxb = 0;
2258	for(i=0; i<=n-1; i++)
2259	{
2260	mxb = Math.Max(mxb, AP.Math.AbsComplex(bc[i]));
2261	}
2262	if( (double)(mxb)==(double)(0) )
2263	{
2264	mxb = 1;
2265	}
2266	scaleright = 1/mxb;
2267
2268	//
2269	// First, non-iterative part of solution process.
2270	// We use separate code for this task because
2271	// XDot is quite slow and we want to save time.
2272	//
2273	for(i_=0; i_<=n-1;i_++)
2274	{
2275	xc[i_] = scaleright*bc[i_];
2276	}
2277	cbasiclusolve(ref lua, ref p, scalea, n, ref xc, ref tx);
2278
2279	//
2280	// Iterative refinement of xc:
2281	// * calculate r = bc-A*xc using extra-precise dot product
2282	// * solve A*y = r
2283	// * update x:=x+r
2284	//
2285	// This cycle is executed until one of two things happens:
2286	// 1. maximum number of iterations reached
2287	// 2. last iteration decreased error to the lower limit
2288	//
2289	if( havea )
2290	{
2291	nrfs = densesolverrfsmax(n, rep.r1, rep.rinf);
2292	terminatenexttime = false;
2293	for(rfs=0; rfs<=nrfs-1; rfs++)
2294	{
2295	if( terminatenexttime )
2296	{
2297	break;
2298	}
2299
2300	//
2301	// generate right part
2302	//
2303	smallerr = true;
2304	for(i_=0; i_<=n-1;i_++)
2305	{
2306	xb[i_] = xc[i_];
2307	}
2308	for(i=0; i<=n-1; i++)
2309	{
2310	for(i_=0; i_<=n-1;i_++)
2311	{
2312	xa[i_] = scalea*a[i,i_];
2313	}
2314	xa[n] = -1;
2315	xb[n] = scaleright*bc[i];
2316	xblas.xcdot(ref xa, ref xb, n+1, ref tmpbuf, ref v, ref verr);
2317	y[i] = -v;
2318	smallerr = smallerr & (double)(AP.Math.AbsComplex(v))<(double)(4*verr);
2319	}
2320	if( smallerr )
2321	{
2322	terminatenexttime = true;
2323	}
2324
2325	//
2326	// solve and update
2327	//
2328	cbasiclusolve(ref lua, ref p, scalea, n, ref y, ref tx);
2329	for(i_=0; i_<=n-1;i_++)
2330	{
2331	xc[i_] = xc[i_] + y[i_];
2332	}
2333	}
2334	}
2335
2336	//
2337	// Store xc.
2338	// Post-scale result.
2339	//
2340	v = scalea*mxb;
2341	for(i_=0; i_<=n-1;i_++)
2342	{
2343	x[i_,k] = v*xc[i_];
2344	}
2345	}
2346	}
2347
2348
2349	/*************************************************************************
2350	Internal Cholesky solver
2351
2352	-- ALGLIB --
2353	Copyright 27.01.2010 by Bochkanov Sergey
2354	*************************************************************************/
2355	private static void hpdmatrixcholeskysolveinternal(ref AP.Complex[,] cha,
2356	double sqrtscalea,
2357	int n,
2358	bool isupper,
2359	ref AP.Complex[,] a,
2360	bool havea,
2361	ref AP.Complex[,] b,
2362	int m,
2363	ref int info,
2364	ref densesolverreport rep,
2365	ref AP.Complex[,] x)
2366	{
2367	int i = 0;
2368	int j = 0;
2369	int k = 0;
2370	int rfs = 0;
2371	int nrfs = 0;
2372	AP.Complex[] xc = new AP.Complex[0];
2373	AP.Complex[] y = new AP.Complex[0];
2374	AP.Complex[] bc = new AP.Complex[0];
2375	AP.Complex[] xa = new AP.Complex[0];
2376	AP.Complex[] xb = new AP.Complex[0];
2377	AP.Complex[] tx = new AP.Complex[0];
2378	double v = 0;
2379	double verr = 0;
2380	double mxb = 0;
2381	double scaleright = 0;
2382	bool smallerr = new bool();
2383	bool terminatenexttime = new bool();
2384	int i_ = 0;
2385
2386	System.Diagnostics.Debug.Assert((double)(sqrtscalea)>(double)(0));
2387
2388	//
2389	// prepare: check inputs, allocate space...
2390	//
2391	if( n<=0 \| m<=0 )
2392	{
2393	info = -1;
2394	return;
2395	}
2396	x = new AP.Complex[n, m];
2397	y = new AP.Complex[n];
2398	xc = new AP.Complex[n];
2399	bc = new AP.Complex[n];
2400	tx = new AP.Complex[n+1];
2401	xa = new AP.Complex[n+1];
2402	xb = new AP.Complex[n+1];
2403
2404	//
2405	// estimate condition number, test for near singularity
2406	//
2407	rep.r1 = rcond.hpdmatrixcholeskyrcond(ref cha, n, isupper);
2408	rep.rinf = rep.r1;
2409	if( (double)(rep.r1)<(double)(rcond.rcondthreshold()) )
2410	{
2411	for(i=0; i<=n-1; i++)
2412	{
2413	for(j=0; j<=m-1; j++)
2414	{
2415	x[i,j] = 0;
2416	}
2417	}
2418	rep.r1 = 0;
2419	rep.rinf = 0;
2420	info = -3;
2421	return;
2422	}
2423	info = 1;
2424
2425	//
2426	// solve
2427	//
2428	for(k=0; k<=m-1; k++)
2429	{
2430
2431	//
2432	// copy B to contiguous storage
2433	//
2434	for(i_=0; i_<=n-1;i_++)
2435	{
2436	bc[i_] = b[i_,k];
2437	}
2438
2439	//
2440	// Scale right part:
2441	// * MX stores max(\|Bi\|)
2442	// * ScaleRight stores actual scaling applied to B when solving systems
2443	// it is chosen to make \|scaleRight*b\| close to 1.
2444	//
2445	mxb = 0;
2446	for(i=0; i<=n-1; i++)
2447	{
2448	mxb = Math.Max(mxb, AP.Math.AbsComplex(bc[i]));
2449	}
2450	if( (double)(mxb)==(double)(0) )
2451	{
2452	mxb = 1;
2453	}
2454	scaleright = 1/mxb;
2455
2456	//
2457	// First, non-iterative part of solution process.
2458	// We use separate code for this task because
2459	// XDot is quite slow and we want to save time.
2460	//
2461	for(i_=0; i_<=n-1;i_++)
2462	{
2463	xc[i_] = scaleright*bc[i_];
2464	}
2465	hpdbasiccholeskysolve(ref cha, sqrtscalea, n, isupper, ref xc, ref tx);
2466
2467	//
2468	// Store xc.
2469	// Post-scale result.
2470	//
2471	v = AP.Math.Sqr(sqrtscalea)*mxb;
2472	for(i_=0; i_<=n-1;i_++)
2473	{
2474	x[i_,k] = v*xc[i_];
2475	}
2476	}
2477	}
2478
2479
2480	/*************************************************************************
2481	Internal subroutine.
2482	Returns maximum count of RFS iterations as function of:
2483	1. machine epsilon
2484	2. task size.
2485	3. condition number
2486
2487	-- ALGLIB --
2488	Copyright 27.01.2010 by Bochkanov Sergey
2489	*************************************************************************/
2490	private static int densesolverrfsmax(int n,
2491	double r1,
2492	double rinf)
2493	{
2494	int result = 0;
2495
2496	result = 5;
2497	return result;
2498	}
2499
2500
2501	/*************************************************************************
2502	Internal subroutine.
2503	Returns maximum count of RFS iterations as function of:
2504	1. machine epsilon
2505	2. task size.
2506	3. norm-2 condition number
2507
2508	-- ALGLIB --
2509	Copyright 27.01.2010 by Bochkanov Sergey
2510	*************************************************************************/
2511	private static int densesolverrfsmaxv2(int n,
2512	double r2)
2513	{
2514	int result = 0;
2515
2516	result = densesolverrfsmax(n, 0, 0);
2517	return result;
2518	}
2519
2520
2521	/*************************************************************************
2522	Basic LU solver for ScaleAPLUx = y.
2523
2524	This subroutine assumes that:
2525	* L is well-scaled, and it is U which needs scaling by ScaleA.
2526	* A=PLU is well-conditioned, so no zero divisions or overflow may occur
2527
2528	-- ALGLIB --
2529	Copyright 27.01.2010 by Bochkanov Sergey
2530	*************************************************************************/
2531	private static void rbasiclusolve(ref double[,] lua,
2532	ref int[] p,
2533	double scalea,
2534	int n,
2535	ref double[] xb,
2536	ref double[] tmp)
2537	{
2538	int i = 0;
2539	double v = 0;
2540	int i_ = 0;
2541
2542	for(i=0; i<=n-1; i++)
2543	{
2544	if( p[i]!=i )
2545	{
2546	v = xb[i];
2547	xb[i] = xb[p[i]];
2548	xb[p[i]] = v;
2549	}
2550	}
2551	for(i=1; i<=n-1; i++)
2552	{
2553	v = 0.0;
2554	for(i_=0; i_<=i-1;i_++)
2555	{
2556	v += lua[i,i_]*xb[i_];
2557	}
2558	xb[i] = xb[i]-v;
2559	}
2560	xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
2561	for(i=n-2; i>=0; i--)
2562	{
2563	for(i_=i+1; i_<=n-1;i_++)
2564	{
2565	tmp[i_] = scalea*lua[i,i_];
2566	}
2567	v = 0.0;
2568	for(i_=i+1; i_<=n-1;i_++)
2569	{
2570	v += tmp[i_]*xb[i_];
2571	}
2572	xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
2573	}
2574	}
2575
2576
2577	/*************************************************************************
2578	Basic Cholesky solver for ScaleACholesky(A)'x = y.
2579
2580	This subroutine assumes that:
2581	* A*ScaleA is well scaled
2582	* A is well-conditioned, so no zero divisions or overflow may occur
2583
2584	-- ALGLIB --
2585	Copyright 27.01.2010 by Bochkanov Sergey
2586	*************************************************************************/
2587	private static void spdbasiccholeskysolve(ref double[,] cha,
2588	double sqrtscalea,
2589	int n,
2590	bool isupper,
2591	ref double[] xb,
2592	ref double[] tmp)
2593	{
2594	int i = 0;
2595	double v = 0;
2596	int i_ = 0;
2597
2598
2599	//
2600	// A = LL' or A=U'U
2601	//
2602	if( isupper )
2603	{
2604
2605	//
2606	// Solve U'*y=b first.
2607	//
2608	for(i=0; i<=n-1; i++)
2609	{
2610	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2611	if( i<n-1 )
2612	{
2613	v = xb[i];
2614	for(i_=i+1; i_<=n-1;i_++)
2615	{
2616	tmp[i_] = sqrtscalea*cha[i,i_];
2617	}
2618	for(i_=i+1; i_<=n-1;i_++)
2619	{
2620	xb[i_] = xb[i_] - v*tmp[i_];
2621	}
2622	}
2623	}
2624
2625	//
2626	// Solve U*x=y then.
2627	//
2628	for(i=n-1; i>=0; i--)
2629	{
2630	if( i<n-1 )
2631	{
2632	for(i_=i+1; i_<=n-1;i_++)
2633	{
2634	tmp[i_] = sqrtscalea*cha[i,i_];
2635	}
2636	v = 0.0;
2637	for(i_=i+1; i_<=n-1;i_++)
2638	{
2639	v += tmp[i_]*xb[i_];
2640	}
2641	xb[i] = xb[i]-v;
2642	}
2643	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2644	}
2645	}
2646	else
2647	{
2648
2649	//
2650	// Solve L*y=b first
2651	//
2652	for(i=0; i<=n-1; i++)
2653	{
2654	if( i>0 )
2655	{
2656	for(i_=0; i_<=i-1;i_++)
2657	{
2658	tmp[i_] = sqrtscalea*cha[i,i_];
2659	}
2660	v = 0.0;
2661	for(i_=0; i_<=i-1;i_++)
2662	{
2663	v += tmp[i_]*xb[i_];
2664	}
2665	xb[i] = xb[i]-v;
2666	}
2667	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2668	}
2669
2670	//
2671	// Solve L'*x=y then.
2672	//
2673	for(i=n-1; i>=0; i--)
2674	{
2675	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2676	if( i>0 )
2677	{
2678	v = xb[i];
2679	for(i_=0; i_<=i-1;i_++)
2680	{
2681	tmp[i_] = sqrtscalea*cha[i,i_];
2682	}
2683	for(i_=0; i_<=i-1;i_++)
2684	{
2685	xb[i_] = xb[i_] - v*tmp[i_];
2686	}
2687	}
2688	}
2689	}
2690	}
2691
2692
2693	/*************************************************************************
2694	Basic LU solver for ScaleAPLUx = y.
2695
2696	This subroutine assumes that:
2697	* L is well-scaled, and it is U which needs scaling by ScaleA.
2698	* A=PLU is well-conditioned, so no zero divisions or overflow may occur
2699
2700	-- ALGLIB --
2701	Copyright 27.01.2010 by Bochkanov Sergey
2702	*************************************************************************/
2703	private static void cbasiclusolve(ref AP.Complex[,] lua,
2704	ref int[] p,
2705	double scalea,
2706	int n,
2707	ref AP.Complex[] xb,
2708	ref AP.Complex[] tmp)
2709	{
2710	int i = 0;
2711	AP.Complex v = 0;
2712	int i_ = 0;
2713
2714	for(i=0; i<=n-1; i++)
2715	{
2716	if( p[i]!=i )
2717	{
2718	v = xb[i];
2719	xb[i] = xb[p[i]];
2720	xb[p[i]] = v;
2721	}
2722	}
2723	for(i=1; i<=n-1; i++)
2724	{
2725	v = 0.0;
2726	for(i_=0; i_<=i-1;i_++)
2727	{
2728	v += lua[i,i_]*xb[i_];
2729	}
2730	xb[i] = xb[i]-v;
2731	}
2732	xb[n-1] = xb[n-1]/(scalea*lua[n-1,n-1]);
2733	for(i=n-2; i>=0; i--)
2734	{
2735	for(i_=i+1; i_<=n-1;i_++)
2736	{
2737	tmp[i_] = scalea*lua[i,i_];
2738	}
2739	v = 0.0;
2740	for(i_=i+1; i_<=n-1;i_++)
2741	{
2742	v += tmp[i_]*xb[i_];
2743	}
2744	xb[i] = (xb[i]-v)/(scalea*lua[i,i]);
2745	}
2746	}
2747
2748
2749	/*************************************************************************
2750	Basic Cholesky solver for ScaleACholesky(A)'x = y.
2751
2752	This subroutine assumes that:
2753	* A*ScaleA is well scaled
2754	* A is well-conditioned, so no zero divisions or overflow may occur
2755
2756	-- ALGLIB --
2757	Copyright 27.01.2010 by Bochkanov Sergey
2758	*************************************************************************/
2759	private static void hpdbasiccholeskysolve(ref AP.Complex[,] cha,
2760	double sqrtscalea,
2761	int n,
2762	bool isupper,
2763	ref AP.Complex[] xb,
2764	ref AP.Complex[] tmp)
2765	{
2766	int i = 0;
2767	AP.Complex v = 0;
2768	int i_ = 0;
2769
2770
2771	//
2772	// A = LL' or A=U'U
2773	//
2774	if( isupper )
2775	{
2776
2777	//
2778	// Solve U'*y=b first.
2779	//
2780	for(i=0; i<=n-1; i++)
2781	{
2782	xb[i] = xb[i]/(sqrtscalea*AP.Math.Conj(cha[i,i]));
2783	if( i<n-1 )
2784	{
2785	v = xb[i];
2786	for(i_=i+1; i_<=n-1;i_++)
2787	{
2788	tmp[i_] = sqrtscalea*AP.Math.Conj(cha[i,i_]);
2789	}
2790	for(i_=i+1; i_<=n-1;i_++)
2791	{
2792	xb[i_] = xb[i_] - v*tmp[i_];
2793	}
2794	}
2795	}
2796
2797	//
2798	// Solve U*x=y then.
2799	//
2800	for(i=n-1; i>=0; i--)
2801	{
2802	if( i<n-1 )
2803	{
2804	for(i_=i+1; i_<=n-1;i_++)
2805	{
2806	tmp[i_] = sqrtscalea*cha[i,i_];
2807	}
2808	v = 0.0;
2809	for(i_=i+1; i_<=n-1;i_++)
2810	{
2811	v += tmp[i_]*xb[i_];
2812	}
2813	xb[i] = xb[i]-v;
2814	}
2815	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2816	}
2817	}
2818	else
2819	{
2820
2821	//
2822	// Solve L*y=b first
2823	//
2824	for(i=0; i<=n-1; i++)
2825	{
2826	if( i>0 )
2827	{
2828	for(i_=0; i_<=i-1;i_++)
2829	{
2830	tmp[i_] = sqrtscalea*cha[i,i_];
2831	}
2832	v = 0.0;
2833	for(i_=0; i_<=i-1;i_++)
2834	{
2835	v += tmp[i_]*xb[i_];
2836	}
2837	xb[i] = xb[i]-v;
2838	}
2839	xb[i] = xb[i]/(sqrtscalea*cha[i,i]);
2840	}
2841
2842	//
2843	// Solve L'*x=y then.
2844	//
2845	for(i=n-1; i>=0; i--)
2846	{
2847	xb[i] = xb[i]/(sqrtscalea*AP.Math.Conj(cha[i,i]));
2848	if( i>0 )
2849	{
2850	v = xb[i];
2851	for(i_=0; i_<=i-1;i_++)
2852	{
2853	tmp[i_] = sqrtscalea*AP.Math.Conj(cha[i,i_]);
2854	}
2855	for(i_=0; i_<=i-1;i_++)
2856	{
2857	xb[i_] = xb[i_] - v*tmp[i_];
2858	}
2859	}
2860	}
2861	}
2862	}
2863	}
2864	}

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

Download in other formats:

Update cookies preferences