Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/rcond.cs @ 2806

Last change on this file since 2806 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 75.0 KB
Line 
1/*************************************************************************
2Copyright (c) 1992-2007 The University of Tennessee.  All rights reserved.
3
4Contributors:
5    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
6      pseudocode.
7
8See subroutines comments for additional copyrights.
9
10>>> SOURCE LICENSE >>>
11This program is free software; you can redistribute it and/or modify
12it under the terms of the GNU General Public License as published by
13the Free Software Foundation (www.fsf.org); either version 2 of the
14License, or (at your option) any later version.
15
16This program is distributed in the hope that it will be useful,
17but WITHOUT ANY WARRANTY; without even the implied warranty of
18MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19GNU General Public License for more details.
20
21A copy of the GNU General Public License is available at
22http://www.fsf.org/licensing/licenses
23
24>>> END OF LICENSE >>>
25*************************************************************************/
26
27using System;
28
29namespace 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)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
48            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)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
103            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)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
160            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        Condition number estimate of a Hermitian positive definite matrix.
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        It should be noted that 1-norm and inf-norm of condition numbers of symmetric
233        matrices are equal, so the algorithm doesn't take into account the
234        differences between these types of norms.
235
236        Input parameters:
237            A       -   Hermitian positive definite matrix which is given by its
238                        upper or lower triangle depending on the value of
239                        IsUpper. Array with elements [0..N-1, 0..N-1].
240            N       -   size of matrix A.
241            IsUpper -   storage format.
242
243        Result:
244            1/LowerBound(cond(A)), if matrix A is positive definite,
245           -1, if matrix A is not positive definite, and its condition number
246            could not be found by this algorithm.
247
248        NOTE:
249            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
250            returned in such cases.
251        *************************************************************************/
252        public static double hpdmatrixrcond(AP.Complex[,] a,
253            int n,
254            bool isupper)
255        {
256            double result = 0;
257            int i = 0;
258            int j = 0;
259            int j1 = 0;
260            int j2 = 0;
261            double v = 0;
262            double nrm = 0;
263            double[] t = new double[0];
264
265            a = (AP.Complex[,])a.Clone();
266
267            t = new double[n];
268            for(i=0; i<=n-1; i++)
269            {
270                t[i] = 0;
271            }
272            for(i=0; i<=n-1; i++)
273            {
274                if( isupper )
275                {
276                    j1 = i;
277                    j2 = n-1;
278                }
279                else
280                {
281                    j1 = 0;
282                    j2 = i;
283                }
284                for(j=j1; j<=j2; j++)
285                {
286                    if( i==j )
287                    {
288                        t[i] = t[i]+AP.Math.AbsComplex(a[i,i]);
289                    }
290                    else
291                    {
292                        t[i] = t[i]+AP.Math.AbsComplex(a[i,j]);
293                        t[j] = t[j]+AP.Math.AbsComplex(a[i,j]);
294                    }
295                }
296            }
297            nrm = 0;
298            for(i=0; i<=n-1; i++)
299            {
300                nrm = Math.Max(nrm, t[i]);
301            }
302            if( trfac.hpdmatrixcholesky(ref a, n, isupper) )
303            {
304                hpdmatrixrcondcholeskyinternal(ref a, n, isupper, true, nrm, ref v);
305                result = v;
306            }
307            else
308            {
309                result = -1;
310            }
311            return result;
312        }
313
314
315        /*************************************************************************
316        Estimate of a matrix condition number (1-norm)
317
318        The algorithm calculates a lower bound of the condition number. In this case,
319        the algorithm does not return a lower bound of the condition number, but an
320        inverse number (to avoid an overflow in case of a singular matrix).
321
322        Input parameters:
323            A   -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
324            N   -   size of matrix A.
325
326        Result: 1/LowerBound(cond(A))
327
328        NOTE:
329            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
330            returned in such cases.
331        *************************************************************************/
332        public static double cmatrixrcond1(AP.Complex[,] a,
333            int n)
334        {
335            double result = 0;
336            int i = 0;
337            int j = 0;
338            double v = 0;
339            double nrm = 0;
340            int[] pivots = new int[0];
341            double[] t = new double[0];
342
343            a = (AP.Complex[,])a.Clone();
344
345            System.Diagnostics.Debug.Assert(n>=1, "CMatrixRCond1: N<1!");
346            t = new double[n];
347            for(i=0; i<=n-1; i++)
348            {
349                t[i] = 0;
350            }
351            for(i=0; i<=n-1; i++)
352            {
353                for(j=0; j<=n-1; j++)
354                {
355                    t[j] = t[j]+AP.Math.AbsComplex(a[i,j]);
356                }
357            }
358            nrm = 0;
359            for(i=0; i<=n-1; i++)
360            {
361                nrm = Math.Max(nrm, t[i]);
362            }
363            trfac.cmatrixlu(ref a, n, n, ref pivots);
364            cmatrixrcondluinternal(ref a, n, true, true, nrm, ref v);
365            result = v;
366            return result;
367        }
368
369
370        /*************************************************************************
371        Estimate of a matrix condition number (infinity-norm).
372
373        The algorithm calculates a lower bound of the condition number. In this case,
374        the algorithm does not return a lower bound of the condition number, but an
375        inverse number (to avoid an overflow in case of a singular matrix).
376
377        Input parameters:
378            A   -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
379            N   -   size of matrix A.
380
381        Result: 1/LowerBound(cond(A))
382
383        NOTE:
384            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
385            returned in such cases.
386        *************************************************************************/
387        public static double cmatrixrcondinf(AP.Complex[,] a,
388            int n)
389        {
390            double result = 0;
391            int i = 0;
392            int j = 0;
393            double v = 0;
394            double nrm = 0;
395            int[] pivots = new int[0];
396
397            a = (AP.Complex[,])a.Clone();
398
399            System.Diagnostics.Debug.Assert(n>=1, "CMatrixRCondInf: N<1!");
400            nrm = 0;
401            for(i=0; i<=n-1; i++)
402            {
403                v = 0;
404                for(j=0; j<=n-1; j++)
405                {
406                    v = v+AP.Math.AbsComplex(a[i,j]);
407                }
408                nrm = Math.Max(nrm, v);
409            }
410            trfac.cmatrixlu(ref a, n, n, ref pivots);
411            cmatrixrcondluinternal(ref a, n, false, true, nrm, ref v);
412            result = v;
413            return result;
414        }
415
416
417        /*************************************************************************
418        Estimate of the condition number of a matrix given by its LU decomposition (1-norm)
419
420        The algorithm calculates a lower bound of the condition number. In this case,
421        the algorithm does not return a lower bound of the condition number, but an
422        inverse number (to avoid an overflow in case of a singular matrix).
423
424        Input parameters:
425            LUA         -   LU decomposition of a matrix in compact form. Output of
426                            the RMatrixLU subroutine.
427            N           -   size of matrix A.
428
429        Result: 1/LowerBound(cond(A))
430
431        NOTE:
432            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
433            returned in such cases.
434        *************************************************************************/
435        public static double rmatrixlurcond1(ref double[,] lua,
436            int n)
437        {
438            double result = 0;
439            double v = 0;
440
441            rmatrixrcondluinternal(ref lua, n, true, false, 0, ref v);
442            result = v;
443            return result;
444        }
445
446
447        /*************************************************************************
448        Estimate of the condition number of a matrix given by its LU decomposition
449        (infinity norm).
450
451        The algorithm calculates a lower bound of the condition number. In this case,
452        the algorithm does not return a lower bound of the condition number, but an
453        inverse number (to avoid an overflow in case of a singular matrix).
454
455        Input parameters:
456            LUA     -   LU decomposition of a matrix in compact form. Output of
457                        the RMatrixLU subroutine.
458            N       -   size of matrix A.
459
460        Result: 1/LowerBound(cond(A))
461
462        NOTE:
463            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
464            returned in such cases.
465        *************************************************************************/
466        public static double rmatrixlurcondinf(ref double[,] lua,
467            int n)
468        {
469            double result = 0;
470            double v = 0;
471
472            rmatrixrcondluinternal(ref lua, n, false, false, 0, ref v);
473            result = v;
474            return result;
475        }
476
477
478        /*************************************************************************
479        Condition number estimate of a symmetric positive definite matrix given by
480        Cholesky decomposition.
481
482        The algorithm calculates a lower bound of the condition number. In this
483        case, the algorithm does not return a lower bound of the condition number,
484        but an inverse number (to avoid an overflow in case of a singular matrix).
485
486        It should be noted that 1-norm and inf-norm condition numbers of symmetric
487        matrices are equal, so the algorithm doesn't take into account the
488        differences between these types of norms.
489
490        Input parameters:
491            CD  - Cholesky decomposition of matrix A,
492                  output of SMatrixCholesky subroutine.
493            N   - size of matrix A.
494
495        Result: 1/LowerBound(cond(A))
496
497        NOTE:
498            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
499            returned in such cases.
500        *************************************************************************/
501        public static double spdmatrixcholeskyrcond(ref double[,] a,
502            int n,
503            bool isupper)
504        {
505            double result = 0;
506            double v = 0;
507
508            spdmatrixrcondcholeskyinternal(ref a, n, isupper, false, 0, ref v);
509            result = v;
510            return result;
511        }
512
513
514        /*************************************************************************
515        Condition number estimate of a Hermitian positive definite matrix given by
516        Cholesky decomposition.
517
518        The algorithm calculates a lower bound of the condition number. In this
519        case, the algorithm does not return a lower bound of the condition number,
520        but an inverse number (to avoid an overflow in case of a singular matrix).
521
522        It should be noted that 1-norm and inf-norm condition numbers of symmetric
523        matrices are equal, so the algorithm doesn't take into account the
524        differences between these types of norms.
525
526        Input parameters:
527            CD  - Cholesky decomposition of matrix A,
528                  output of SMatrixCholesky subroutine.
529            N   - size of matrix A.
530
531        Result: 1/LowerBound(cond(A))
532
533        NOTE:
534            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
535            returned in such cases.
536        *************************************************************************/
537        public static double hpdmatrixcholeskyrcond(ref AP.Complex[,] a,
538            int n,
539            bool isupper)
540        {
541            double result = 0;
542            double v = 0;
543
544            hpdmatrixrcondcholeskyinternal(ref a, n, isupper, false, 0, ref v);
545            result = v;
546            return result;
547        }
548
549
550        /*************************************************************************
551        Estimate of the condition number of a matrix given by its LU decomposition (1-norm)
552
553        The algorithm calculates a lower bound of the condition number. In this case,
554        the algorithm does not return a lower bound of the condition number, but an
555        inverse number (to avoid an overflow in case of a singular matrix).
556
557        Input parameters:
558            LUA         -   LU decomposition of a matrix in compact form. Output of
559                            the CMatrixLU subroutine.
560            N           -   size of matrix A.
561
562        Result: 1/LowerBound(cond(A))
563
564        NOTE:
565            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
566            returned in such cases.
567        *************************************************************************/
568        public static double cmatrixlurcond1(ref AP.Complex[,] lua,
569            int n)
570        {
571            double result = 0;
572            double v = 0;
573
574            System.Diagnostics.Debug.Assert(n>=1, "CMatrixLURCond1: N<1!");
575            cmatrixrcondluinternal(ref lua, n, true, false, 0.0, ref v);
576            result = v;
577            return result;
578        }
579
580
581        /*************************************************************************
582        Estimate of the condition number of a matrix given by its LU decomposition
583        (infinity norm).
584
585        The algorithm calculates a lower bound of the condition number. In this case,
586        the algorithm does not return a lower bound of the condition number, but an
587        inverse number (to avoid an overflow in case of a singular matrix).
588
589        Input parameters:
590            LUA     -   LU decomposition of a matrix in compact form. Output of
591                        the CMatrixLU subroutine.
592            N       -   size of matrix A.
593
594        Result: 1/LowerBound(cond(A))
595
596        NOTE:
597            if k(A)>1/epsilon, then matrix is assumed degenerate, k(A)=INF, 0.0 is
598            returned in such cases.
599        *************************************************************************/
600        public static double cmatrixlurcondinf(ref AP.Complex[,] lua,
601            int n)
602        {
603            double result = 0;
604            double v = 0;
605
606            System.Diagnostics.Debug.Assert(n>=1, "CMatrixLURCondInf: N<1!");
607            cmatrixrcondluinternal(ref lua, n, false, false, 0.0, ref v);
608            result = v;
609            return result;
610        }
611
612
613        /*************************************************************************
614        Threshold for rcond: matrices with condition number beyond this  threshold
615        are considered singular.
616
617        Threshold must be far enough from underflow, at least Sqr(Threshold)  must
618        be greater than underflow.
619        *************************************************************************/
620        public static double rcondthreshold()
621        {
622            double result = 0;
623
624            result = Math.Sqrt(Math.Sqrt(AP.Math.MinRealNumber));
625            return result;
626        }
627
628
629        /*************************************************************************
630        Internal subroutine for condition number estimation
631
632          -- LAPACK routine (version 3.0) --
633             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
634             Courant Institute, Argonne National Lab, and Rice University
635             February 29, 1992
636        *************************************************************************/
637        private static void spdmatrixrcondcholeskyinternal(ref double[,] cha,
638            int n,
639            bool isupper,
640            bool isnormprovided,
641            double anorm,
642            ref double rc)
643        {
644            int i = 0;
645            int j = 0;
646            int kase = 0;
647            double ainvnm = 0;
648            double[] ex = new double[0];
649            double[] ev = new double[0];
650            double[] tmp = new double[0];
651            int[] iwork = new int[0];
652            double sa = 0;
653            double v = 0;
654            double maxgrowth = 0;
655            int i_ = 0;
656            int i1_ = 0;
657
658            System.Diagnostics.Debug.Assert(n>=1);
659            tmp = new double[n];
660           
661            //
662            // RC=0 if something happens
663            //
664            rc = 0;
665           
666            //
667            // prepare parameters for triangular solver
668            //
669            maxgrowth = 1/rcondthreshold();
670            sa = 0;
671            if( isupper )
672            {
673                for(i=0; i<=n-1; i++)
674                {
675                    for(j=i; j<=n-1; j++)
676                    {
677                        sa = Math.Max(sa, AP.Math.AbsComplex(cha[i,j]));
678                    }
679                }
680            }
681            else
682            {
683                for(i=0; i<=n-1; i++)
684                {
685                    for(j=0; j<=i; j++)
686                    {
687                        sa = Math.Max(sa, AP.Math.AbsComplex(cha[i,j]));
688                    }
689                }
690            }
691            if( (double)(sa)==(double)(0) )
692            {
693                sa = 1;
694            }
695            sa = 1/sa;
696           
697            //
698            // Estimate the norm of A.
699            //
700            if( !isnormprovided )
701            {
702                kase = 0;
703                anorm = 0;
704                while( true )
705                {
706                    rmatrixestimatenorm(n, ref ev, ref ex, ref iwork, ref anorm, ref kase);
707                    if( kase==0 )
708                    {
709                        break;
710                    }
711                    if( isupper )
712                    {
713                       
714                        //
715                        // Multiply by U
716                        //
717                        for(i=1; i<=n; i++)
718                        {
719                            i1_ = (i)-(i-1);
720                            v = 0.0;
721                            for(i_=i-1; i_<=n-1;i_++)
722                            {
723                                v += cha[i-1,i_]*ex[i_+i1_];
724                            }
725                            ex[i] = v;
726                        }
727                        for(i_=1; i_<=n;i_++)
728                        {
729                            ex[i_] = sa*ex[i_];
730                        }
731                       
732                        //
733                        // Multiply by U'
734                        //
735                        for(i=0; i<=n-1; i++)
736                        {
737                            tmp[i] = 0;
738                        }
739                        for(i=0; i<=n-1; i++)
740                        {
741                            v = ex[i+1];
742                            for(i_=i; i_<=n-1;i_++)
743                            {
744                                tmp[i_] = tmp[i_] + v*cha[i,i_];
745                            }
746                        }
747                        i1_ = (0) - (1);
748                        for(i_=1; i_<=n;i_++)
749                        {
750                            ex[i_] = tmp[i_+i1_];
751                        }
752                        for(i_=1; i_<=n;i_++)
753                        {
754                            ex[i_] = sa*ex[i_];
755                        }
756                    }
757                    else
758                    {
759                       
760                        //
761                        // Multiply by L'
762                        //
763                        for(i=0; i<=n-1; i++)
764                        {
765                            tmp[i] = 0;
766                        }
767                        for(i=0; i<=n-1; i++)
768                        {
769                            v = ex[i+1];
770                            for(i_=0; i_<=i;i_++)
771                            {
772                                tmp[i_] = tmp[i_] + v*cha[i,i_];
773                            }
774                        }
775                        i1_ = (0) - (1);
776                        for(i_=1; i_<=n;i_++)
777                        {
778                            ex[i_] = tmp[i_+i1_];
779                        }
780                        for(i_=1; i_<=n;i_++)
781                        {
782                            ex[i_] = sa*ex[i_];
783                        }
784                       
785                        //
786                        // Multiply by L
787                        //
788                        for(i=n; i>=1; i--)
789                        {
790                            i1_ = (1)-(0);
791                            v = 0.0;
792                            for(i_=0; i_<=i-1;i_++)
793                            {
794                                v += cha[i-1,i_]*ex[i_+i1_];
795                            }
796                            ex[i] = v;
797                        }
798                        for(i_=1; i_<=n;i_++)
799                        {
800                            ex[i_] = sa*ex[i_];
801                        }
802                    }
803                }
804            }
805           
806            //
807            // Quick return if possible
808            //
809            if( (double)(anorm)==(double)(0) )
810            {
811                return;
812            }
813            if( n==1 )
814            {
815                rc = 1;
816                return;
817            }
818           
819            //
820            // Estimate the 1-norm of inv(A).
821            //
822            kase = 0;
823            while( true )
824            {
825                rmatrixestimatenorm(n, ref ev, ref ex, ref iwork, ref ainvnm, ref kase);
826                if( kase==0 )
827                {
828                    break;
829                }
830                for(i=0; i<=n-1; i++)
831                {
832                    ex[i] = ex[i+1];
833                }
834                if( isupper )
835                {
836                   
837                    //
838                    // Multiply by inv(U').
839                    //
840                    if( !safesolve.rmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 1, false, maxgrowth) )
841                    {
842                        return;
843                    }
844                   
845                    //
846                    // Multiply by inv(U).
847                    //
848                    if( !safesolve.rmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 0, false, maxgrowth) )
849                    {
850                        return;
851                    }
852                }
853                else
854                {
855                   
856                    //
857                    // Multiply by inv(L).
858                    //
859                    if( !safesolve.rmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 0, false, maxgrowth) )
860                    {
861                        return;
862                    }
863                   
864                    //
865                    // Multiply by inv(L').
866                    //
867                    if( !safesolve.rmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 1, false, maxgrowth) )
868                    {
869                        return;
870                    }
871                }
872                for(i=n-1; i>=0; i--)
873                {
874                    ex[i+1] = ex[i];
875                }
876            }
877           
878            //
879            // Compute the estimate of the reciprocal condition number.
880            //
881            if( (double)(ainvnm)!=(double)(0) )
882            {
883                v = 1/ainvnm;
884                rc = v/anorm;
885                if( (double)(rc)<(double)(rcondthreshold()) )
886                {
887                    rc = 0;
888                }
889            }
890        }
891
892
893        /*************************************************************************
894        Internal subroutine for condition number estimation
895
896          -- LAPACK routine (version 3.0) --
897             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
898             Courant Institute, Argonne National Lab, and Rice University
899             February 29, 1992
900        *************************************************************************/
901        private static void hpdmatrixrcondcholeskyinternal(ref AP.Complex[,] cha,
902            int n,
903            bool isupper,
904            bool isnormprovided,
905            double anorm,
906            ref double rc)
907        {
908            int[] isave = new int[0];
909            double[] rsave = new double[0];
910            AP.Complex[] ex = new AP.Complex[0];
911            AP.Complex[] ev = new AP.Complex[0];
912            AP.Complex[] tmp = new AP.Complex[0];
913            int kase = 0;
914            double ainvnm = 0;
915            AP.Complex v = 0;
916            int i = 0;
917            int j = 0;
918            double sa = 0;
919            double maxgrowth = 0;
920            int i_ = 0;
921            int i1_ = 0;
922
923            System.Diagnostics.Debug.Assert(n>=1);
924            tmp = new AP.Complex[n];
925           
926            //
927            // RC=0 if something happens
928            //
929            rc = 0;
930           
931            //
932            // prepare parameters for triangular solver
933            //
934            maxgrowth = 1/rcondthreshold();
935            sa = 0;
936            if( isupper )
937            {
938                for(i=0; i<=n-1; i++)
939                {
940                    for(j=i; j<=n-1; j++)
941                    {
942                        sa = Math.Max(sa, AP.Math.AbsComplex(cha[i,j]));
943                    }
944                }
945            }
946            else
947            {
948                for(i=0; i<=n-1; i++)
949                {
950                    for(j=0; j<=i; j++)
951                    {
952                        sa = Math.Max(sa, AP.Math.AbsComplex(cha[i,j]));
953                    }
954                }
955            }
956            if( (double)(sa)==(double)(0) )
957            {
958                sa = 1;
959            }
960            sa = 1/sa;
961           
962            //
963            // Estimate the norm of A
964            //
965            if( !isnormprovided )
966            {
967                anorm = 0;
968                kase = 0;
969                while( true )
970                {
971                    cmatrixestimatenorm(n, ref ev, ref ex, ref anorm, ref kase, ref isave, ref rsave);
972                    if( kase==0 )
973                    {
974                        break;
975                    }
976                    if( isupper )
977                    {
978                       
979                        //
980                        // Multiply by U
981                        //
982                        for(i=1; i<=n; i++)
983                        {
984                            i1_ = (i)-(i-1);
985                            v = 0.0;
986                            for(i_=i-1; i_<=n-1;i_++)
987                            {
988                                v += cha[i-1,i_]*ex[i_+i1_];
989                            }
990                            ex[i] = v;
991                        }
992                        for(i_=1; i_<=n;i_++)
993                        {
994                            ex[i_] = sa*ex[i_];
995                        }
996                       
997                        //
998                        // Multiply by U'
999                        //
1000                        for(i=0; i<=n-1; i++)
1001                        {
1002                            tmp[i] = 0;
1003                        }
1004                        for(i=0; i<=n-1; i++)
1005                        {
1006                            v = ex[i+1];
1007                            for(i_=i; i_<=n-1;i_++)
1008                            {
1009                                tmp[i_] = tmp[i_] + v*AP.Math.Conj(cha[i,i_]);
1010                            }
1011                        }
1012                        i1_ = (0) - (1);
1013                        for(i_=1; i_<=n;i_++)
1014                        {
1015                            ex[i_] = tmp[i_+i1_];
1016                        }
1017                        for(i_=1; i_<=n;i_++)
1018                        {
1019                            ex[i_] = sa*ex[i_];
1020                        }
1021                    }
1022                    else
1023                    {
1024                       
1025                        //
1026                        // Multiply by L'
1027                        //
1028                        for(i=0; i<=n-1; i++)
1029                        {
1030                            tmp[i] = 0;
1031                        }
1032                        for(i=0; i<=n-1; i++)
1033                        {
1034                            v = ex[i+1];
1035                            for(i_=0; i_<=i;i_++)
1036                            {
1037                                tmp[i_] = tmp[i_] + v*AP.Math.Conj(cha[i,i_]);
1038                            }
1039                        }
1040                        i1_ = (0) - (1);
1041                        for(i_=1; i_<=n;i_++)
1042                        {
1043                            ex[i_] = tmp[i_+i1_];
1044                        }
1045                        for(i_=1; i_<=n;i_++)
1046                        {
1047                            ex[i_] = sa*ex[i_];
1048                        }
1049                       
1050                        //
1051                        // Multiply by L
1052                        //
1053                        for(i=n; i>=1; i--)
1054                        {
1055                            i1_ = (1)-(0);
1056                            v = 0.0;
1057                            for(i_=0; i_<=i-1;i_++)
1058                            {
1059                                v += cha[i-1,i_]*ex[i_+i1_];
1060                            }
1061                            ex[i] = v;
1062                        }
1063                        for(i_=1; i_<=n;i_++)
1064                        {
1065                            ex[i_] = sa*ex[i_];
1066                        }
1067                    }
1068                }
1069            }
1070           
1071            //
1072            // Quick return if possible
1073            // After this block we assume that ANORM<>0
1074            //
1075            if( (double)(anorm)==(double)(0) )
1076            {
1077                return;
1078            }
1079            if( n==1 )
1080            {
1081                rc = 1;
1082                return;
1083            }
1084           
1085            //
1086            // Estimate the norm of inv(A).
1087            //
1088            ainvnm = 0;
1089            kase = 0;
1090            while( true )
1091            {
1092                cmatrixestimatenorm(n, ref ev, ref ex, ref ainvnm, ref kase, ref isave, ref rsave);
1093                if( kase==0 )
1094                {
1095                    break;
1096                }
1097                for(i=0; i<=n-1; i++)
1098                {
1099                    ex[i] = ex[i+1];
1100                }
1101                if( isupper )
1102                {
1103                   
1104                    //
1105                    // Multiply by inv(U').
1106                    //
1107                    if( !safesolve.cmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 2, false, maxgrowth) )
1108                    {
1109                        return;
1110                    }
1111                   
1112                    //
1113                    // Multiply by inv(U).
1114                    //
1115                    if( !safesolve.cmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 0, false, maxgrowth) )
1116                    {
1117                        return;
1118                    }
1119                }
1120                else
1121                {
1122                   
1123                    //
1124                    // Multiply by inv(L).
1125                    //
1126                    if( !safesolve.cmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 0, false, maxgrowth) )
1127                    {
1128                        return;
1129                    }
1130                   
1131                    //
1132                    // Multiply by inv(L').
1133                    //
1134                    if( !safesolve.cmatrixscaledtrsafesolve(ref cha, sa, n, ref ex, isupper, 2, false, maxgrowth) )
1135                    {
1136                        return;
1137                    }
1138                }
1139                for(i=n-1; i>=0; i--)
1140                {
1141                    ex[i+1] = ex[i];
1142                }
1143            }
1144           
1145            //
1146            // Compute the estimate of the reciprocal condition number.
1147            //
1148            if( (double)(ainvnm)!=(double)(0) )
1149            {
1150                rc = 1/ainvnm;
1151                rc = rc/anorm;
1152                if( (double)(rc)<(double)(rcondthreshold()) )
1153                {
1154                    rc = 0;
1155                }
1156            }
1157        }
1158
1159
1160        /*************************************************************************
1161        Internal subroutine for condition number estimation
1162
1163          -- LAPACK routine (version 3.0) --
1164             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1165             Courant Institute, Argonne National Lab, and Rice University
1166             February 29, 1992
1167        *************************************************************************/
1168        private static void rmatrixrcondluinternal(ref double[,] lua,
1169            int n,
1170            bool onenorm,
1171            bool isanormprovided,
1172            double anorm,
1173            ref double rc)
1174        {
1175            double[] ex = new double[0];
1176            double[] ev = new double[0];
1177            int[] iwork = new int[0];
1178            double[] tmp = new double[0];
1179            double v = 0;
1180            int i = 0;
1181            int j = 0;
1182            int kase = 0;
1183            int kase1 = 0;
1184            double ainvnm = 0;
1185            double maxgrowth = 0;
1186            double su = 0;
1187            double sl = 0;
1188            bool mupper = new bool();
1189            bool mtrans = new bool();
1190            bool munit = new bool();
1191            int i_ = 0;
1192            int i1_ = 0;
1193
1194           
1195            //
1196            // RC=0 if something happens
1197            //
1198            rc = 0;
1199           
1200            //
1201            // init
1202            //
1203            if( onenorm )
1204            {
1205                kase1 = 1;
1206            }
1207            else
1208            {
1209                kase1 = 2;
1210            }
1211            mupper = true;
1212            mtrans = true;
1213            munit = true;
1214            iwork = new int[n+1];
1215            tmp = new double[n];
1216           
1217            //
1218            // prepare parameters for triangular solver
1219            //
1220            maxgrowth = 1/rcondthreshold();
1221            su = 0;
1222            sl = 1;
1223            for(i=0; i<=n-1; i++)
1224            {
1225                for(j=0; j<=i-1; j++)
1226                {
1227                    sl = Math.Max(sl, Math.Abs(lua[i,j]));
1228                }
1229                for(j=i; j<=n-1; j++)
1230                {
1231                    su = Math.Max(su, Math.Abs(lua[i,j]));
1232                }
1233            }
1234            if( (double)(su)==(double)(0) )
1235            {
1236                su = 1;
1237            }
1238            su = 1/su;
1239            sl = 1/sl;
1240           
1241            //
1242            // Estimate the norm of A.
1243            //
1244            if( !isanormprovided )
1245            {
1246                kase = 0;
1247                anorm = 0;
1248                while( true )
1249                {
1250                    rmatrixestimatenorm(n, ref ev, ref ex, ref iwork, ref anorm, ref kase);
1251                    if( kase==0 )
1252                    {
1253                        break;
1254                    }
1255                    if( kase==kase1 )
1256                    {
1257                       
1258                        //
1259                        // Multiply by U
1260                        //
1261                        for(i=1; i<=n; i++)
1262                        {
1263                            i1_ = (i)-(i-1);
1264                            v = 0.0;
1265                            for(i_=i-1; i_<=n-1;i_++)
1266                            {
1267                                v += lua[i-1,i_]*ex[i_+i1_];
1268                            }
1269                            ex[i] = v;
1270                        }
1271                       
1272                        //
1273                        // Multiply by L
1274                        //
1275                        for(i=n; i>=1; i--)
1276                        {
1277                            if( i>1 )
1278                            {
1279                                i1_ = (1)-(0);
1280                                v = 0.0;
1281                                for(i_=0; i_<=i-2;i_++)
1282                                {
1283                                    v += lua[i-1,i_]*ex[i_+i1_];
1284                                }
1285                            }
1286                            else
1287                            {
1288                                v = 0;
1289                            }
1290                            ex[i] = ex[i]+v;
1291                        }
1292                    }
1293                    else
1294                    {
1295                       
1296                        //
1297                        // Multiply by L'
1298                        //
1299                        for(i=0; i<=n-1; i++)
1300                        {
1301                            tmp[i] = 0;
1302                        }
1303                        for(i=0; i<=n-1; i++)
1304                        {
1305                            v = ex[i+1];
1306                            if( i>=1 )
1307                            {
1308                                for(i_=0; i_<=i-1;i_++)
1309                                {
1310                                    tmp[i_] = tmp[i_] + v*lua[i,i_];
1311                                }
1312                            }
1313                            tmp[i] = tmp[i]+v;
1314                        }
1315                        i1_ = (0) - (1);
1316                        for(i_=1; i_<=n;i_++)
1317                        {
1318                            ex[i_] = tmp[i_+i1_];
1319                        }
1320                       
1321                        //
1322                        // Multiply by U'
1323                        //
1324                        for(i=0; i<=n-1; i++)
1325                        {
1326                            tmp[i] = 0;
1327                        }
1328                        for(i=0; i<=n-1; i++)
1329                        {
1330                            v = ex[i+1];
1331                            for(i_=i; i_<=n-1;i_++)
1332                            {
1333                                tmp[i_] = tmp[i_] + v*lua[i,i_];
1334                            }
1335                        }
1336                        i1_ = (0) - (1);
1337                        for(i_=1; i_<=n;i_++)
1338                        {
1339                            ex[i_] = tmp[i_+i1_];
1340                        }
1341                    }
1342                }
1343            }
1344           
1345            //
1346            // Scale according to SU/SL
1347            //
1348            anorm = anorm*su*sl;
1349           
1350            //
1351            // Quick return if possible
1352            // We assume that ANORM<>0 after this block
1353            //
1354            if( (double)(anorm)==(double)(0) )
1355            {
1356                return;
1357            }
1358            if( n==1 )
1359            {
1360                rc = 1;
1361                return;
1362            }
1363           
1364            //
1365            // Estimate the norm of inv(A).
1366            //
1367            ainvnm = 0;
1368            kase = 0;
1369            while( true )
1370            {
1371                rmatrixestimatenorm(n, ref ev, ref ex, ref iwork, ref ainvnm, ref kase);
1372                if( kase==0 )
1373                {
1374                    break;
1375                }
1376               
1377                //
1378                // from 1-based array to 0-based
1379                //
1380                for(i=0; i<=n-1; i++)
1381                {
1382                    ex[i] = ex[i+1];
1383                }
1384               
1385                //
1386                // multiply by inv(A) or inv(A')
1387                //
1388                if( kase==kase1 )
1389                {
1390                   
1391                    //
1392                    // Multiply by inv(L).
1393                    //
1394                    if( !safesolve.rmatrixscaledtrsafesolve(ref lua, sl, n, ref ex, !mupper, 0, munit, maxgrowth) )
1395                    {
1396                        return;
1397                    }
1398                   
1399                    //
1400                    // Multiply by inv(U).
1401                    //
1402                    if( !safesolve.rmatrixscaledtrsafesolve(ref lua, su, n, ref ex, mupper, 0, !munit, maxgrowth) )
1403                    {
1404                        return;
1405                    }
1406                }
1407                else
1408                {
1409                   
1410                    //
1411                    // Multiply by inv(U').
1412                    //
1413                    if( !safesolve.rmatrixscaledtrsafesolve(ref lua, su, n, ref ex, mupper, 1, !munit, maxgrowth) )
1414                    {
1415                        return;
1416                    }
1417                   
1418                    //
1419                    // Multiply by inv(L').
1420                    //
1421                    if( !safesolve.rmatrixscaledtrsafesolve(ref lua, sl, n, ref ex, !mupper, 1, munit, maxgrowth) )
1422                    {
1423                        return;
1424                    }
1425                }
1426               
1427                //
1428                // from 0-based array to 1-based
1429                //
1430                for(i=n-1; i>=0; i--)
1431                {
1432                    ex[i+1] = ex[i];
1433                }
1434            }
1435           
1436            //
1437            // Compute the estimate of the reciprocal condition number.
1438            //
1439            if( (double)(ainvnm)!=(double)(0) )
1440            {
1441                rc = 1/ainvnm;
1442                rc = rc/anorm;
1443                if( (double)(rc)<(double)(rcondthreshold()) )
1444                {
1445                    rc = 0;
1446                }
1447            }
1448        }
1449
1450
1451        /*************************************************************************
1452        Condition number estimation
1453
1454          -- LAPACK routine (version 3.0) --
1455             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1456             Courant Institute, Argonne National Lab, and Rice University
1457             March 31, 1993
1458        *************************************************************************/
1459        private static void cmatrixrcondluinternal(ref AP.Complex[,] lua,
1460            int n,
1461            bool onenorm,
1462            bool isanormprovided,
1463            double anorm,
1464            ref double rc)
1465        {
1466            AP.Complex[] ex = new AP.Complex[0];
1467            AP.Complex[] cwork2 = new AP.Complex[0];
1468            AP.Complex[] cwork3 = new AP.Complex[0];
1469            AP.Complex[] cwork4 = new AP.Complex[0];
1470            int[] isave = new int[0];
1471            double[] rsave = new double[0];
1472            int kase = 0;
1473            int kase1 = 0;
1474            double ainvnm = 0;
1475            AP.Complex v = 0;
1476            int i = 0;
1477            int j = 0;
1478            double su = 0;
1479            double sl = 0;
1480            double maxgrowth = 0;
1481            int i_ = 0;
1482            int i1_ = 0;
1483
1484            if( n<=0 )
1485            {
1486                return;
1487            }
1488            cwork2 = new AP.Complex[n+1];
1489            rc = 0;
1490            if( n==0 )
1491            {
1492                rc = 1;
1493                return;
1494            }
1495           
1496            //
1497            // prepare parameters for triangular solver
1498            //
1499            maxgrowth = 1/rcondthreshold();
1500            su = 0;
1501            sl = 1;
1502            for(i=0; i<=n-1; i++)
1503            {
1504                for(j=0; j<=i-1; j++)
1505                {
1506                    sl = Math.Max(sl, AP.Math.AbsComplex(lua[i,j]));
1507                }
1508                for(j=i; j<=n-1; j++)
1509                {
1510                    su = Math.Max(su, AP.Math.AbsComplex(lua[i,j]));
1511                }
1512            }
1513            if( (double)(su)==(double)(0) )
1514            {
1515                su = 1;
1516            }
1517            su = 1/su;
1518            sl = 1/sl;
1519           
1520            //
1521            // Estimate the norm of SU*SL*A.
1522            //
1523            if( !isanormprovided )
1524            {
1525                anorm = 0;
1526                if( onenorm )
1527                {
1528                    kase1 = 1;
1529                }
1530                else
1531                {
1532                    kase1 = 2;
1533                }
1534                kase = 0;
1535                do
1536                {
1537                    cmatrixestimatenorm(n, ref cwork4, ref ex, ref anorm, ref kase, ref isave, ref rsave);
1538                    if( kase!=0 )
1539                    {
1540                        if( kase==kase1 )
1541                        {
1542                           
1543                            //
1544                            // Multiply by U
1545                            //
1546                            for(i=1; i<=n; i++)
1547                            {
1548                                i1_ = (i)-(i-1);
1549                                v = 0.0;
1550                                for(i_=i-1; i_<=n-1;i_++)
1551                                {
1552                                    v += lua[i-1,i_]*ex[i_+i1_];
1553                                }
1554                                ex[i] = v;
1555                            }
1556                           
1557                            //
1558                            // Multiply by L
1559                            //
1560                            for(i=n; i>=1; i--)
1561                            {
1562                                v = 0;
1563                                if( i>1 )
1564                                {
1565                                    i1_ = (1)-(0);
1566                                    v = 0.0;
1567                                    for(i_=0; i_<=i-2;i_++)
1568                                    {
1569                                        v += lua[i-1,i_]*ex[i_+i1_];
1570                                    }
1571                                }
1572                                ex[i] = v+ex[i];
1573                            }
1574                        }
1575                        else
1576                        {
1577                           
1578                            //
1579                            // Multiply by L'
1580                            //
1581                            for(i=1; i<=n; i++)
1582                            {
1583                                cwork2[i] = 0;
1584                            }
1585                            for(i=1; i<=n; i++)
1586                            {
1587                                v = ex[i];
1588                                if( i>1 )
1589                                {
1590                                    i1_ = (0) - (1);
1591                                    for(i_=1; i_<=i-1;i_++)
1592                                    {
1593                                        cwork2[i_] = cwork2[i_] + v*AP.Math.Conj(lua[i-1,i_+i1_]);
1594                                    }
1595                                }
1596                                cwork2[i] = cwork2[i]+v;
1597                            }
1598                           
1599                            //
1600                            // Multiply by U'
1601                            //
1602                            for(i=1; i<=n; i++)
1603                            {
1604                                ex[i] = 0;
1605                            }
1606                            for(i=1; i<=n; i++)
1607                            {
1608                                v = cwork2[i];
1609                                i1_ = (i-1) - (i);
1610                                for(i_=i; i_<=n;i_++)
1611                                {
1612                                    ex[i_] = ex[i_] + v*AP.Math.Conj(lua[i-1,i_+i1_]);
1613                                }
1614                            }
1615                        }
1616                    }
1617                }
1618                while( kase!=0 );
1619            }
1620           
1621            //
1622            // Scale according to SU/SL
1623            //
1624            anorm = anorm*su*sl;
1625           
1626            //
1627            // Quick return if possible
1628            //
1629            if( (double)(anorm)==(double)(0) )
1630            {
1631                return;
1632            }
1633           
1634            //
1635            // Estimate the norm of inv(A).
1636            //
1637            ainvnm = 0;
1638            if( onenorm )
1639            {
1640                kase1 = 1;
1641            }
1642            else
1643            {
1644                kase1 = 2;
1645            }
1646            kase = 0;
1647            while( true )
1648            {
1649                cmatrixestimatenorm(n, ref cwork4, ref ex, ref ainvnm, ref kase, ref isave, ref rsave);
1650                if( kase==0 )
1651                {
1652                    break;
1653                }
1654               
1655                //
1656                // From 1-based to 0-based
1657                //
1658                for(i=0; i<=n-1; i++)
1659                {
1660                    ex[i] = ex[i+1];
1661                }
1662               
1663                //
1664                // multiply by inv(A) or inv(A')
1665                //
1666                if( kase==kase1 )
1667                {
1668                   
1669                    //
1670                    // Multiply by inv(L).
1671                    //
1672                    if( !safesolve.cmatrixscaledtrsafesolve(ref lua, sl, n, ref ex, false, 0, true, maxgrowth) )
1673                    {
1674                        rc = 0;
1675                        return;
1676                    }
1677                   
1678                    //
1679                    // Multiply by inv(U).
1680                    //
1681                    if( !safesolve.cmatrixscaledtrsafesolve(ref lua, su, n, ref ex, true, 0, false, maxgrowth) )
1682                    {
1683                        rc = 0;
1684                        return;
1685                    }
1686                }
1687                else
1688                {
1689                   
1690                    //
1691                    // Multiply by inv(U').
1692                    //
1693                    if( !safesolve.cmatrixscaledtrsafesolve(ref lua, su, n, ref ex, true, 2, false, maxgrowth) )
1694                    {
1695                        rc = 0;
1696                        return;
1697                    }
1698                   
1699                    //
1700                    // Multiply by inv(L').
1701                    //
1702                    if( !safesolve.cmatrixscaledtrsafesolve(ref lua, sl, n, ref ex, false, 2, true, maxgrowth) )
1703                    {
1704                        rc = 0;
1705                        return;
1706                    }
1707                }
1708               
1709                //
1710                // from 0-based to 1-based
1711                //
1712                for(i=n-1; i>=0; i--)
1713                {
1714                    ex[i+1] = ex[i];
1715                }
1716            }
1717           
1718            //
1719            // Compute the estimate of the reciprocal condition number.
1720            //
1721            if( (double)(ainvnm)!=(double)(0) )
1722            {
1723                rc = 1/ainvnm;
1724                rc = rc/anorm;
1725                if( (double)(rc)<(double)(rcondthreshold()) )
1726                {
1727                    rc = 0;
1728                }
1729            }
1730        }
1731
1732
1733        /*************************************************************************
1734        Internal subroutine for matrix norm estimation
1735
1736          -- LAPACK auxiliary routine (version 3.0) --
1737             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1738             Courant Institute, Argonne National Lab, and Rice University
1739             February 29, 1992
1740        *************************************************************************/
1741        private static void rmatrixestimatenorm(int n,
1742            ref double[] v,
1743            ref double[] x,
1744            ref int[] isgn,
1745            ref double est,
1746            ref int kase)
1747        {
1748            int itmax = 0;
1749            int i = 0;
1750            double t = 0;
1751            bool flg = new bool();
1752            int positer = 0;
1753            int posj = 0;
1754            int posjlast = 0;
1755            int posjump = 0;
1756            int posaltsgn = 0;
1757            int posestold = 0;
1758            int postemp = 0;
1759            int i_ = 0;
1760
1761            itmax = 5;
1762            posaltsgn = n+1;
1763            posestold = n+2;
1764            postemp = n+3;
1765            positer = n+1;
1766            posj = n+2;
1767            posjlast = n+3;
1768            posjump = n+4;
1769            if( kase==0 )
1770            {
1771                v = new double[n+4];
1772                x = new double[n+1];
1773                isgn = new int[n+5];
1774                t = (double)(1)/(double)(n);
1775                for(i=1; i<=n; i++)
1776                {
1777                    x[i] = t;
1778                }
1779                kase = 1;
1780                isgn[posjump] = 1;
1781                return;
1782            }
1783           
1784            //
1785            //     ................ ENTRY   (JUMP = 1)
1786            //     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
1787            //
1788            if( isgn[posjump]==1 )
1789            {
1790                if( n==1 )
1791                {
1792                    v[1] = x[1];
1793                    est = Math.Abs(v[1]);
1794                    kase = 0;
1795                    return;
1796                }
1797                est = 0;
1798                for(i=1; i<=n; i++)
1799                {
1800                    est = est+Math.Abs(x[i]);
1801                }
1802                for(i=1; i<=n; i++)
1803                {
1804                    if( (double)(x[i])>=(double)(0) )
1805                    {
1806                        x[i] = 1;
1807                    }
1808                    else
1809                    {
1810                        x[i] = -1;
1811                    }
1812                    isgn[i] = Math.Sign(x[i]);
1813                }
1814                kase = 2;
1815                isgn[posjump] = 2;
1816                return;
1817            }
1818           
1819            //
1820            //     ................ ENTRY   (JUMP = 2)
1821            //     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
1822            //
1823            if( isgn[posjump]==2 )
1824            {
1825                isgn[posj] = 1;
1826                for(i=2; i<=n; i++)
1827                {
1828                    if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
1829                    {
1830                        isgn[posj] = i;
1831                    }
1832                }
1833                isgn[positer] = 2;
1834               
1835                //
1836                // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
1837                //
1838                for(i=1; i<=n; i++)
1839                {
1840                    x[i] = 0;
1841                }
1842                x[isgn[posj]] = 1;
1843                kase = 1;
1844                isgn[posjump] = 3;
1845                return;
1846            }
1847           
1848            //
1849            //     ................ ENTRY   (JUMP = 3)
1850            //     X HAS BEEN OVERWRITTEN BY A*X.
1851            //
1852            if( isgn[posjump]==3 )
1853            {
1854                for(i_=1; i_<=n;i_++)
1855                {
1856                    v[i_] = x[i_];
1857                }
1858                v[posestold] = est;
1859                est = 0;
1860                for(i=1; i<=n; i++)
1861                {
1862                    est = est+Math.Abs(v[i]);
1863                }
1864                flg = false;
1865                for(i=1; i<=n; i++)
1866                {
1867                    if( (double)(x[i])>=(double)(0) & isgn[i]<0 | (double)(x[i])<(double)(0) & isgn[i]>=0 )
1868                    {
1869                        flg = true;
1870                    }
1871                }
1872               
1873                //
1874                // REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
1875                // OR MAY BE CYCLING.
1876                //
1877                if( !flg | (double)(est)<=(double)(v[posestold]) )
1878                {
1879                    v[posaltsgn] = 1;
1880                    for(i=1; i<=n; i++)
1881                    {
1882                        x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
1883                        v[posaltsgn] = -v[posaltsgn];
1884                    }
1885                    kase = 1;
1886                    isgn[posjump] = 5;
1887                    return;
1888                }
1889                for(i=1; i<=n; i++)
1890                {
1891                    if( (double)(x[i])>=(double)(0) )
1892                    {
1893                        x[i] = 1;
1894                        isgn[i] = 1;
1895                    }
1896                    else
1897                    {
1898                        x[i] = -1;
1899                        isgn[i] = -1;
1900                    }
1901                }
1902                kase = 2;
1903                isgn[posjump] = 4;
1904                return;
1905            }
1906           
1907            //
1908            //     ................ ENTRY   (JUMP = 4)
1909            //     X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
1910            //
1911            if( isgn[posjump]==4 )
1912            {
1913                isgn[posjlast] = isgn[posj];
1914                isgn[posj] = 1;
1915                for(i=2; i<=n; i++)
1916                {
1917                    if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
1918                    {
1919                        isgn[posj] = i;
1920                    }
1921                }
1922                if( (double)(x[isgn[posjlast]])!=(double)(Math.Abs(x[isgn[posj]])) & isgn[positer]<itmax )
1923                {
1924                    isgn[positer] = isgn[positer]+1;
1925                    for(i=1; i<=n; i++)
1926                    {
1927                        x[i] = 0;
1928                    }
1929                    x[isgn[posj]] = 1;
1930                    kase = 1;
1931                    isgn[posjump] = 3;
1932                    return;
1933                }
1934               
1935                //
1936                // ITERATION COMPLETE.  FINAL STAGE.
1937                //
1938                v[posaltsgn] = 1;
1939                for(i=1; i<=n; i++)
1940                {
1941                    x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
1942                    v[posaltsgn] = -v[posaltsgn];
1943                }
1944                kase = 1;
1945                isgn[posjump] = 5;
1946                return;
1947            }
1948           
1949            //
1950            //     ................ ENTRY   (JUMP = 5)
1951            //     X HAS BEEN OVERWRITTEN BY A*X.
1952            //
1953            if( isgn[posjump]==5 )
1954            {
1955                v[postemp] = 0;
1956                for(i=1; i<=n; i++)
1957                {
1958                    v[postemp] = v[postemp]+Math.Abs(x[i]);
1959                }
1960                v[postemp] = 2*v[postemp]/(3*n);
1961                if( (double)(v[postemp])>(double)(est) )
1962                {
1963                    for(i_=1; i_<=n;i_++)
1964                    {
1965                        v[i_] = x[i_];
1966                    }
1967                    est = v[postemp];
1968                }
1969                kase = 0;
1970                return;
1971            }
1972        }
1973
1974
1975        private static void cmatrixestimatenorm(int n,
1976            ref AP.Complex[] v,
1977            ref AP.Complex[] x,
1978            ref double est,
1979            ref int kase,
1980            ref int[] isave,
1981            ref double[] rsave)
1982        {
1983            int itmax = 0;
1984            int i = 0;
1985            int iter = 0;
1986            int j = 0;
1987            int jlast = 0;
1988            int jump = 0;
1989            double absxi = 0;
1990            double altsgn = 0;
1991            double estold = 0;
1992            double safmin = 0;
1993            double temp = 0;
1994            int i_ = 0;
1995
1996           
1997            //
1998            //Executable Statements ..
1999            //
2000            itmax = 5;
2001            safmin = AP.Math.MinRealNumber;
2002            if( kase==0 )
2003            {
2004                v = new AP.Complex[n+1];
2005                x = new AP.Complex[n+1];
2006                isave = new int[5];
2007                rsave = new double[4];
2008                for(i=1; i<=n; i++)
2009                {
2010                    x[i] = (double)(1)/(double)(n);
2011                }
2012                kase = 1;
2013                jump = 1;
2014                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2015                return;
2016            }
2017            internalcomplexrcondloadall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2018           
2019            //
2020            // ENTRY   (JUMP = 1)
2021            // FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
2022            //
2023            if( jump==1 )
2024            {
2025                if( n==1 )
2026                {
2027                    v[1] = x[1];
2028                    est = AP.Math.AbsComplex(v[1]);
2029                    kase = 0;
2030                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2031                    return;
2032                }
2033                est = internalcomplexrcondscsum1(ref x, n);
2034                for(i=1; i<=n; i++)
2035                {
2036                    absxi = AP.Math.AbsComplex(x[i]);
2037                    if( (double)(absxi)>(double)(safmin) )
2038                    {
2039                        x[i] = x[i]/absxi;
2040                    }
2041                    else
2042                    {
2043                        x[i] = 1;
2044                    }
2045                }
2046                kase = 2;
2047                jump = 2;
2048                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2049                return;
2050            }
2051           
2052            //
2053            // ENTRY   (JUMP = 2)
2054            // FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
2055            //
2056            if( jump==2 )
2057            {
2058                j = internalcomplexrcondicmax1(ref x, n);
2059                iter = 2;
2060               
2061                //
2062                // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
2063                //
2064                for(i=1; i<=n; i++)
2065                {
2066                    x[i] = 0;
2067                }
2068                x[j] = 1;
2069                kase = 1;
2070                jump = 3;
2071                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2072                return;
2073            }
2074           
2075            //
2076            // ENTRY   (JUMP = 3)
2077            // X HAS BEEN OVERWRITTEN BY A*X.
2078            //
2079            if( jump==3 )
2080            {
2081                for(i_=1; i_<=n;i_++)
2082                {
2083                    v[i_] = x[i_];
2084                }
2085                estold = est;
2086                est = internalcomplexrcondscsum1(ref v, n);
2087               
2088                //
2089                // TEST FOR CYCLING.
2090                //
2091                if( (double)(est)<=(double)(estold) )
2092                {
2093                   
2094                    //
2095                    // ITERATION COMPLETE.  FINAL STAGE.
2096                    //
2097                    altsgn = 1;
2098                    for(i=1; i<=n; i++)
2099                    {
2100                        x[i] = altsgn*(1+((double)(i-1))/((double)(n-1)));
2101                        altsgn = -altsgn;
2102                    }
2103                    kase = 1;
2104                    jump = 5;
2105                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2106                    return;
2107                }
2108                for(i=1; i<=n; i++)
2109                {
2110                    absxi = AP.Math.AbsComplex(x[i]);
2111                    if( (double)(absxi)>(double)(safmin) )
2112                    {
2113                        x[i] = x[i]/absxi;
2114                    }
2115                    else
2116                    {
2117                        x[i] = 1;
2118                    }
2119                }
2120                kase = 2;
2121                jump = 4;
2122                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2123                return;
2124            }
2125           
2126            //
2127            // ENTRY   (JUMP = 4)
2128            // X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
2129            //
2130            if( jump==4 )
2131            {
2132                jlast = j;
2133                j = internalcomplexrcondicmax1(ref x, n);
2134                if( (double)(AP.Math.AbsComplex(x[jlast]))!=(double)(AP.Math.AbsComplex(x[j])) & iter<itmax )
2135                {
2136                    iter = iter+1;
2137                   
2138                    //
2139                    // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
2140                    //
2141                    for(i=1; i<=n; i++)
2142                    {
2143                        x[i] = 0;
2144                    }
2145                    x[j] = 1;
2146                    kase = 1;
2147                    jump = 3;
2148                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2149                    return;
2150                }
2151               
2152                //
2153                // ITERATION COMPLETE.  FINAL STAGE.
2154                //
2155                altsgn = 1;
2156                for(i=1; i<=n; i++)
2157                {
2158                    x[i] = altsgn*(1+((double)(i-1))/((double)(n-1)));
2159                    altsgn = -altsgn;
2160                }
2161                kase = 1;
2162                jump = 5;
2163                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2164                return;
2165            }
2166           
2167            //
2168            // ENTRY   (JUMP = 5)
2169            // X HAS BEEN OVERWRITTEN BY A*X.
2170            //
2171            if( jump==5 )
2172            {
2173                temp = 2*(internalcomplexrcondscsum1(ref x, n)/(3*n));
2174                if( (double)(temp)>(double)(est) )
2175                {
2176                    for(i_=1; i_<=n;i_++)
2177                    {
2178                        v[i_] = x[i_];
2179                    }
2180                    est = temp;
2181                }
2182                kase = 0;
2183                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
2184                return;
2185            }
2186        }
2187
2188
2189        private static double internalcomplexrcondscsum1(ref AP.Complex[] x,
2190            int n)
2191        {
2192            double result = 0;
2193            int i = 0;
2194
2195            result = 0;
2196            for(i=1; i<=n; i++)
2197            {
2198                result = result+AP.Math.AbsComplex(x[i]);
2199            }
2200            return result;
2201        }
2202
2203
2204        private static int internalcomplexrcondicmax1(ref AP.Complex[] x,
2205            int n)
2206        {
2207            int result = 0;
2208            int i = 0;
2209            double m = 0;
2210
2211            result = 1;
2212            m = AP.Math.AbsComplex(x[1]);
2213            for(i=2; i<=n; i++)
2214            {
2215                if( (double)(AP.Math.AbsComplex(x[i]))>(double)(m) )
2216                {
2217                    result = i;
2218                    m = AP.Math.AbsComplex(x[i]);
2219                }
2220            }
2221            return result;
2222        }
2223
2224
2225        private static void internalcomplexrcondsaveall(ref int[] isave,
2226            ref double[] rsave,
2227            ref int i,
2228            ref int iter,
2229            ref int j,
2230            ref int jlast,
2231            ref int jump,
2232            ref double absxi,
2233            ref double altsgn,
2234            ref double estold,
2235            ref double temp)
2236        {
2237            isave[0] = i;
2238            isave[1] = iter;
2239            isave[2] = j;
2240            isave[3] = jlast;
2241            isave[4] = jump;
2242            rsave[0] = absxi;
2243            rsave[1] = altsgn;
2244            rsave[2] = estold;
2245            rsave[3] = temp;
2246        }
2247
2248
2249        private static void internalcomplexrcondloadall(ref int[] isave,
2250            ref double[] rsave,
2251            ref int i,
2252            ref int iter,
2253            ref int j,
2254            ref int jlast,
2255            ref int jump,
2256            ref double absxi,
2257            ref double altsgn,
2258            ref double estold,
2259            ref double temp)
2260        {
2261            i = isave[0];
2262            iter = isave[1];
2263            j = isave[2];
2264            jlast = isave[3];
2265            jump = isave[4];
2266            absxi = rsave[0];
2267            altsgn = rsave[1];
2268            estold = rsave[2];
2269            temp = rsave[3];
2270        }
2271    }
2272}
Note: See TracBrowser for help on using the repository browser.