Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/srcond.cs @ 2563

Last change on this file since 2563 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: 22.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 srcond
32    {
33        /*************************************************************************
34        Condition number estimate of a symmetric matrix
35
36        The algorithm calculates a lower bound of the condition number. In this
37        case, the algorithm does not return a lower bound of the condition number,
38        but an inverse number (to avoid an overflow in case of a singular matrix).
39
40        It should be noted that 1-norm and inf-norm condition numbers of symmetric
41        matrices are equal, so the algorithm doesn't take into account the
42        differences between these types of norms.
43
44        Input parameters:
45            A       -   symmetric definite matrix which is given by its upper or
46                        lower triangle depending on IsUpper.
47                        Array with elements [0..N-1, 0..N-1].
48            N       -   size of matrix A.
49            IsUpper -   storage format.
50
51        Result:
52            1/LowerBound(cond(A))
53        *************************************************************************/
54        public static double smatrixrcond(ref double[,] a,
55            int n,
56            bool isupper)
57        {
58            double result = 0;
59            int i = 0;
60            int j = 0;
61            double[,] a1 = new double[0,0];
62
63            a1 = new double[n+1, n+1];
64            for(i=1; i<=n; i++)
65            {
66                if( isupper )
67                {
68                    for(j=i; j<=n; j++)
69                    {
70                        a1[i,j] = a[i-1,j-1];
71                    }
72                }
73                else
74                {
75                    for(j=1; j<=i; j++)
76                    {
77                        a1[i,j] = a[i-1,j-1];
78                    }
79                }
80            }
81            result = rcondsymmetric(a1, n, isupper);
82            return result;
83        }
84
85
86        /*************************************************************************
87        Condition number estimate of a matrix given by LDLT-decomposition
88
89        The algorithm calculates a lower bound of the condition number. In this
90        case, the algorithm does not return a lower bound of the condition number,
91        but an inverse number (to avoid an overflow in case of a singular matrix).
92
93        It should be noted that 1-norm and inf-norm condition numbers of symmetric
94        matrices are equal, so the algorithm doesn't take into account the
95        differences between these types of norms.
96
97        Input parameters:
98            L       -   LDLT-decomposition of matrix A given by the upper or lower
99                        triangle depending on IsUpper.
100                        Output of SMatrixLDLT subroutine.
101            Pivots  -   table of permutations which were made during LDLT-decomposition,
102                        Output of SMatrixLDLT subroutine.
103            N       -   size of matrix A.
104            IsUpper -   storage format.
105
106        Result:
107            1/LowerBound(cond(A))
108        *************************************************************************/
109        public static double smatrixldltrcond(ref double[,] l,
110            ref int[] pivots,
111            int n,
112            bool isupper)
113        {
114            double result = 0;
115            int i = 0;
116            int j = 0;
117            double[,] l1 = new double[0,0];
118            int[] p1 = new int[0];
119
120            l1 = new double[n+1, n+1];
121            for(i=1; i<=n; i++)
122            {
123                if( isupper )
124                {
125                    for(j=i; j<=n; j++)
126                    {
127                        l1[i,j] = l[i-1,j-1];
128                    }
129                }
130                else
131                {
132                    for(j=1; j<=i; j++)
133                    {
134                        l1[i,j] = l[i-1,j-1];
135                    }
136                }
137            }
138            p1 = new int[n+1];
139            for(i=1; i<=n; i++)
140            {
141                if( pivots[i-1]>=0 )
142                {
143                    p1[i] = pivots[i-1]+1;
144                }
145                else
146                {
147                    p1[i] = -(pivots[i-1]+n+1);
148                }
149            }
150            result = rcondldlt(ref l1, ref p1, n, isupper);
151            return result;
152        }
153
154
155        public static double rcondsymmetric(double[,] a,
156            int n,
157            bool isupper)
158        {
159            double result = 0;
160            int i = 0;
161            int j = 0;
162            int im = 0;
163            int jm = 0;
164            double v = 0;
165            double nrm = 0;
166            int[] pivots = new int[0];
167
168            a = (double[,])a.Clone();
169
170            nrm = 0;
171            for(j=1; j<=n; j++)
172            {
173                v = 0;
174                for(i=1; i<=n; i++)
175                {
176                    im = i;
177                    jm = j;
178                    if( isupper & j<i )
179                    {
180                        im = j;
181                        jm = i;
182                    }
183                    if( !isupper & j>i )
184                    {
185                        im = j;
186                        jm = i;
187                    }
188                    v = v+Math.Abs(a[im,jm]);
189                }
190                nrm = Math.Max(nrm, v);
191            }
192            ldlt.ldltdecomposition(ref a, n, isupper, ref pivots);
193            internalldltrcond(ref a, ref pivots, n, isupper, true, nrm, ref v);
194            result = v;
195            return result;
196        }
197
198
199        public static double rcondldlt(ref double[,] l,
200            ref int[] pivots,
201            int n,
202            bool isupper)
203        {
204            double result = 0;
205            double v = 0;
206
207            internalldltrcond(ref l, ref pivots, n, isupper, false, 0, ref v);
208            result = v;
209            return result;
210        }
211
212
213        public static void internalldltrcond(ref double[,] l,
214            ref int[] pivots,
215            int n,
216            bool isupper,
217            bool isnormprovided,
218            double anorm,
219            ref double rcond)
220        {
221            int i = 0;
222            int ix = 0;
223            int kase = 0;
224            int k = 0;
225            int km1 = 0;
226            int km2 = 0;
227            int kp1 = 0;
228            int kp2 = 0;
229            double ainvnm = 0;
230            double[] work0 = new double[0];
231            double[] work1 = new double[0];
232            double[] work2 = new double[0];
233            int[] iwork = new int[0];
234            double v = 0;
235            int i_ = 0;
236
237            System.Diagnostics.Debug.Assert(n>=0);
238           
239            //
240            // Check that the diagonal matrix D is nonsingular.
241            //
242            rcond = 0;
243            if( isupper )
244            {
245                for(i=n; i>=1; i--)
246                {
247                    if( pivots[i]>0 & (double)(l[i,i])==(double)(0) )
248                    {
249                        return;
250                    }
251                }
252            }
253            else
254            {
255                for(i=1; i<=n; i++)
256                {
257                    if( pivots[i]>0 & (double)(l[i,i])==(double)(0) )
258                    {
259                        return;
260                    }
261                }
262            }
263           
264            //
265            // Estimate the norm of A.
266            //
267            if( !isnormprovided )
268            {
269                kase = 0;
270                anorm = 0;
271                while( true )
272                {
273                    estnorm.iterativeestimate1norm(n, ref work1, ref work0, ref iwork, ref anorm, ref kase);
274                    if( kase==0 )
275                    {
276                        break;
277                    }
278                    if( isupper )
279                    {
280                       
281                        //
282                        // Multiply by U'
283                        //
284                        k = n;
285                        while( k>=1 )
286                        {
287                            if( pivots[k]>0 )
288                            {
289                               
290                                //
291                                // P(k)
292                                //
293                                v = work0[k];
294                                work0[k] = work0[pivots[k]];
295                                work0[pivots[k]] = v;
296                               
297                                //
298                                // U(k)
299                                //
300                                km1 = k-1;
301                                v = 0.0;
302                                for(i_=1; i_<=km1;i_++)
303                                {
304                                    v += work0[i_]*l[i_,k];
305                                }
306                                work0[k] = work0[k]+v;
307                               
308                                //
309                                // Next k
310                                //
311                                k = k-1;
312                            }
313                            else
314                            {
315                               
316                                //
317                                // P(k)
318                                //
319                                v = work0[k-1];
320                                work0[k-1] = work0[-pivots[k-1]];
321                                work0[-pivots[k-1]] = v;
322                               
323                                //
324                                // U(k)
325                                //
326                                km1 = k-1;
327                                km2 = k-2;
328                                v = 0.0;
329                                for(i_=1; i_<=km2;i_++)
330                                {
331                                    v += work0[i_]*l[i_,km1];
332                                }
333                                work0[km1] = work0[km1]+v;
334                                v = 0.0;
335                                for(i_=1; i_<=km2;i_++)
336                                {
337                                    v += work0[i_]*l[i_,k];
338                                }
339                                work0[k] = work0[k]+v;
340                               
341                                //
342                                // Next k
343                                //
344                                k = k-2;
345                            }
346                        }
347                       
348                        //
349                        // Multiply by D
350                        //
351                        k = n;
352                        while( k>=1 )
353                        {
354                            if( pivots[k]>0 )
355                            {
356                                work0[k] = work0[k]*l[k,k];
357                                k = k-1;
358                            }
359                            else
360                            {
361                                v = work0[k-1];
362                                work0[k-1] = l[k-1,k-1]*work0[k-1]+l[k-1,k]*work0[k];
363                                work0[k] = l[k-1,k]*v+l[k,k]*work0[k];
364                                k = k-2;
365                            }
366                        }
367                       
368                        //
369                        // Multiply by U
370                        //
371                        k = 1;
372                        while( k<=n )
373                        {
374                            if( pivots[k]>0 )
375                            {
376                               
377                                //
378                                // U(k)
379                                //
380                                km1 = k-1;
381                                v = work0[k];
382                                for(i_=1; i_<=km1;i_++)
383                                {
384                                    work0[i_] = work0[i_] + v*l[i_,k];
385                                }
386                               
387                                //
388                                // P(k)
389                                //
390                                v = work0[k];
391                                work0[k] = work0[pivots[k]];
392                                work0[pivots[k]] = v;
393                               
394                                //
395                                // Next k
396                                //
397                                k = k+1;
398                            }
399                            else
400                            {
401                               
402                                //
403                                // U(k)
404                                //
405                                km1 = k-1;
406                                kp1 = k+1;
407                                v = work0[k];
408                                for(i_=1; i_<=km1;i_++)
409                                {
410                                    work0[i_] = work0[i_] + v*l[i_,k];
411                                }
412                                v = work0[kp1];
413                                for(i_=1; i_<=km1;i_++)
414                                {
415                                    work0[i_] = work0[i_] + v*l[i_,kp1];
416                                }
417                               
418                                //
419                                // P(k)
420                                //
421                                v = work0[k];
422                                work0[k] = work0[-pivots[k]];
423                                work0[-pivots[k]] = v;
424                               
425                                //
426                                // Next k
427                                //
428                                k = k+2;
429                            }
430                        }
431                    }
432                    else
433                    {
434                       
435                        //
436                        // Multiply by L'
437                        //
438                        k = 1;
439                        while( k<=n )
440                        {
441                            if( pivots[k]>0 )
442                            {
443                               
444                                //
445                                // P(k)
446                                //
447                                v = work0[k];
448                                work0[k] = work0[pivots[k]];
449                                work0[pivots[k]] = v;
450                               
451                                //
452                                // L(k)
453                                //
454                                kp1 = k+1;
455                                v = 0.0;
456                                for(i_=kp1; i_<=n;i_++)
457                                {
458                                    v += work0[i_]*l[i_,k];
459                                }
460                                work0[k] = work0[k]+v;
461                               
462                                //
463                                // Next k
464                                //
465                                k = k+1;
466                            }
467                            else
468                            {
469                               
470                                //
471                                // P(k)
472                                //
473                                v = work0[k+1];
474                                work0[k+1] = work0[-pivots[k+1]];
475                                work0[-pivots[k+1]] = v;
476                               
477                                //
478                                // L(k)
479                                //
480                                kp1 = k+1;
481                                kp2 = k+2;
482                                v = 0.0;
483                                for(i_=kp2; i_<=n;i_++)
484                                {
485                                    v += work0[i_]*l[i_,k];
486                                }
487                                work0[k] = work0[k]+v;
488                                v = 0.0;
489                                for(i_=kp2; i_<=n;i_++)
490                                {
491                                    v += work0[i_]*l[i_,kp1];
492                                }
493                                work0[kp1] = work0[kp1]+v;
494                               
495                                //
496                                // Next k
497                                //
498                                k = k+2;
499                            }
500                        }
501                       
502                        //
503                        // Multiply by D
504                        //
505                        k = n;
506                        while( k>=1 )
507                        {
508                            if( pivots[k]>0 )
509                            {
510                                work0[k] = work0[k]*l[k,k];
511                                k = k-1;
512                            }
513                            else
514                            {
515                                v = work0[k-1];
516                                work0[k-1] = l[k-1,k-1]*work0[k-1]+l[k,k-1]*work0[k];
517                                work0[k] = l[k,k-1]*v+l[k,k]*work0[k];
518                                k = k-2;
519                            }
520                        }
521                       
522                        //
523                        // Multiply by L
524                        //
525                        k = n;
526                        while( k>=1 )
527                        {
528                            if( pivots[k]>0 )
529                            {
530                               
531                                //
532                                // L(k)
533                                //
534                                kp1 = k+1;
535                                v = work0[k];
536                                for(i_=kp1; i_<=n;i_++)
537                                {
538                                    work0[i_] = work0[i_] + v*l[i_,k];
539                                }
540                               
541                                //
542                                // P(k)
543                                //
544                                v = work0[k];
545                                work0[k] = work0[pivots[k]];
546                                work0[pivots[k]] = v;
547                               
548                                //
549                                // Next k
550                                //
551                                k = k-1;
552                            }
553                            else
554                            {
555                               
556                                //
557                                // L(k)
558                                //
559                                kp1 = k+1;
560                                km1 = k-1;
561                                v = work0[k];
562                                for(i_=kp1; i_<=n;i_++)
563                                {
564                                    work0[i_] = work0[i_] + v*l[i_,k];
565                                }
566                                v = work0[km1];
567                                for(i_=kp1; i_<=n;i_++)
568                                {
569                                    work0[i_] = work0[i_] + v*l[i_,km1];
570                                }
571                               
572                                //
573                                // P(k)
574                                //
575                                v = work0[k];
576                                work0[k] = work0[-pivots[k]];
577                                work0[-pivots[k]] = v;
578                               
579                                //
580                                // Next k
581                                //
582                                k = k-2;
583                            }
584                        }
585                    }
586                }
587            }
588           
589            //
590            // Quick return if possible
591            //
592            rcond = 0;
593            if( n==0 )
594            {
595                rcond = 1;
596                return;
597            }
598            if( (double)(anorm)==(double)(0) )
599            {
600                return;
601            }
602           
603            //
604            // Estimate the 1-norm of inv(A).
605            //
606            kase = 0;
607            while( true )
608            {
609                estnorm.iterativeestimate1norm(n, ref work1, ref work0, ref iwork, ref ainvnm, ref kase);
610                if( kase==0 )
611                {
612                    break;
613                }
614                ssolve.solvesystemldlt(ref l, ref pivots, work0, n, isupper, ref work2);
615                for(i_=1; i_<=n;i_++)
616                {
617                    work0[i_] = work2[i_];
618                }
619            }
620           
621            //
622            // Compute the estimate of the reciprocal condition number.
623            //
624            if( (double)(ainvnm)!=(double)(0) )
625            {
626                v = 1/ainvnm;
627                rcond = v/anorm;
628            }
629        }
630    }
631}
Note: See TracBrowser for help on using the repository browser.