Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/rcond.cs @ 3031

Last change on this file since 3031 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)
File size: 23.2 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        public static double rmatrixrcond1(ref double[,] a,
47            int n)
48        {
49            double result = 0;
50            int i = 0;
51            double[,] a1 = new double[0,0];
52            int i_ = 0;
53            int i1_ = 0;
54
55            System.Diagnostics.Debug.Assert(n>=1, "RMatrixRCond1: N<1!");
56            a1 = new double[n+1, n+1];
57            for(i=1; i<=n; i++)
58            {
59                i1_ = (0) - (1);
60                for(i_=1; i_<=n;i_++)
61                {
62                    a1[i,i_] = a[i-1,i_+i1_];
63                }
64            }
65            result = rcond1(a1, n);
66            return result;
67        }
68
69
70        /*************************************************************************
71        Estimate of the condition number of a matrix given by its LU decomposition (1-norm)
72
73        The algorithm calculates a lower bound of the condition number. In this case,
74        the algorithm does not return a lower bound of the condition number, but an
75        inverse number (to avoid an overflow in case of a singular matrix).
76
77        Input parameters:
78            LUDcmp      -   LU decomposition of a matrix in compact form. Output of
79                            the RMatrixLU subroutine.
80            N           -   size of matrix A.
81
82        Result: 1/LowerBound(cond(A))
83        *************************************************************************/
84        public static double rmatrixlurcond1(ref double[,] ludcmp,
85            int n)
86        {
87            double result = 0;
88            int i = 0;
89            double[,] a1 = new double[0,0];
90            int i_ = 0;
91            int i1_ = 0;
92
93            System.Diagnostics.Debug.Assert(n>=1, "RMatrixLURCond1: N<1!");
94            a1 = new double[n+1, n+1];
95            for(i=1; i<=n; i++)
96            {
97                i1_ = (0) - (1);
98                for(i_=1; i_<=n;i_++)
99                {
100                    a1[i,i_] = ludcmp[i-1,i_+i1_];
101                }
102            }
103            result = rcond1lu(ref a1, n);
104            return result;
105        }
106
107
108        /*************************************************************************
109        Estimate of a matrix condition number (infinity-norm).
110
111        The algorithm calculates a lower bound of the condition number. In this case,
112        the algorithm does not return a lower bound of the condition number, but an
113        inverse number (to avoid an overflow in case of a singular matrix).
114
115        Input parameters:
116            A   -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
117            N   -   size of matrix A.
118
119        Result: 1/LowerBound(cond(A))
120        *************************************************************************/
121        public static double rmatrixrcondinf(ref double[,] a,
122            int n)
123        {
124            double result = 0;
125            int i = 0;
126            double[,] a1 = new double[0,0];
127            int i_ = 0;
128            int i1_ = 0;
129
130            System.Diagnostics.Debug.Assert(n>=1, "RMatrixRCondInf: N<1!");
131            a1 = new double[n+1, n+1];
132            for(i=1; i<=n; i++)
133            {
134                i1_ = (0) - (1);
135                for(i_=1; i_<=n;i_++)
136                {
137                    a1[i,i_] = a[i-1,i_+i1_];
138                }
139            }
140            result = rcondinf(a1, n);
141            return result;
142        }
143
144
145        /*************************************************************************
146        Estimate of the condition number of a matrix given by its LU decomposition
147        (infinity norm).
148
149        The algorithm calculates a lower bound of the condition number. In this case,
150        the algorithm does not return a lower bound of the condition number, but an
151        inverse number (to avoid an overflow in case of a singular matrix).
152
153        Input parameters:
154            LUDcmp  -   LU decomposition of a matrix in compact form. Output of
155                        the RMatrixLU subroutine.
156            N       -   size of matrix A.
157
158        Result: 1/LowerBound(cond(A))
159        *************************************************************************/
160        public static double rmatrixlurcondinf(ref double[,] ludcmp,
161            int n)
162        {
163            double result = 0;
164            int i = 0;
165            double[,] a1 = new double[0,0];
166            int i_ = 0;
167            int i1_ = 0;
168
169            System.Diagnostics.Debug.Assert(n>=1, "RMatrixLURCondInf: N<1!");
170            a1 = new double[n+1, n+1];
171            for(i=1; i<=n; i++)
172            {
173                i1_ = (0) - (1);
174                for(i_=1; i_<=n;i_++)
175                {
176                    a1[i,i_] = ludcmp[i-1,i_+i1_];
177                }
178            }
179            result = rcondinflu(ref a1, n);
180            return result;
181        }
182
183
184        public static double rcond1(double[,] a,
185            int n)
186        {
187            double result = 0;
188            int i = 0;
189            int j = 0;
190            double v = 0;
191            double nrm = 0;
192            int[] pivots = new int[0];
193
194            a = (double[,])a.Clone();
195
196            nrm = 0;
197            for(j=1; j<=n; j++)
198            {
199                v = 0;
200                for(i=1; i<=n; i++)
201                {
202                    v = v+Math.Abs(a[i,j]);
203                }
204                nrm = Math.Max(nrm, v);
205            }
206            lu.ludecomposition(ref a, n, n, ref pivots);
207            internalestimatercondlu(ref a, n, true, true, nrm, ref v);
208            result = v;
209            return result;
210        }
211
212
213        public static double rcond1lu(ref double[,] ludcmp,
214            int n)
215        {
216            double result = 0;
217            double v = 0;
218
219            internalestimatercondlu(ref ludcmp, n, true, false, 0, ref v);
220            result = v;
221            return result;
222        }
223
224
225        public static double rcondinf(double[,] a,
226            int n)
227        {
228            double result = 0;
229            int i = 0;
230            int j = 0;
231            double v = 0;
232            double nrm = 0;
233            int[] pivots = new int[0];
234
235            a = (double[,])a.Clone();
236
237            nrm = 0;
238            for(i=1; i<=n; i++)
239            {
240                v = 0;
241                for(j=1; j<=n; j++)
242                {
243                    v = v+Math.Abs(a[i,j]);
244                }
245                nrm = Math.Max(nrm, v);
246            }
247            lu.ludecomposition(ref a, n, n, ref pivots);
248            internalestimatercondlu(ref a, n, false, true, nrm, ref v);
249            result = v;
250            return result;
251        }
252
253
254        public static double rcondinflu(ref double[,] ludcmp,
255            int n)
256        {
257            double result = 0;
258            double v = 0;
259
260            internalestimatercondlu(ref ludcmp, n, false, false, 0, ref v);
261            result = v;
262            return result;
263        }
264
265
266        private static void internalestimatercondlu(ref double[,] ludcmp,
267            int n,
268            bool onenorm,
269            bool isanormprovided,
270            double anorm,
271            ref double rc)
272        {
273            double[] work0 = new double[0];
274            double[] work1 = new double[0];
275            double[] work2 = new double[0];
276            double[] work3 = new double[0];
277            int[] iwork = new int[0];
278            double v = 0;
279            bool normin = new bool();
280            int i = 0;
281            int im1 = 0;
282            int ip1 = 0;
283            int ix = 0;
284            int kase = 0;
285            int kase1 = 0;
286            double ainvnm = 0;
287            double ascale = 0;
288            double sl = 0;
289            double smlnum = 0;
290            double su = 0;
291            bool mupper = new bool();
292            bool mtrans = new bool();
293            bool munit = new bool();
294            int i_ = 0;
295
296           
297            //
298            // Quick return if possible
299            //
300            if( n==0 )
301            {
302                rc = 1;
303                return;
304            }
305           
306            //
307            // init
308            //
309            if( onenorm )
310            {
311                kase1 = 1;
312            }
313            else
314            {
315                kase1 = 2;
316            }
317            mupper = true;
318            mtrans = true;
319            munit = true;
320            work0 = new double[n+1];
321            work1 = new double[n+1];
322            work2 = new double[n+1];
323            work3 = new double[n+1];
324            iwork = new int[n+1];
325           
326            //
327            // Estimate the norm of A.
328            //
329            if( !isanormprovided )
330            {
331                kase = 0;
332                anorm = 0;
333                while( true )
334                {
335                    internalestimatenorm(n, ref work1, ref work0, ref iwork, ref anorm, ref kase);
336                    if( kase==0 )
337                    {
338                        break;
339                    }
340                    if( kase==kase1 )
341                    {
342                       
343                        //
344                        // Multiply by U
345                        //
346                        for(i=1; i<=n; i++)
347                        {
348                            v = 0.0;
349                            for(i_=i; i_<=n;i_++)
350                            {
351                                v += ludcmp[i,i_]*work0[i_];
352                            }
353                            work0[i] = v;
354                        }
355                       
356                        //
357                        // Multiply by L
358                        //
359                        for(i=n; i>=1; i--)
360                        {
361                            im1 = i-1;
362                            if( i>1 )
363                            {
364                                v = 0.0;
365                                for(i_=1; i_<=im1;i_++)
366                                {
367                                    v += ludcmp[i,i_]*work0[i_];
368                                }
369                            }
370                            else
371                            {
372                                v = 0;
373                            }
374                            work0[i] = work0[i]+v;
375                        }
376                    }
377                    else
378                    {
379                       
380                        //
381                        // Multiply by L'
382                        //
383                        for(i=1; i<=n; i++)
384                        {
385                            ip1 = i+1;
386                            v = 0.0;
387                            for(i_=ip1; i_<=n;i_++)
388                            {
389                                v += ludcmp[i_,i]*work0[i_];
390                            }
391                            work0[i] = work0[i]+v;
392                        }
393                       
394                        //
395                        // Multiply by U'
396                        //
397                        for(i=n; i>=1; i--)
398                        {
399                            v = 0.0;
400                            for(i_=1; i_<=i;i_++)
401                            {
402                                v += ludcmp[i_,i]*work0[i_];
403                            }
404                            work0[i] = v;
405                        }
406                    }
407                }
408            }
409           
410            //
411            // Quick return if possible
412            //
413            rc = 0;
414            if( (double)(anorm)==(double)(0) )
415            {
416                return;
417            }
418           
419            //
420            // Estimate the norm of inv(A).
421            //
422            smlnum = AP.Math.MinRealNumber;
423            ainvnm = 0;
424            normin = false;
425            kase = 0;
426            while( true )
427            {
428                internalestimatenorm(n, ref work1, ref work0, ref iwork, ref ainvnm, ref kase);
429                if( kase==0 )
430                {
431                    break;
432                }
433                if( kase==kase1 )
434                {
435                   
436                    //
437                    // Multiply by inv(L).
438                    //
439                    trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref sl, !mupper, !mtrans, munit, normin, ref work2);
440                   
441                    //
442                    // Multiply by inv(U).
443                    //
444                    trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref su, mupper, !mtrans, !munit, normin, ref work3);
445                }
446                else
447                {
448                   
449                    //
450                    // Multiply by inv(U').
451                    //
452                    trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref su, mupper, mtrans, !munit, normin, ref work3);
453                   
454                    //
455                    // Multiply by inv(L').
456                    //
457                    trlinsolve.safesolvetriangular(ref ludcmp, n, ref work0, ref sl, !mupper, mtrans, munit, normin, ref work2);
458                }
459               
460                //
461                // Divide X by 1/(SL*SU) if doing so will not cause overflow.
462                //
463                ascale = sl*su;
464                normin = true;
465                if( (double)(ascale)!=(double)(1) )
466                {
467                    ix = 1;
468                    for(i=2; i<=n; i++)
469                    {
470                        if( (double)(Math.Abs(work0[i]))>(double)(Math.Abs(work0[ix])) )
471                        {
472                            ix = i;
473                        }
474                    }
475                    if( (double)(ascale)<(double)(Math.Abs(work0[ix])*smlnum) | (double)(ascale)==(double)(0) )
476                    {
477                        return;
478                    }
479                    for(i=1; i<=n; i++)
480                    {
481                        work0[i] = work0[i]/ascale;
482                    }
483                }
484            }
485           
486            //
487            // Compute the estimate of the reciprocal condition number.
488            //
489            if( (double)(ainvnm)!=(double)(0) )
490            {
491                rc = 1/ainvnm;
492                rc = rc/anorm;
493            }
494        }
495
496
497        private static void internalestimatenorm(int n,
498            ref double[] v,
499            ref double[] x,
500            ref int[] isgn,
501            ref double est,
502            ref int kase)
503        {
504            int itmax = 0;
505            int i = 0;
506            double t = 0;
507            bool flg = new bool();
508            int positer = 0;
509            int posj = 0;
510            int posjlast = 0;
511            int posjump = 0;
512            int posaltsgn = 0;
513            int posestold = 0;
514            int postemp = 0;
515            int i_ = 0;
516
517            itmax = 5;
518            posaltsgn = n+1;
519            posestold = n+2;
520            postemp = n+3;
521            positer = n+1;
522            posj = n+2;
523            posjlast = n+3;
524            posjump = n+4;
525            if( kase==0 )
526            {
527                v = new double[n+3+1];
528                x = new double[n+1];
529                isgn = new int[n+4+1];
530                t = (double)(1)/(double)(n);
531                for(i=1; i<=n; i++)
532                {
533                    x[i] = t;
534                }
535                kase = 1;
536                isgn[posjump] = 1;
537                return;
538            }
539           
540            //
541            //     ................ ENTRY   (JUMP = 1)
542            //     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
543            //
544            if( isgn[posjump]==1 )
545            {
546                if( n==1 )
547                {
548                    v[1] = x[1];
549                    est = Math.Abs(v[1]);
550                    kase = 0;
551                    return;
552                }
553                est = 0;
554                for(i=1; i<=n; i++)
555                {
556                    est = est+Math.Abs(x[i]);
557                }
558                for(i=1; i<=n; i++)
559                {
560                    if( (double)(x[i])>=(double)(0) )
561                    {
562                        x[i] = 1;
563                    }
564                    else
565                    {
566                        x[i] = -1;
567                    }
568                    isgn[i] = Math.Sign(x[i]);
569                }
570                kase = 2;
571                isgn[posjump] = 2;
572                return;
573            }
574           
575            //
576            //     ................ ENTRY   (JUMP = 2)
577            //     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
578            //
579            if( isgn[posjump]==2 )
580            {
581                isgn[posj] = 1;
582                for(i=2; i<=n; i++)
583                {
584                    if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
585                    {
586                        isgn[posj] = i;
587                    }
588                }
589                isgn[positer] = 2;
590               
591                //
592                // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
593                //
594                for(i=1; i<=n; i++)
595                {
596                    x[i] = 0;
597                }
598                x[isgn[posj]] = 1;
599                kase = 1;
600                isgn[posjump] = 3;
601                return;
602            }
603           
604            //
605            //     ................ ENTRY   (JUMP = 3)
606            //     X HAS BEEN OVERWRITTEN BY A*X.
607            //
608            if( isgn[posjump]==3 )
609            {
610                for(i_=1; i_<=n;i_++)
611                {
612                    v[i_] = x[i_];
613                }
614                v[posestold] = est;
615                est = 0;
616                for(i=1; i<=n; i++)
617                {
618                    est = est+Math.Abs(v[i]);
619                }
620                flg = false;
621                for(i=1; i<=n; i++)
622                {
623                    if( (double)(x[i])>=(double)(0) & isgn[i]<0 | (double)(x[i])<(double)(0) & isgn[i]>=0 )
624                    {
625                        flg = true;
626                    }
627                }
628               
629                //
630                // REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
631                // OR MAY BE CYCLING.
632                //
633                if( !flg | (double)(est)<=(double)(v[posestold]) )
634                {
635                    v[posaltsgn] = 1;
636                    for(i=1; i<=n; i++)
637                    {
638                        x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
639                        v[posaltsgn] = -v[posaltsgn];
640                    }
641                    kase = 1;
642                    isgn[posjump] = 5;
643                    return;
644                }
645                for(i=1; i<=n; i++)
646                {
647                    if( (double)(x[i])>=(double)(0) )
648                    {
649                        x[i] = 1;
650                        isgn[i] = 1;
651                    }
652                    else
653                    {
654                        x[i] = -1;
655                        isgn[i] = -1;
656                    }
657                }
658                kase = 2;
659                isgn[posjump] = 4;
660                return;
661            }
662           
663            //
664            //     ................ ENTRY   (JUMP = 4)
665            //     X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
666            //
667            if( isgn[posjump]==4 )
668            {
669                isgn[posjlast] = isgn[posj];
670                isgn[posj] = 1;
671                for(i=2; i<=n; i++)
672                {
673                    if( (double)(Math.Abs(x[i]))>(double)(Math.Abs(x[isgn[posj]])) )
674                    {
675                        isgn[posj] = i;
676                    }
677                }
678                if( (double)(x[isgn[posjlast]])!=(double)(Math.Abs(x[isgn[posj]])) & isgn[positer]<itmax )
679                {
680                    isgn[positer] = isgn[positer]+1;
681                    for(i=1; i<=n; i++)
682                    {
683                        x[i] = 0;
684                    }
685                    x[isgn[posj]] = 1;
686                    kase = 1;
687                    isgn[posjump] = 3;
688                    return;
689                }
690               
691                //
692                // ITERATION COMPLETE.  FINAL STAGE.
693                //
694                v[posaltsgn] = 1;
695                for(i=1; i<=n; i++)
696                {
697                    x[i] = v[posaltsgn]*(1+((double)(i-1))/((double)(n-1)));
698                    v[posaltsgn] = -v[posaltsgn];
699                }
700                kase = 1;
701                isgn[posjump] = 5;
702                return;
703            }
704           
705            //
706            //     ................ ENTRY   (JUMP = 5)
707            //     X HAS BEEN OVERWRITTEN BY A*X.
708            //
709            if( isgn[posjump]==5 )
710            {
711                v[postemp] = 0;
712                for(i=1; i<=n; i++)
713                {
714                    v[postemp] = v[postemp]+Math.Abs(x[i]);
715                }
716                v[postemp] = 2*v[postemp]/(3*n);
717                if( (double)(v[postemp])>(double)(est) )
718                {
719                    for(i_=1; i_<=n;i_++)
720                    {
721                        v[i_] = x[i_];
722                    }
723                    est = v[postemp];
724                }
725                kase = 0;
726                return;
727            }
728        }
729    }
730}
Note: See TracBrowser for help on using the repository browser.