Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/rcond.cs @ 3932

Visit:

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

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

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

Download in other formats:

Update cookies preferences