Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/crcond.cs @ 2575

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

Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 26.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 crcond
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 cmatrixrcond1(ref AP.Complex[,] a,
47            int n)
48        {
49            double result = 0;
50            int i = 0;
51            AP.Complex[,] a1 = new AP.Complex[0,0];
52            int i_ = 0;
53            int i1_ = 0;
54
55            System.Diagnostics.Debug.Assert(n>=1, "CMatrixRCond1: N<1!");
56            a1 = new AP.Complex[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 = complexrcond1(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 CMatrixLU subroutine.
80            N           -   size of matrix A.
81
82        Result: 1/LowerBound(cond(A))
83        *************************************************************************/
84        public static double cmatrixlurcond1(ref AP.Complex[,] ludcmp,
85            int n)
86        {
87            double result = 0;
88            int i = 0;
89            AP.Complex[,] a1 = new AP.Complex[0,0];
90            int i_ = 0;
91            int i1_ = 0;
92
93            System.Diagnostics.Debug.Assert(n>=1, "CMatrixLURCond1: N<1!");
94            a1 = new AP.Complex[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 = complexrcond1lu(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 cmatrixrcondinf(ref AP.Complex[,] a,
122            int n)
123        {
124            double result = 0;
125            int i = 0;
126            AP.Complex[,] a1 = new AP.Complex[0,0];
127            int i_ = 0;
128            int i1_ = 0;
129
130            System.Diagnostics.Debug.Assert(n>=1, "CMatrixRCondInf: N<1!");
131            a1 = new AP.Complex[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 = complexrcondinf(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 CMatrixLU subroutine.
156            N       -   size of matrix A.
157
158        Result: 1/LowerBound(cond(A))
159        *************************************************************************/
160        public static double cmatrixlurcondinf(ref AP.Complex[,] ludcmp,
161            int n)
162        {
163            double result = 0;
164            int i = 0;
165            AP.Complex[,] a1 = new AP.Complex[0,0];
166            int i_ = 0;
167            int i1_ = 0;
168
169            System.Diagnostics.Debug.Assert(n>=1, "CMatrixLURCondInf: N<1!");
170            a1 = new AP.Complex[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 = complexrcondinflu(ref a1, n);
180            return result;
181        }
182
183
184        public static double complexrcond1(AP.Complex[,] 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 = (AP.Complex[,])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+AP.Math.AbsComplex(a[i,j]);
203                }
204                nrm = Math.Max(nrm, v);
205            }
206            clu.complexludecomposition(ref a, n, n, ref pivots);
207            internalestimatecomplexrcondlu(ref a, n, true, true, nrm, ref v);
208            result = v;
209            return result;
210        }
211
212
213        public static double complexrcond1lu(ref AP.Complex[,] lu,
214            int n)
215        {
216            double result = 0;
217            double v = 0;
218
219            internalestimatecomplexrcondlu(ref lu, n, true, false, 0, ref v);
220            result = v;
221            return result;
222        }
223
224
225        public static double complexrcondinf(AP.Complex[,] 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 = (AP.Complex[,])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+AP.Math.AbsComplex(a[i,j]);
244                }
245                nrm = Math.Max(nrm, v);
246            }
247            clu.complexludecomposition(ref a, n, n, ref pivots);
248            internalestimatecomplexrcondlu(ref a, n, false, true, nrm, ref v);
249            result = v;
250            return result;
251        }
252
253
254        public static double complexrcondinflu(ref AP.Complex[,] lu,
255            int n)
256        {
257            double result = 0;
258            double v = 0;
259
260            internalestimatecomplexrcondlu(ref lu, n, false, false, 0, ref v);
261            result = v;
262            return result;
263        }
264
265
266        public static void internalestimatecomplexrcondlu(ref AP.Complex[,] lu,
267            int n,
268            bool onenorm,
269            bool isanormprovided,
270            double anorm,
271            ref double rcond)
272        {
273            AP.Complex[] cwork1 = new AP.Complex[0];
274            AP.Complex[] cwork2 = new AP.Complex[0];
275            AP.Complex[] cwork3 = new AP.Complex[0];
276            AP.Complex[] cwork4 = new AP.Complex[0];
277            int[] isave = new int[0];
278            double[] rsave = new double[0];
279            int kase = 0;
280            int kase1 = 0;
281            double ainvnm = 0;
282            double smlnum = 0;
283            bool cw = new bool();
284            AP.Complex v = 0;
285            int i = 0;
286            int i_ = 0;
287
288            if( n<=0 )
289            {
290                return;
291            }
292            cwork1 = new AP.Complex[n+1];
293            cwork2 = new AP.Complex[n+1];
294            cwork3 = new AP.Complex[n+1];
295            cwork4 = new AP.Complex[n+1];
296            isave = new int[4+1];
297            rsave = new double[3+1];
298            rcond = 0;
299            if( n==0 )
300            {
301                rcond = 1;
302                return;
303            }
304            smlnum = AP.Math.MinRealNumber;
305           
306            //
307            // Estimate the norm of inv(A).
308            //
309            if( !isanormprovided )
310            {
311                anorm = 0;
312                if( onenorm )
313                {
314                    kase1 = 1;
315                }
316                else
317                {
318                    kase1 = 2;
319                }
320                kase = 0;
321                do
322                {
323                    internalcomplexrcondestimatenorm(n, ref cwork4, ref cwork1, ref anorm, ref kase, ref isave, ref rsave);
324                    if( kase!=0 )
325                    {
326                        if( kase==kase1 )
327                        {
328                           
329                            //
330                            // Multiply by U
331                            //
332                            for(i=1; i<=n; i++)
333                            {
334                                v = 0.0;
335                                for(i_=i; i_<=n;i_++)
336                                {
337                                    v += lu[i,i_]*cwork1[i_];
338                                }
339                                cwork1[i] = v;
340                            }
341                           
342                            //
343                            // Multiply by L
344                            //
345                            for(i=n; i>=1; i--)
346                            {
347                                v = 0;
348                                if( i>1 )
349                                {
350                                    v = 0.0;
351                                    for(i_=1; i_<=i-1;i_++)
352                                    {
353                                        v += lu[i,i_]*cwork1[i_];
354                                    }
355                                }
356                                cwork1[i] = v+cwork1[i];
357                            }
358                        }
359                        else
360                        {
361                           
362                            //
363                            // Multiply by L'
364                            //
365                            for(i=1; i<=n; i++)
366                            {
367                                cwork2[i] = 0;
368                            }
369                            for(i=1; i<=n; i++)
370                            {
371                                v = cwork1[i];
372                                if( i>1 )
373                                {
374                                    for(i_=1; i_<=i-1;i_++)
375                                    {
376                                        cwork2[i_] = cwork2[i_] + v*AP.Math.Conj(lu[i,i_]);
377                                    }
378                                }
379                                cwork2[i] = cwork2[i]+v;
380                            }
381                           
382                            //
383                            // Multiply by U'
384                            //
385                            for(i=1; i<=n; i++)
386                            {
387                                cwork1[i] = 0;
388                            }
389                            for(i=1; i<=n; i++)
390                            {
391                                v = cwork2[i];
392                                for(i_=i; i_<=n;i_++)
393                                {
394                                    cwork1[i_] = cwork1[i_] + v*AP.Math.Conj(lu[i,i_]);
395                                }
396                            }
397                        }
398                    }
399                }
400                while( kase!=0 );
401            }
402           
403            //
404            // Quick return if possible
405            //
406            if( (double)(anorm)==(double)(0) )
407            {
408                return;
409            }
410           
411            //
412            // Estimate the norm of inv(A).
413            //
414            ainvnm = 0;
415            if( onenorm )
416            {
417                kase1 = 1;
418            }
419            else
420            {
421                kase1 = 2;
422            }
423            kase = 0;
424            do
425            {
426                internalcomplexrcondestimatenorm(n, ref cwork4, ref cwork1, ref ainvnm, ref kase, ref isave, ref rsave);
427                if( kase!=0 )
428                {
429                    if( kase==kase1 )
430                    {
431                       
432                        //
433                        // Multiply by inv(L).
434                        //
435                        cw = ctrlinsolve.complexsafesolvetriangular(ref lu, n, ref cwork1, false, 0, true, ref cwork2, ref cwork3);
436                        if( !cw )
437                        {
438                            rcond = 0;
439                            return;
440                        }
441                       
442                        //
443                        // Multiply by inv(U).
444                        //
445                        cw = ctrlinsolve.complexsafesolvetriangular(ref lu, n, ref cwork1, true, 0, false, ref cwork2, ref cwork3);
446                        if( !cw )
447                        {
448                            rcond = 0;
449                            return;
450                        }
451                    }
452                    else
453                    {
454                       
455                        //
456                        // Multiply by inv(U').
457                        //
458                        cw = ctrlinsolve.complexsafesolvetriangular(ref lu, n, ref cwork1, true, 2, false, ref cwork2, ref cwork3);
459                        if( !cw )
460                        {
461                            rcond = 0;
462                            return;
463                        }
464                       
465                        //
466                        // Multiply by inv(L').
467                        //
468                        cw = ctrlinsolve.complexsafesolvetriangular(ref lu, n, ref cwork1, false, 2, true, ref cwork2, ref cwork3);
469                        if( !cw )
470                        {
471                            rcond = 0;
472                            return;
473                        }
474                    }
475                }
476            }
477            while( kase!=0 );
478           
479            //
480            // Compute the estimate of the reciprocal condition number.
481            //
482            if( (double)(ainvnm)!=(double)(0) )
483            {
484                rcond = 1/ainvnm;
485                rcond = rcond/anorm;
486            }
487        }
488
489
490        private static void internalcomplexrcondestimatenorm(int n,
491            ref AP.Complex[] v,
492            ref AP.Complex[] x,
493            ref double est,
494            ref int kase,
495            ref int[] isave,
496            ref double[] rsave)
497        {
498            int itmax = 0;
499            int i = 0;
500            int iter = 0;
501            int j = 0;
502            int jlast = 0;
503            int jump = 0;
504            double absxi = 0;
505            double altsgn = 0;
506            double estold = 0;
507            double safmin = 0;
508            double temp = 0;
509            int i_ = 0;
510
511           
512            //
513            //Executable Statements ..
514            //
515            itmax = 5;
516            safmin = AP.Math.MinRealNumber;
517            if( kase==0 )
518            {
519                for(i=1; i<=n; i++)
520                {
521                    x[i] = (double)(1)/(double)(n);
522                }
523                kase = 1;
524                jump = 1;
525                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
526                return;
527            }
528            internalcomplexrcondloadall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
529           
530            //
531            // ENTRY   (JUMP = 1)
532            // FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
533            //
534            if( jump==1 )
535            {
536                if( n==1 )
537                {
538                    v[1] = x[1];
539                    est = AP.Math.AbsComplex(v[1]);
540                    kase = 0;
541                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
542                    return;
543                }
544                est = internalcomplexrcondscsum1(ref x, n);
545                for(i=1; i<=n; i++)
546                {
547                    absxi = AP.Math.AbsComplex(x[i]);
548                    if( (double)(absxi)>(double)(safmin) )
549                    {
550                        x[i] = x[i]/absxi;
551                    }
552                    else
553                    {
554                        x[i] = 1;
555                    }
556                }
557                kase = 2;
558                jump = 2;
559                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
560                return;
561            }
562           
563            //
564            // ENTRY   (JUMP = 2)
565            // FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
566            //
567            if( jump==2 )
568            {
569                j = internalcomplexrcondicmax1(ref x, n);
570                iter = 2;
571               
572                //
573                // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
574                //
575                for(i=1; i<=n; i++)
576                {
577                    x[i] = 0;
578                }
579                x[j] = 1;
580                kase = 1;
581                jump = 3;
582                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
583                return;
584            }
585           
586            //
587            // ENTRY   (JUMP = 3)
588            // X HAS BEEN OVERWRITTEN BY A*X.
589            //
590            if( jump==3 )
591            {
592                for(i_=1; i_<=n;i_++)
593                {
594                    v[i_] = x[i_];
595                }
596                estold = est;
597                est = internalcomplexrcondscsum1(ref v, n);
598               
599                //
600                // TEST FOR CYCLING.
601                //
602                if( (double)(est)<=(double)(estold) )
603                {
604                   
605                    //
606                    // ITERATION COMPLETE.  FINAL STAGE.
607                    //
608                    altsgn = 1;
609                    for(i=1; i<=n; i++)
610                    {
611                        x[i] = altsgn*(1+((double)(i-1))/((double)(n-1)));
612                        altsgn = -altsgn;
613                    }
614                    kase = 1;
615                    jump = 5;
616                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
617                    return;
618                }
619                for(i=1; i<=n; i++)
620                {
621                    absxi = AP.Math.AbsComplex(x[i]);
622                    if( (double)(absxi)>(double)(safmin) )
623                    {
624                        x[i] = x[i]/absxi;
625                    }
626                    else
627                    {
628                        x[i] = 1;
629                    }
630                }
631                kase = 2;
632                jump = 4;
633                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
634                return;
635            }
636           
637            //
638            // ENTRY   (JUMP = 4)
639            // X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
640            //
641            if( jump==4 )
642            {
643                jlast = j;
644                j = internalcomplexrcondicmax1(ref x, n);
645                if( (double)(AP.Math.AbsComplex(x[jlast]))!=(double)(AP.Math.AbsComplex(x[j])) & iter<itmax )
646                {
647                    iter = iter+1;
648                   
649                    //
650                    // MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
651                    //
652                    for(i=1; i<=n; i++)
653                    {
654                        x[i] = 0;
655                    }
656                    x[j] = 1;
657                    kase = 1;
658                    jump = 3;
659                    internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
660                    return;
661                }
662               
663                //
664                // ITERATION COMPLETE.  FINAL STAGE.
665                //
666                altsgn = 1;
667                for(i=1; i<=n; i++)
668                {
669                    x[i] = altsgn*(1+((double)(i-1))/((double)(n-1)));
670                    altsgn = -altsgn;
671                }
672                kase = 1;
673                jump = 5;
674                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
675                return;
676            }
677           
678            //
679            // ENTRY   (JUMP = 5)
680            // X HAS BEEN OVERWRITTEN BY A*X.
681            //
682            if( jump==5 )
683            {
684                temp = 2*(internalcomplexrcondscsum1(ref x, n)/(3*n));
685                if( (double)(temp)>(double)(est) )
686                {
687                    for(i_=1; i_<=n;i_++)
688                    {
689                        v[i_] = x[i_];
690                    }
691                    est = temp;
692                }
693                kase = 0;
694                internalcomplexrcondsaveall(ref isave, ref rsave, ref i, ref iter, ref j, ref jlast, ref jump, ref absxi, ref altsgn, ref estold, ref temp);
695                return;
696            }
697        }
698
699
700        private static double internalcomplexrcondscsum1(ref AP.Complex[] x,
701            int n)
702        {
703            double result = 0;
704            int i = 0;
705
706            result = 0;
707            for(i=1; i<=n; i++)
708            {
709                result = result+AP.Math.AbsComplex(x[i]);
710            }
711            return result;
712        }
713
714
715        private static int internalcomplexrcondicmax1(ref AP.Complex[] x,
716            int n)
717        {
718            int result = 0;
719            int i = 0;
720            double m = 0;
721
722            result = 1;
723            m = AP.Math.AbsComplex(x[1]);
724            for(i=2; i<=n; i++)
725            {
726                if( (double)(AP.Math.AbsComplex(x[i]))>(double)(m) )
727                {
728                    result = i;
729                    m = AP.Math.AbsComplex(x[i]);
730                }
731            }
732            return result;
733        }
734
735
736        private static void internalcomplexrcondsaveall(ref int[] isave,
737            ref double[] rsave,
738            ref int i,
739            ref int iter,
740            ref int j,
741            ref int jlast,
742            ref int jump,
743            ref double absxi,
744            ref double altsgn,
745            ref double estold,
746            ref double temp)
747        {
748            isave[0] = i;
749            isave[1] = iter;
750            isave[2] = j;
751            isave[3] = jlast;
752            isave[4] = jump;
753            rsave[0] = absxi;
754            rsave[1] = altsgn;
755            rsave[2] = estold;
756            rsave[3] = temp;
757        }
758
759
760        private static void internalcomplexrcondloadall(ref int[] isave,
761            ref double[] rsave,
762            ref int i,
763            ref int iter,
764            ref int j,
765            ref int jlast,
766            ref int jump,
767            ref double absxi,
768            ref double altsgn,
769            ref double estold,
770            ref double temp)
771        {
772            i = isave[0];
773            iter = isave[1];
774            j = isave[2];
775            jlast = isave[3];
776            jump = isave[4];
777            absxi = rsave[0];
778            altsgn = rsave[1];
779            estold = rsave[2];
780            temp = rsave[3];
781        }
782    }
783}
Note: See TracBrowser for help on using the repository browser.