Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/safesolve.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: 22.8 KB
Line 
1/*************************************************************************
2This file is a part of ALGLIB project.
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class safesolve
26    {
27        /*************************************************************************
28        Real implementation of CMatrixScaledTRSafeSolve
29
30          -- ALGLIB routine --
31             21.01.2010
32             Bochkanov Sergey
33        *************************************************************************/
34        public static bool rmatrixscaledtrsafesolve(ref double[,] a,
35            double sa,
36            int n,
37            ref double[] x,
38            bool isupper,
39            int trans,
40            bool isunit,
41            double maxgrowth)
42        {
43            bool result = new bool();
44            double lnmax = 0;
45            double nrmb = 0;
46            double nrmx = 0;
47            int i = 0;
48            AP.Complex alpha = 0;
49            AP.Complex beta = 0;
50            double vr = 0;
51            AP.Complex cx = 0;
52            double[] tmp = new double[0];
53            int i_ = 0;
54
55            System.Diagnostics.Debug.Assert(n>0, "RMatrixTRSafeSolve: incorrect N!");
56            System.Diagnostics.Debug.Assert(trans==0 | trans==1, "RMatrixTRSafeSolve: incorrect Trans!");
57            result = true;
58            lnmax = Math.Log(AP.Math.MaxRealNumber);
59           
60            //
61            // Quick return if possible
62            //
63            if( n<=0 )
64            {
65                return result;
66            }
67           
68            //
69            // Load norms: right part and X
70            //
71            nrmb = 0;
72            for(i=0; i<=n-1; i++)
73            {
74                nrmb = Math.Max(nrmb, Math.Abs(x[i]));
75            }
76            nrmx = 0;
77           
78            //
79            // Solve
80            //
81            tmp = new double[n];
82            result = true;
83            if( isupper & trans==0 )
84            {
85               
86                //
87                // U*x = b
88                //
89                for(i=n-1; i>=0; i--)
90                {
91                   
92                    //
93                    // Task is reduced to alpha*x[i] = beta
94                    //
95                    if( isunit )
96                    {
97                        alpha = sa;
98                    }
99                    else
100                    {
101                        alpha = a[i,i]*sa;
102                    }
103                    if( i<n-1 )
104                    {
105                        for(i_=i+1; i_<=n-1;i_++)
106                        {
107                            tmp[i_] = sa*a[i,i_];
108                        }
109                        vr = 0.0;
110                        for(i_=i+1; i_<=n-1;i_++)
111                        {
112                            vr += tmp[i_]*x[i_];
113                        }
114                        beta = x[i]-vr;
115                    }
116                    else
117                    {
118                        beta = x[i];
119                    }
120                   
121                    //
122                    // solve alpha*x[i] = beta
123                    //
124                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref cx);
125                    if( !result )
126                    {
127                        return result;
128                    }
129                    x[i] = cx.x;
130                }
131                return result;
132            }
133            if( !isupper & trans==0 )
134            {
135               
136                //
137                // L*x = b
138                //
139                for(i=0; i<=n-1; i++)
140                {
141                   
142                    //
143                    // Task is reduced to alpha*x[i] = beta
144                    //
145                    if( isunit )
146                    {
147                        alpha = sa;
148                    }
149                    else
150                    {
151                        alpha = a[i,i]*sa;
152                    }
153                    if( i>0 )
154                    {
155                        for(i_=0; i_<=i-1;i_++)
156                        {
157                            tmp[i_] = sa*a[i,i_];
158                        }
159                        vr = 0.0;
160                        for(i_=0; i_<=i-1;i_++)
161                        {
162                            vr += tmp[i_]*x[i_];
163                        }
164                        beta = x[i]-vr;
165                    }
166                    else
167                    {
168                        beta = x[i];
169                    }
170                   
171                    //
172                    // solve alpha*x[i] = beta
173                    //
174                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref cx);
175                    if( !result )
176                    {
177                        return result;
178                    }
179                    x[i] = cx.x;
180                }
181                return result;
182            }
183            if( isupper & trans==1 )
184            {
185               
186                //
187                // U^T*x = b
188                //
189                for(i=0; i<=n-1; i++)
190                {
191                   
192                    //
193                    // Task is reduced to alpha*x[i] = beta
194                    //
195                    if( isunit )
196                    {
197                        alpha = sa;
198                    }
199                    else
200                    {
201                        alpha = a[i,i]*sa;
202                    }
203                    beta = x[i];
204                   
205                    //
206                    // solve alpha*x[i] = beta
207                    //
208                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref cx);
209                    if( !result )
210                    {
211                        return result;
212                    }
213                    x[i] = cx.x;
214                   
215                    //
216                    // update the rest of right part
217                    //
218                    if( i<n-1 )
219                    {
220                        vr = cx.x;
221                        for(i_=i+1; i_<=n-1;i_++)
222                        {
223                            tmp[i_] = sa*a[i,i_];
224                        }
225                        for(i_=i+1; i_<=n-1;i_++)
226                        {
227                            x[i_] = x[i_] - vr*tmp[i_];
228                        }
229                    }
230                }
231                return result;
232            }
233            if( !isupper & trans==1 )
234            {
235               
236                //
237                // L^T*x = b
238                //
239                for(i=n-1; i>=0; i--)
240                {
241                   
242                    //
243                    // Task is reduced to alpha*x[i] = beta
244                    //
245                    if( isunit )
246                    {
247                        alpha = sa;
248                    }
249                    else
250                    {
251                        alpha = a[i,i]*sa;
252                    }
253                    beta = x[i];
254                   
255                    //
256                    // solve alpha*x[i] = beta
257                    //
258                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref cx);
259                    if( !result )
260                    {
261                        return result;
262                    }
263                    x[i] = cx.x;
264                   
265                    //
266                    // update the rest of right part
267                    //
268                    if( i>0 )
269                    {
270                        vr = cx.x;
271                        for(i_=0; i_<=i-1;i_++)
272                        {
273                            tmp[i_] = sa*a[i,i_];
274                        }
275                        for(i_=0; i_<=i-1;i_++)
276                        {
277                            x[i_] = x[i_] - vr*tmp[i_];
278                        }
279                    }
280                }
281                return result;
282            }
283            result = false;
284            return result;
285        }
286
287
288        /*************************************************************************
289        Internal subroutine for safe solution of
290
291            SA*op(A)=b
292           
293        where  A  is  NxN  upper/lower  triangular/unitriangular  matrix, op(A) is
294        either identity transform, transposition or Hermitian transposition, SA is
295        a scaling factor such that max(|SA*A[i,j]|) is close to 1.0 in magnutude.
296
297        This subroutine  limits  relative  growth  of  solution  (in inf-norm)  by
298        MaxGrowth,  returning  False  if  growth  exceeds MaxGrowth. Degenerate or
299        near-degenerate matrices are handled correctly (False is returned) as long
300        as MaxGrowth is significantly less than MaxRealNumber/norm(b).
301
302          -- ALGLIB routine --
303             21.01.2010
304             Bochkanov Sergey
305        *************************************************************************/
306        public static bool cmatrixscaledtrsafesolve(ref AP.Complex[,] a,
307            double sa,
308            int n,
309            ref AP.Complex[] x,
310            bool isupper,
311            int trans,
312            bool isunit,
313            double maxgrowth)
314        {
315            bool result = new bool();
316            double lnmax = 0;
317            double nrmb = 0;
318            double nrmx = 0;
319            int i = 0;
320            AP.Complex alpha = 0;
321            AP.Complex beta = 0;
322            AP.Complex vc = 0;
323            AP.Complex[] tmp = new AP.Complex[0];
324            int i_ = 0;
325
326            System.Diagnostics.Debug.Assert(n>0, "CMatrixTRSafeSolve: incorrect N!");
327            System.Diagnostics.Debug.Assert(trans==0 | trans==1 | trans==2, "CMatrixTRSafeSolve: incorrect Trans!");
328            result = true;
329            lnmax = Math.Log(AP.Math.MaxRealNumber);
330           
331            //
332            // Quick return if possible
333            //
334            if( n<=0 )
335            {
336                return result;
337            }
338           
339            //
340            // Load norms: right part and X
341            //
342            nrmb = 0;
343            for(i=0; i<=n-1; i++)
344            {
345                nrmb = Math.Max(nrmb, AP.Math.AbsComplex(x[i]));
346            }
347            nrmx = 0;
348           
349            //
350            // Solve
351            //
352            tmp = new AP.Complex[n];
353            result = true;
354            if( isupper & trans==0 )
355            {
356               
357                //
358                // U*x = b
359                //
360                for(i=n-1; i>=0; i--)
361                {
362                   
363                    //
364                    // Task is reduced to alpha*x[i] = beta
365                    //
366                    if( isunit )
367                    {
368                        alpha = sa;
369                    }
370                    else
371                    {
372                        alpha = a[i,i]*sa;
373                    }
374                    if( i<n-1 )
375                    {
376                        for(i_=i+1; i_<=n-1;i_++)
377                        {
378                            tmp[i_] = sa*a[i,i_];
379                        }
380                        vc = 0.0;
381                        for(i_=i+1; i_<=n-1;i_++)
382                        {
383                            vc += tmp[i_]*x[i_];
384                        }
385                        beta = x[i]-vc;
386                    }
387                    else
388                    {
389                        beta = x[i];
390                    }
391                   
392                    //
393                    // solve alpha*x[i] = beta
394                    //
395                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
396                    if( !result )
397                    {
398                        return result;
399                    }
400                    x[i] = vc;
401                }
402                return result;
403            }
404            if( !isupper & trans==0 )
405            {
406               
407                //
408                // L*x = b
409                //
410                for(i=0; i<=n-1; i++)
411                {
412                   
413                    //
414                    // Task is reduced to alpha*x[i] = beta
415                    //
416                    if( isunit )
417                    {
418                        alpha = sa;
419                    }
420                    else
421                    {
422                        alpha = a[i,i]*sa;
423                    }
424                    if( i>0 )
425                    {
426                        for(i_=0; i_<=i-1;i_++)
427                        {
428                            tmp[i_] = sa*a[i,i_];
429                        }
430                        vc = 0.0;
431                        for(i_=0; i_<=i-1;i_++)
432                        {
433                            vc += tmp[i_]*x[i_];
434                        }
435                        beta = x[i]-vc;
436                    }
437                    else
438                    {
439                        beta = x[i];
440                    }
441                   
442                    //
443                    // solve alpha*x[i] = beta
444                    //
445                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
446                    if( !result )
447                    {
448                        return result;
449                    }
450                    x[i] = vc;
451                }
452                return result;
453            }
454            if( isupper & trans==1 )
455            {
456               
457                //
458                // U^T*x = b
459                //
460                for(i=0; i<=n-1; i++)
461                {
462                   
463                    //
464                    // Task is reduced to alpha*x[i] = beta
465                    //
466                    if( isunit )
467                    {
468                        alpha = sa;
469                    }
470                    else
471                    {
472                        alpha = a[i,i]*sa;
473                    }
474                    beta = x[i];
475                   
476                    //
477                    // solve alpha*x[i] = beta
478                    //
479                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
480                    if( !result )
481                    {
482                        return result;
483                    }
484                    x[i] = vc;
485                   
486                    //
487                    // update the rest of right part
488                    //
489                    if( i<n-1 )
490                    {
491                        for(i_=i+1; i_<=n-1;i_++)
492                        {
493                            tmp[i_] = sa*a[i,i_];
494                        }
495                        for(i_=i+1; i_<=n-1;i_++)
496                        {
497                            x[i_] = x[i_] - vc*tmp[i_];
498                        }
499                    }
500                }
501                return result;
502            }
503            if( !isupper & trans==1 )
504            {
505               
506                //
507                // L^T*x = b
508                //
509                for(i=n-1; i>=0; i--)
510                {
511                   
512                    //
513                    // Task is reduced to alpha*x[i] = beta
514                    //
515                    if( isunit )
516                    {
517                        alpha = sa;
518                    }
519                    else
520                    {
521                        alpha = a[i,i]*sa;
522                    }
523                    beta = x[i];
524                   
525                    //
526                    // solve alpha*x[i] = beta
527                    //
528                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
529                    if( !result )
530                    {
531                        return result;
532                    }
533                    x[i] = vc;
534                   
535                    //
536                    // update the rest of right part
537                    //
538                    if( i>0 )
539                    {
540                        for(i_=0; i_<=i-1;i_++)
541                        {
542                            tmp[i_] = sa*a[i,i_];
543                        }
544                        for(i_=0; i_<=i-1;i_++)
545                        {
546                            x[i_] = x[i_] - vc*tmp[i_];
547                        }
548                    }
549                }
550                return result;
551            }
552            if( isupper & trans==2 )
553            {
554               
555                //
556                // U^H*x = b
557                //
558                for(i=0; i<=n-1; i++)
559                {
560                   
561                    //
562                    // Task is reduced to alpha*x[i] = beta
563                    //
564                    if( isunit )
565                    {
566                        alpha = sa;
567                    }
568                    else
569                    {
570                        alpha = AP.Math.Conj(a[i,i])*sa;
571                    }
572                    beta = x[i];
573                   
574                    //
575                    // solve alpha*x[i] = beta
576                    //
577                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
578                    if( !result )
579                    {
580                        return result;
581                    }
582                    x[i] = vc;
583                   
584                    //
585                    // update the rest of right part
586                    //
587                    if( i<n-1 )
588                    {
589                        for(i_=i+1; i_<=n-1;i_++)
590                        {
591                            tmp[i_] = sa*AP.Math.Conj(a[i,i_]);
592                        }
593                        for(i_=i+1; i_<=n-1;i_++)
594                        {
595                            x[i_] = x[i_] - vc*tmp[i_];
596                        }
597                    }
598                }
599                return result;
600            }
601            if( !isupper & trans==2 )
602            {
603               
604                //
605                // L^T*x = b
606                //
607                for(i=n-1; i>=0; i--)
608                {
609                   
610                    //
611                    // Task is reduced to alpha*x[i] = beta
612                    //
613                    if( isunit )
614                    {
615                        alpha = sa;
616                    }
617                    else
618                    {
619                        alpha = AP.Math.Conj(a[i,i])*sa;
620                    }
621                    beta = x[i];
622                   
623                    //
624                    // solve alpha*x[i] = beta
625                    //
626                    result = cbasicsolveandupdate(alpha, beta, lnmax, nrmb, maxgrowth, ref nrmx, ref vc);
627                    if( !result )
628                    {
629                        return result;
630                    }
631                    x[i] = vc;
632                   
633                    //
634                    // update the rest of right part
635                    //
636                    if( i>0 )
637                    {
638                        for(i_=0; i_<=i-1;i_++)
639                        {
640                            tmp[i_] = sa*AP.Math.Conj(a[i,i_]);
641                        }
642                        for(i_=0; i_<=i-1;i_++)
643                        {
644                            x[i_] = x[i_] - vc*tmp[i_];
645                        }
646                    }
647                }
648                return result;
649            }
650            result = false;
651            return result;
652        }
653
654
655        /*************************************************************************
656        complex basic solver-updater for reduced linear system
657
658            alpha*x[i] = beta
659
660        solves this equation and updates it in overlfow-safe manner (keeping track
661        of relative growth of solution).
662
663        Parameters:
664            Alpha   -   alpha
665            Beta    -   beta
666            LnMax   -   precomputed Ln(MaxRealNumber)
667            BNorm   -   inf-norm of b (right part of original system)
668            MaxGrowth-  maximum growth of norm(x) relative to norm(b)
669            XNorm   -   inf-norm of other components of X (which are already processed)
670                        it is updated by CBasicSolveAndUpdate.
671            X       -   solution
672
673          -- ALGLIB routine --
674             26.01.2009
675             Bochkanov Sergey
676        *************************************************************************/
677        private static bool cbasicsolveandupdate(AP.Complex alpha,
678            AP.Complex beta,
679            double lnmax,
680            double bnorm,
681            double maxgrowth,
682            ref double xnorm,
683            ref AP.Complex x)
684        {
685            bool result = new bool();
686            double v = 0;
687
688            result = false;
689            if( alpha==0 )
690            {
691                return result;
692            }
693            if( beta!=0 )
694            {
695               
696                //
697                // alpha*x[i]=beta
698                //
699                v = Math.Log(AP.Math.AbsComplex(beta))-Math.Log(AP.Math.AbsComplex(alpha));
700                if( (double)(v)>(double)(lnmax) )
701                {
702                    return result;
703                }
704                x = beta/alpha;
705            }
706            else
707            {
708               
709                //
710                // alpha*x[i]=0
711                //
712                x = 0;
713            }
714           
715            //
716            // update NrmX, test growth limit
717            //
718            xnorm = Math.Max(xnorm, AP.Math.AbsComplex(x));
719            if( (double)(xnorm)>(double)(maxgrowth*bnorm) )
720            {
721                return result;
722            }
723            result = true;
724            return result;
725        }
726    }
727}
Note: See TracBrowser for help on using the repository browser.