Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/densesolver.cs @ 2578

Last change on this file since 2578 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.0 KB
Line 
1/*************************************************************************
2Copyright (c) 2007-2008, Sergey Bochkanov (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 densesolver
26    {
27        public struct densesolverreport
28        {
29            public double r1;
30            public double rinf;
31        };
32
33
34        public struct densesolverlsreport
35        {
36            public double r2;
37            public double[,] cx;
38            public int n;
39            public int k;
40        };
41
42
43
44
45        /*************************************************************************
46        Dense solver.
47
48        This  subroutine  solves  a  system  A*X=B,  where A is NxN non-denegerate
49        real matrix, X and B are NxM real matrices.
50
51        Additional features include:
52        * automatic detection of degenerate cases
53        * iterative improvement
54
55        INPUT PARAMETERS
56            A       -   array[0..N-1,0..N-1], system matrix
57            N       -   size of A
58            B       -   array[0..N-1,0..M-1], right part
59            M       -   size of right part
60           
61        OUTPUT PARAMETERS
62            Info    -   return code:
63                        * -3    if A is singular, or VERY close to singular.
64                                X is filled by zeros in such cases.
65                        * -1    if N<=0 or M<=0 was passed
66                        *  1    if task is solved (matrix A may be near  singular,
67                                check R1/RInf parameters for condition numbers).
68            Rep     -   solver report, see below for more info
69            X       -   array[0..N-1,0..M-1], it contains:
70                        * solution of A*X=B if A is non-singular (well-conditioned
71                          or ill-conditioned, but not very close to singular)
72                        * zeros,  if  A  is  singular  or  VERY  close to singular
73                          (in this case Info=-3).
74
75        SOLVER REPORT
76
77        Subroutine sets following fields of the Rep structure:
78        * R1        reciprocal of condition number: 1/cond(A), 1-norm.
79        * RInf      reciprocal of condition number: 1/cond(A), inf-norm.
80
81        SEE ALSO:
82            DenseSolverR() - solves A*x = b, where x and b are Nx1 matrices.
83
84          -- ALGLIB --
85             Copyright 24.08.2009 by Bochkanov Sergey
86        *************************************************************************/
87        public static void rmatrixsolvem(ref double[,] a,
88            int n,
89            ref double[,] b,
90            int m,
91            ref int info,
92            ref densesolverreport rep,
93            ref double[,] x)
94        {
95            int i = 0;
96            int j = 0;
97            int k = 0;
98            int rfs = 0;
99            int nrfs = 0;
100            int[] p = new int[0];
101            double[] xc = new double[0];
102            double[] y = new double[0];
103            double[] bc = new double[0];
104            double[] xa = new double[0];
105            double[] xb = new double[0];
106            double[] tx = new double[0];
107            double[,] da = new double[0,0];
108            double v = 0;
109            double verr = 0;
110            bool smallerr = new bool();
111            bool terminatenexttime = new bool();
112            int i_ = 0;
113
114           
115            //
116            // prepare: check inputs, allocate space...
117            //
118            if( n<=0 | m<=0 )
119            {
120                info = -1;
121                return;
122            }
123            da = new double[n, n];
124            x = new double[n, m];
125            y = new double[n];
126            xc = new double[n];
127            bc = new double[n];
128            tx = new double[n+1];
129            xa = new double[n+1];
130            xb = new double[n+1];
131           
132            //
133            // factorize matrix, test for exact/near singularity
134            //
135            for(i=0; i<=n-1; i++)
136            {
137                for(i_=0; i_<=n-1;i_++)
138                {
139                    da[i,i_] = a[i,i_];
140                }
141            }
142            lu.rmatrixlu(ref da, n, n, ref p);
143            rep.r1 = rcond.rmatrixlurcond1(ref da, n);
144            rep.rinf = rcond.rmatrixlurcondinf(ref da, n);
145            if( (double)(rep.r1)<(double)(10*AP.Math.MachineEpsilon) | (double)(rep.rinf)<(double)(10*AP.Math.MachineEpsilon) )
146            {
147                for(i=0; i<=n-1; i++)
148                {
149                    for(j=0; j<=m-1; j++)
150                    {
151                        x[i,j] = 0;
152                    }
153                }
154                rep.r1 = 0;
155                rep.rinf = 0;
156                info = -3;
157                return;
158            }
159            info = 1;
160           
161            //
162            // solve
163            //
164            for(k=0; k<=m-1; k++)
165            {
166               
167                //
168                // First, non-iterative part of solution process:
169                // * pivots
170                // * L*y = b
171                // * U*x = y
172                //
173                for(i_=0; i_<=n-1;i_++)
174                {
175                    bc[i_] = b[i_,k];
176                }
177                for(i=0; i<=n-1; i++)
178                {
179                    if( p[i]!=i )
180                    {
181                        v = bc[i];
182                        bc[i] = bc[p[i]];
183                        bc[p[i]] = v;
184                    }
185                }
186                y[0] = bc[0];
187                for(i=1; i<=n-1; i++)
188                {
189                    v = 0.0;
190                    for(i_=0; i_<=i-1;i_++)
191                    {
192                        v += da[i,i_]*y[i_];
193                    }
194                    y[i] = bc[i]-v;
195                }
196                xc[n-1] = y[n-1]/da[n-1,n-1];
197                for(i=n-2; i>=0; i--)
198                {
199                    v = 0.0;
200                    for(i_=i+1; i_<=n-1;i_++)
201                    {
202                        v += da[i,i_]*xc[i_];
203                    }
204                    xc[i] = (y[i]-v)/da[i,i];
205                }
206               
207                //
208                // Iterative improvement of xc:
209                // * calculate r = bc-A*xc using extra-precise dot product
210                // * solve A*y = r
211                // * update x:=x+r
212                //
213                // This cycle is executed until one of two things happens:
214                // 1. maximum number of iterations reached
215                // 2. last iteration decreased error to the lower limit
216                //
217                nrfs = densesolverrfsmax(n, rep.r1, rep.rinf);
218                terminatenexttime = false;
219                for(rfs=0; rfs<=nrfs-1; rfs++)
220                {
221                    if( terminatenexttime )
222                    {
223                        break;
224                    }
225                   
226                    //
227                    // generate right part
228                    //
229                    smallerr = true;
230                    for(i=0; i<=n-1; i++)
231                    {
232                        for(i_=0; i_<=n-1;i_++)
233                        {
234                            xa[i_] = a[i,i_];
235                        }
236                        xa[n] = -1;
237                        for(i_=0; i_<=n-1;i_++)
238                        {
239                            xb[i_] = xc[i_];
240                        }
241                        xb[n] = b[i,k];
242                        xblas.xdot(ref xa, ref xb, n+1, ref tx, ref v, ref verr);
243                        bc[i] = -v;
244                        smallerr = smallerr & (double)(Math.Abs(v))<(double)(4*verr);
245                    }
246                    if( smallerr )
247                    {
248                        terminatenexttime = true;
249                    }
250                   
251                    //
252                    // solve
253                    //
254                    for(i=0; i<=n-1; i++)
255                    {
256                        if( p[i]!=i )
257                        {
258                            v = bc[i];
259                            bc[i] = bc[p[i]];
260                            bc[p[i]] = v;
261                        }
262                    }
263                    y[0] = bc[0];
264                    for(i=1; i<=n-1; i++)
265                    {
266                        v = 0.0;
267                        for(i_=0; i_<=i-1;i_++)
268                        {
269                            v += da[i,i_]*y[i_];
270                        }
271                        y[i] = bc[i]-v;
272                    }
273                    tx[n-1] = y[n-1]/da[n-1,n-1];
274                    for(i=n-2; i>=0; i--)
275                    {
276                        v = 0.0;
277                        for(i_=i+1; i_<=n-1;i_++)
278                        {
279                            v += da[i,i_]*tx[i_];
280                        }
281                        tx[i] = (y[i]-v)/da[i,i];
282                    }
283                   
284                    //
285                    // update
286                    //
287                    for(i_=0; i_<=n-1;i_++)
288                    {
289                        xc[i_] = xc[i_] + tx[i_];
290                    }
291                }
292               
293                //
294                // Store xc
295                //
296                for(i_=0; i_<=n-1;i_++)
297                {
298                    x[i_,k] = xc[i_];
299                }
300            }
301        }
302
303
304        /*************************************************************************
305        Dense solver.
306
307        This subroutine finds solution of the linear system A*X=B with non-square,
308        possibly degenerate A.  System  is  solved in the least squares sense, and
309        general least squares solution  X = X0 + CX*y  which  minimizes |A*X-B| is
310        returned. If A is non-degenerate, solution in the  usual sense is returned
311
312        Additional features include:
313        * iterative improvement
314
315        INPUT PARAMETERS
316            A       -   array[0..NRows-1,0..NCols-1], system matrix
317            NRows   -   vertical size of A
318            NCols   -   horizontal size of A
319            B       -   array[0..NCols-1], right part
320            Threshold-  a number in [0,1]. Singular values  beyond  Threshold  are
321                        considered  zero.  Set  it to 0.0, if you don't understand
322                        what it means, so the solver will choose good value on its
323                        own.
324                       
325        OUTPUT PARAMETERS
326            Info    -   return code:
327                        * -4    SVD subroutine failed
328                        * -1    if NRows<=0 or NCols<=0 or Threshold<0 was passed
329                        *  1    if task is solved
330            Rep     -   solver report, see below for more info
331            X       -   array[0..N-1,0..M-1], it contains:
332                        * solution of A*X=B if A is non-singular (well-conditioned
333                          or ill-conditioned, but not very close to singular)
334                        * zeros,  if  A  is  singular  or  VERY  close to singular
335                          (in this case Info=-3).
336
337        SOLVER REPORT
338
339        Subroutine sets following fields of the Rep structure:
340        * R2        reciprocal of condition number: 1/cond(A), 2-norm.
341        * N         = NCols
342        * K         dim(Null(A))
343        * CX        array[0..N-1,0..K-1], kernel of A.
344                    Columns of CX store such vectors that A*CX[i]=0.
345
346          -- ALGLIB --
347             Copyright 24.08.2009 by Bochkanov Sergey
348        *************************************************************************/
349        public static void rmatrixsolvels(ref double[,] a,
350            int nrows,
351            int ncols,
352            ref double[] b,
353            double threshold,
354            ref int info,
355            ref densesolverlsreport rep,
356            ref double[] x)
357        {
358            double[] sv = new double[0];
359            double[,] u = new double[0,0];
360            double[,] vt = new double[0,0];
361            double[] rp = new double[0];
362            double[] utb = new double[0];
363            double[] sutb = new double[0];
364            double[] tmp = new double[0];
365            double[] ta = new double[0];
366            double[] tx = new double[0];
367            double[] buf = new double[0];
368            double[] w = new double[0];
369            int i = 0;
370            int j = 0;
371            int nsv = 0;
372            int kernelidx = 0;
373            double v = 0;
374            double verr = 0;
375            bool svdfailed = new bool();
376            bool zeroa = new bool();
377            int rfs = 0;
378            int nrfs = 0;
379            bool terminatenexttime = new bool();
380            bool smallerr = new bool();
381            int i_ = 0;
382
383            if( nrows<=0 | ncols<=0 | (double)(threshold)<(double)(0) )
384            {
385                info = -1;
386                return;
387            }
388            if( (double)(threshold)==(double)(0) )
389            {
390                threshold = 1000*AP.Math.MachineEpsilon;
391            }
392           
393            //
394            // Factorize A first
395            //
396            svdfailed = !svd.rmatrixsvd(a, nrows, ncols, 1, 2, 2, ref sv, ref u, ref vt);
397            zeroa = (double)(sv[0])==(double)(0);
398            if( svdfailed | zeroa )
399            {
400                if( svdfailed )
401                {
402                    info = -4;
403                }
404                else
405                {
406                    info = 1;
407                }
408                x = new double[ncols];
409                for(i=0; i<=ncols-1; i++)
410                {
411                    x[i] = 0;
412                }
413                rep.n = ncols;
414                rep.k = ncols;
415                rep.cx = new double[ncols, ncols];
416                for(i=0; i<=ncols-1; i++)
417                {
418                    for(j=0; j<=ncols-1; j++)
419                    {
420                        if( i==j )
421                        {
422                            rep.cx[i,j] = 1;
423                        }
424                        else
425                        {
426                            rep.cx[i,j] = 0;
427                        }
428                    }
429                }
430                rep.r2 = 0;
431                return;
432            }
433            nsv = Math.Min(ncols, nrows);
434            if( nsv==ncols )
435            {
436                rep.r2 = sv[nsv-1]/sv[0];
437            }
438            else
439            {
440                rep.r2 = 0;
441            }
442            rep.n = ncols;
443            info = 1;
444           
445            //
446            // Iterative improvement of xc combined with solution:
447            // 1. xc = 0
448            // 2. calculate r = bc-A*xc using extra-precise dot product
449            // 3. solve A*y = r
450            // 4. update x:=x+r
451            // 5. goto 2
452            //
453            // This cycle is executed until one of two things happens:
454            // 1. maximum number of iterations reached
455            // 2. last iteration decreased error to the lower limit
456            //
457            utb = new double[nsv];
458            sutb = new double[nsv];
459            x = new double[ncols];
460            tmp = new double[ncols];
461            ta = new double[ncols+1];
462            tx = new double[ncols+1];
463            buf = new double[ncols+1];
464            for(i=0; i<=ncols-1; i++)
465            {
466                x[i] = 0;
467            }
468            kernelidx = nsv;
469            for(i=0; i<=nsv-1; i++)
470            {
471                if( (double)(sv[i])<=(double)(threshold*sv[0]) )
472                {
473                    kernelidx = i;
474                    break;
475                }
476            }
477            rep.k = ncols-kernelidx;
478            nrfs = densesolverrfsmaxv2(ncols, rep.r2);
479            terminatenexttime = false;
480            rp = new double[nrows];
481            for(rfs=0; rfs<=nrfs; rfs++)
482            {
483                if( terminatenexttime )
484                {
485                    break;
486                }
487               
488                //
489                // calculate right part
490                //
491                if( rfs==0 )
492                {
493                    for(i_=0; i_<=nrows-1;i_++)
494                    {
495                        rp[i_] = b[i_];
496                    }
497                }
498                else
499                {
500                    smallerr = true;
501                    for(i=0; i<=nrows-1; i++)
502                    {
503                        for(i_=0; i_<=ncols-1;i_++)
504                        {
505                            ta[i_] = a[i,i_];
506                        }
507                        ta[ncols] = -1;
508                        for(i_=0; i_<=ncols-1;i_++)
509                        {
510                            tx[i_] = x[i_];
511                        }
512                        tx[ncols] = b[i];
513                        xblas.xdot(ref ta, ref tx, ncols+1, ref buf, ref v, ref verr);
514                        rp[i] = -v;
515                        smallerr = smallerr & (double)(Math.Abs(v))<(double)(4*verr);
516                    }
517                    if( smallerr )
518                    {
519                        terminatenexttime = true;
520                    }
521                }
522               
523                //
524                // solve A*dx = rp
525                //
526                for(i=0; i<=ncols-1; i++)
527                {
528                    tmp[i] = 0;
529                }
530                for(i=0; i<=nsv-1; i++)
531                {
532                    utb[i] = 0;
533                }
534                for(i=0; i<=nrows-1; i++)
535                {
536                    v = rp[i];
537                    for(i_=0; i_<=nsv-1;i_++)
538                    {
539                        utb[i_] = utb[i_] + v*u[i,i_];
540                    }
541                }
542                for(i=0; i<=nsv-1; i++)
543                {
544                    if( i<kernelidx )
545                    {
546                        sutb[i] = utb[i]/sv[i];
547                    }
548                    else
549                    {
550                        sutb[i] = 0;
551                    }
552                }
553                for(i=0; i<=nsv-1; i++)
554                {
555                    v = sutb[i];
556                    for(i_=0; i_<=ncols-1;i_++)
557                    {
558                        tmp[i_] = tmp[i_] + v*vt[i,i_];
559                    }
560                }
561               
562                //
563                // update x:  x:=x+dx
564                //
565                for(i_=0; i_<=ncols-1;i_++)
566                {
567                    x[i_] = x[i_] + tmp[i_];
568                }
569            }
570           
571            //
572            // fill CX
573            //
574            if( rep.k>0 )
575            {
576                rep.cx = new double[ncols, rep.k];
577                for(i=0; i<=rep.k-1; i++)
578                {
579                    for(i_=0; i_<=ncols-1;i_++)
580                    {
581                        rep.cx[i_,i] = vt[kernelidx+i,i_];
582                    }
583                }
584            }
585        }
586
587
588        /*************************************************************************
589        Dense solver.
590
591        Similar to RMatrixSolveM() but solves task with one right part  (where b/x
592        are vectors, not matrices).
593
594        See RMatrixSolveM()  description  for  more  information  about subroutine
595        parameters.
596
597          -- ALGLIB --
598             Copyright 24.08.2009 by Bochkanov Sergey
599        *************************************************************************/
600        public static void rmatrixsolve(ref double[,] a,
601            int n,
602            ref double[] b,
603            ref int info,
604            ref densesolverreport rep,
605            ref double[] x)
606        {
607            double[,] bm = new double[0,0];
608            double[,] xm = new double[0,0];
609            int i_ = 0;
610
611            if( n<=0 )
612            {
613                info = -1;
614                return;
615            }
616            bm = new double[n, 1];
617            for(i_=0; i_<=n-1;i_++)
618            {
619                bm[i_,0] = b[i_];
620            }
621            rmatrixsolvem(ref a, n, ref bm, 1, ref info, ref rep, ref xm);
622            x = new double[n];
623            for(i_=0; i_<=n-1;i_++)
624            {
625                x[i_] = xm[i_,0];
626            }
627        }
628
629
630        /*************************************************************************
631        Internal subroutine.
632        Returns maximum count of RFS iterations as function of:
633        1. machine epsilon
634        2. task size.
635        3. condition number
636        *************************************************************************/
637        private static int densesolverrfsmax(int n,
638            double r1,
639            double rinf)
640        {
641            int result = 0;
642
643            result = 2;
644            return result;
645        }
646
647
648        /*************************************************************************
649        Internal subroutine.
650        Returns maximum count of RFS iterations as function of:
651        1. machine epsilon
652        2. task size.
653        3. norm-2 condition number
654        *************************************************************************/
655        private static int densesolverrfsmaxv2(int n,
656            double r2)
657        {
658            int result = 0;
659
660            result = densesolverrfsmax(n, 0, 0);
661            return result;
662        }
663    }
664}
Note: See TracBrowser for help on using the repository browser.