Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/ssolve.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: 34.7 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 ssolve
32    {
33        /*************************************************************************
34        Solving  a system  of linear equations  with a system matrix  given by its
35        LDLT decomposition
36
37        The algorithm solves systems with a square matrix only.
38
39        Input parameters:
40            A       -   LDLT decomposition of the matrix (the result of the
41                        SMatrixLDLT subroutine).
42            Pivots  -   row permutation table (the result of the SMatrixLDLT subroutine).
43            B       -   right side of a system.
44                        Array whose index ranges within [0..N-1].
45            N       -   size of matrix A.
46            IsUpper -   points to the triangle of matrix A in which the LDLT
47                        decomposition is stored.
48                        If IsUpper=True, the decomposition has the form of U*D*U',
49                        matrix U is stored in the upper triangle of  matrix A  (in
50                        that case, the lower triangle isn't used and isn't changed
51                        by the subroutine).
52                        Similarly, if IsUpper=False, the decomposition has the form
53                        of L*D*L' and the lower triangle stores matrix L.
54
55        Output parameters:
56            X       -   solution of a system.
57                        Array whose index ranges within [0..N-1].
58
59        Result:
60            True, if the matrix is not singular. X contains the solution.
61            False, if the matrix is singular (the determinant of matrix D is equal
62        to 0). In this case, X doesn't contain a solution.
63        *************************************************************************/
64        public static bool smatrixldltsolve(ref double[,] a,
65            ref int[] pivots,
66            double[] b,
67            int n,
68            bool isupper,
69            ref double[] x)
70        {
71            bool result = new bool();
72            int i = 0;
73            int j = 0;
74            int k = 0;
75            int kp = 0;
76            double ak = 0;
77            double akm1 = 0;
78            double akm1k = 0;
79            double bk = 0;
80            double bkm1 = 0;
81            double denom = 0;
82            double v = 0;
83            int i_ = 0;
84
85            b = (double[])b.Clone();
86
87           
88            //
89            // Quick return if possible
90            //
91            result = true;
92            if( n==0 )
93            {
94                return result;
95            }
96           
97            //
98            // Check that the diagonal matrix D is nonsingular
99            //
100            if( isupper )
101            {
102               
103                //
104                // Upper triangular storage: examine D from bottom to top
105                //
106                for(i=n-1; i>=0; i--)
107                {
108                    if( pivots[i]>=0 & (double)(a[i,i])==(double)(0) )
109                    {
110                        result = false;
111                        return result;
112                    }
113                }
114            }
115            else
116            {
117               
118                //
119                // Lower triangular storage: examine D from top to bottom.
120                //
121                for(i=0; i<=n-1; i++)
122                {
123                    if( pivots[i]>=0 & (double)(a[i,i])==(double)(0) )
124                    {
125                        result = false;
126                        return result;
127                    }
128                }
129            }
130           
131            //
132            // Solve Ax = b
133            //
134            if( isupper )
135            {
136               
137                //
138                // Solve A*X = B, where A = U*D*U'.
139                //
140                // First solve U*D*X = B, overwriting B with X.
141                //
142                // K+1 is the main loop index, decreasing from N to 1 in steps of
143                // 1 or 2, depending on the size of the diagonal blocks.
144                //
145                k = n-1;
146                while( k>=0 )
147                {
148                    if( pivots[k]>=0 )
149                    {
150                       
151                        //
152                        // 1 x 1 diagonal block
153                        //
154                        // Interchange rows K+1 and IPIV(K+1).
155                        //
156                        kp = pivots[k];
157                        if( kp!=k )
158                        {
159                            v = b[k];
160                            b[k] = b[kp];
161                            b[kp] = v;
162                        }
163                       
164                        //
165                        // Multiply by inv(U(K+1)), where U(K+1) is the transformation
166                        // stored in column K+1 of A.
167                        //
168                        v = b[k];
169                        for(i_=0; i_<=k-1;i_++)
170                        {
171                            b[i_] = b[i_] - v*a[i_,k];
172                        }
173                       
174                        //
175                        // Multiply by the inverse of the diagonal block.
176                        //
177                        b[k] = b[k]/a[k,k];
178                        k = k-1;
179                    }
180                    else
181                    {
182                       
183                        //
184                        // 2 x 2 diagonal block
185                        //
186                        // Interchange rows K+1-1 and -IPIV(K+1).
187                        //
188                        kp = pivots[k]+n;
189                        if( kp!=k-1 )
190                        {
191                            v = b[k-1];
192                            b[k-1] = b[kp];
193                            b[kp] = v;
194                        }
195                       
196                        //
197                        // Multiply by inv(U(K+1)), where U(K+1) is the transformation
198                        // stored in columns K+1-1 and K+1 of A.
199                        //
200                        v = b[k];
201                        for(i_=0; i_<=k-2;i_++)
202                        {
203                            b[i_] = b[i_] - v*a[i_,k];
204                        }
205                        v = b[k-1];
206                        for(i_=0; i_<=k-2;i_++)
207                        {
208                            b[i_] = b[i_] - v*a[i_,k-1];
209                        }
210                       
211                        //
212                        // Multiply by the inverse of the diagonal block.
213                        //
214                        akm1k = a[k-1,k];
215                        akm1 = a[k-1,k-1]/akm1k;
216                        ak = a[k,k]/akm1k;
217                        denom = akm1*ak-1;
218                        bkm1 = b[k-1]/akm1k;
219                        bk = b[k]/akm1k;
220                        b[k-1] = (ak*bkm1-bk)/denom;
221                        b[k] = (akm1*bk-bkm1)/denom;
222                        k = k-2;
223                    }
224                }
225               
226                //
227                // Next solve U'*X = B, overwriting B with X.
228                //
229                // K+1 is the main loop index, increasing from 1 to N in steps of
230                // 1 or 2, depending on the size of the diagonal blocks.
231                //
232                k = 0;
233                while( k<=n-1 )
234                {
235                    if( pivots[k]>=0 )
236                    {
237                       
238                        //
239                        // 1 x 1 diagonal block
240                        //
241                        // Multiply by inv(U'(K+1)), where U(K+1) is the transformation
242                        // stored in column K+1 of A.
243                        //
244                        v = 0.0;
245                        for(i_=0; i_<=k-1;i_++)
246                        {
247                            v += b[i_]*a[i_,k];
248                        }
249                        b[k] = b[k]-v;
250                       
251                        //
252                        // Interchange rows K+1 and IPIV(K+1).
253                        //
254                        kp = pivots[k];
255                        if( kp!=k )
256                        {
257                            v = b[k];
258                            b[k] = b[kp];
259                            b[kp] = v;
260                        }
261                        k = k+1;
262                    }
263                    else
264                    {
265                       
266                        //
267                        // 2 x 2 diagonal block
268                        //
269                        // Multiply by inv(U'(K+1+1)), where U(K+1+1) is the transformation
270                        // stored in columns K+1 and K+1+1 of A.
271                        //
272                        v = 0.0;
273                        for(i_=0; i_<=k-1;i_++)
274                        {
275                            v += b[i_]*a[i_,k];
276                        }
277                        b[k] = b[k]-v;
278                        v = 0.0;
279                        for(i_=0; i_<=k-1;i_++)
280                        {
281                            v += b[i_]*a[i_,k+1];
282                        }
283                        b[k+1] = b[k+1]-v;
284                       
285                        //
286                        // Interchange rows K+1 and -IPIV(K+1).
287                        //
288                        kp = pivots[k]+n;
289                        if( kp!=k )
290                        {
291                            v = b[k];
292                            b[k] = b[kp];
293                            b[kp] = v;
294                        }
295                        k = k+2;
296                    }
297                }
298            }
299            else
300            {
301               
302                //
303                // Solve A*X = B, where A = L*D*L'.
304                //
305                // First solve L*D*X = B, overwriting B with X.
306                //
307                // K+1 is the main loop index, increasing from 1 to N in steps of
308                // 1 or 2, depending on the size of the diagonal blocks.
309                //
310                k = 0;
311                while( k<=n-1 )
312                {
313                    if( pivots[k]>=0 )
314                    {
315                       
316                        //
317                        // 1 x 1 diagonal block
318                        //
319                        // Interchange rows K+1 and IPIV(K+1).
320                        //
321                        kp = pivots[k];
322                        if( kp!=k )
323                        {
324                            v = b[k];
325                            b[k] = b[kp];
326                            b[kp] = v;
327                        }
328                       
329                        //
330                        // Multiply by inv(L(K+1)), where L(K+1) is the transformation
331                        // stored in column K+1 of A.
332                        //
333                        if( k+1<n )
334                        {
335                            v = b[k];
336                            for(i_=k+1; i_<=n-1;i_++)
337                            {
338                                b[i_] = b[i_] - v*a[i_,k];
339                            }
340                        }
341                       
342                        //
343                        // Multiply by the inverse of the diagonal block.
344                        //
345                        b[k] = b[k]/a[k,k];
346                        k = k+1;
347                    }
348                    else
349                    {
350                       
351                        //
352                        // 2 x 2 diagonal block
353                        //
354                        // Interchange rows K+1+1 and -IPIV(K+1).
355                        //
356                        kp = pivots[k]+n;
357                        if( kp!=k+1 )
358                        {
359                            v = b[k+1];
360                            b[k+1] = b[kp];
361                            b[kp] = v;
362                        }
363                       
364                        //
365                        // Multiply by inv(L(K+1)), where L(K+1) is the transformation
366                        // stored in columns K+1 and K+1+1 of A.
367                        //
368                        if( k+1<n-1 )
369                        {
370                            v = b[k];
371                            for(i_=k+2; i_<=n-1;i_++)
372                            {
373                                b[i_] = b[i_] - v*a[i_,k];
374                            }
375                            v = b[k+1];
376                            for(i_=k+2; i_<=n-1;i_++)
377                            {
378                                b[i_] = b[i_] - v*a[i_,k+1];
379                            }
380                        }
381                       
382                        //
383                        // Multiply by the inverse of the diagonal block.
384                        //
385                        akm1k = a[k+1,k];
386                        akm1 = a[k,k]/akm1k;
387                        ak = a[k+1,k+1]/akm1k;
388                        denom = akm1*ak-1;
389                        bkm1 = b[k]/akm1k;
390                        bk = b[k+1]/akm1k;
391                        b[k] = (ak*bkm1-bk)/denom;
392                        b[k+1] = (akm1*bk-bkm1)/denom;
393                        k = k+2;
394                    }
395                }
396               
397                //
398                // Next solve L'*X = B, overwriting B with X.
399                //
400                // K+1 is the main loop index, decreasing from N to 1 in steps of
401                // 1 or 2, depending on the size of the diagonal blocks.
402                //
403                k = n-1;
404                while( k>=0 )
405                {
406                    if( pivots[k]>=0 )
407                    {
408                       
409                        //
410                        // 1 x 1 diagonal block
411                        //
412                        // Multiply by inv(L'(K+1)), where L(K+1) is the transformation
413                        // stored in column K+1 of A.
414                        //
415                        if( k+1<n )
416                        {
417                            v = 0.0;
418                            for(i_=k+1; i_<=n-1;i_++)
419                            {
420                                v += b[i_]*a[i_,k];
421                            }
422                            b[k] = b[k]-v;
423                        }
424                       
425                        //
426                        // Interchange rows K+1 and IPIV(K+1).
427                        //
428                        kp = pivots[k];
429                        if( kp!=k )
430                        {
431                            v = b[k];
432                            b[k] = b[kp];
433                            b[kp] = v;
434                        }
435                        k = k-1;
436                    }
437                    else
438                    {
439                       
440                        //
441                        // 2 x 2 diagonal block
442                        //
443                        // Multiply by inv(L'(K+1-1)), where L(K+1-1) is the transformation
444                        // stored in columns K+1-1 and K+1 of A.
445                        //
446                        if( k+1<n )
447                        {
448                            v = 0.0;
449                            for(i_=k+1; i_<=n-1;i_++)
450                            {
451                                v += b[i_]*a[i_,k];
452                            }
453                            b[k] = b[k]-v;
454                            v = 0.0;
455                            for(i_=k+1; i_<=n-1;i_++)
456                            {
457                                v += b[i_]*a[i_,k-1];
458                            }
459                            b[k-1] = b[k-1]-v;
460                        }
461                       
462                        //
463                        // Interchange rows K+1 and -IPIV(K+1).
464                        //
465                        kp = pivots[k]+n;
466                        if( kp!=k )
467                        {
468                            v = b[k];
469                            b[k] = b[kp];
470                            b[kp] = v;
471                        }
472                        k = k-2;
473                    }
474                }
475            }
476            x = new double[n-1+1];
477            for(i_=0; i_<=n-1;i_++)
478            {
479                x[i_] = b[i_];
480            }
481            return result;
482        }
483
484
485        /*************************************************************************
486        Solving a system of linear equations with a symmetric system matrix
487
488        Input parameters:
489            A       -   system matrix (upper or lower triangle).
490                        Array whose indexes range within [0..N-1, 0..N-1].
491            B       -   right side of a system.
492                        Array whose index ranges within [0..N-1].
493            N       -   size of matrix A.
494            IsUpper -   If IsUpper = True, A contains the upper triangle,
495                        otherwise A contains the lower triangle.
496
497        Output parameters:
498            X       -   solution of a system.
499                        Array whose index ranges within [0..N-1].
500
501        Result:
502            True, if the matrix is not singular. X contains the solution.
503            False, if the matrix is singular (the determinant of the matrix is equal
504        to 0). In this case, X doesn't contain a solution.
505
506          -- ALGLIB --
507             Copyright 2005 by Bochkanov Sergey
508        *************************************************************************/
509        public static bool smatrixsolve(double[,] a,
510            ref double[] b,
511            int n,
512            bool isupper,
513            ref double[] x)
514        {
515            bool result = new bool();
516            int[] pivots = new int[0];
517
518            a = (double[,])a.Clone();
519
520            ldlt.smatrixldlt(ref a, n, isupper, ref pivots);
521            result = smatrixldltsolve(ref a, ref pivots, b, n, isupper, ref x);
522            return result;
523        }
524
525
526        public static bool solvesystemldlt(ref double[,] a,
527            ref int[] pivots,
528            double[] b,
529            int n,
530            bool isupper,
531            ref double[] x)
532        {
533            bool result = new bool();
534            int i = 0;
535            int j = 0;
536            int k = 0;
537            int kp = 0;
538            int km1 = 0;
539            int km2 = 0;
540            int kp1 = 0;
541            int kp2 = 0;
542            double ak = 0;
543            double akm1 = 0;
544            double akm1k = 0;
545            double bk = 0;
546            double bkm1 = 0;
547            double denom = 0;
548            double v = 0;
549            int i_ = 0;
550
551            b = (double[])b.Clone();
552
553           
554            //
555            // Quick return if possible
556            //
557            result = true;
558            if( n==0 )
559            {
560                return result;
561            }
562           
563            //
564            // Check that the diagonal matrix D is nonsingular
565            //
566            if( isupper )
567            {
568               
569                //
570                // Upper triangular storage: examine D from bottom to top
571                //
572                for(i=n; i>=1; i--)
573                {
574                    if( pivots[i]>0 & (double)(a[i,i])==(double)(0) )
575                    {
576                        result = false;
577                        return result;
578                    }
579                }
580            }
581            else
582            {
583               
584                //
585                // Lower triangular storage: examine D from top to bottom.
586                //
587                for(i=1; i<=n; i++)
588                {
589                    if( pivots[i]>0 & (double)(a[i,i])==(double)(0) )
590                    {
591                        result = false;
592                        return result;
593                    }
594                }
595            }
596           
597            //
598            // Solve Ax = b
599            //
600            if( isupper )
601            {
602               
603                //
604                // Solve A*X = B, where A = U*D*U'.
605                //
606                // First solve U*D*X = B, overwriting B with X.
607                //
608                // K is the main loop index, decreasing from N to 1 in steps of
609                // 1 or 2, depending on the size of the diagonal blocks.
610                //
611                k = n;
612                while( k>=1 )
613                {
614                    if( pivots[k]>0 )
615                    {
616                       
617                        //
618                        // 1 x 1 diagonal block
619                        //
620                        // Interchange rows K and IPIV(K).
621                        //
622                        kp = pivots[k];
623                        if( kp!=k )
624                        {
625                            v = b[k];
626                            b[k] = b[kp];
627                            b[kp] = v;
628                        }
629                       
630                        //
631                        // Multiply by inv(U(K)), where U(K) is the transformation
632                        // stored in column K of A.
633                        //
634                        km1 = k-1;
635                        v = b[k];
636                        for(i_=1; i_<=km1;i_++)
637                        {
638                            b[i_] = b[i_] - v*a[i_,k];
639                        }
640                       
641                        //
642                        // Multiply by the inverse of the diagonal block.
643                        //
644                        b[k] = b[k]/a[k,k];
645                        k = k-1;
646                    }
647                    else
648                    {
649                       
650                        //
651                        // 2 x 2 diagonal block
652                        //
653                        // Interchange rows K-1 and -IPIV(K).
654                        //
655                        kp = -pivots[k];
656                        if( kp!=k-1 )
657                        {
658                            v = b[k-1];
659                            b[k-1] = b[kp];
660                            b[kp] = v;
661                        }
662                       
663                        //
664                        // Multiply by inv(U(K)), where U(K) is the transformation
665                        // stored in columns K-1 and K of A.
666                        //
667                        km2 = k-2;
668                        km1 = k-1;
669                        v = b[k];
670                        for(i_=1; i_<=km2;i_++)
671                        {
672                            b[i_] = b[i_] - v*a[i_,k];
673                        }
674                        v = b[k-1];
675                        for(i_=1; i_<=km2;i_++)
676                        {
677                            b[i_] = b[i_] - v*a[i_,km1];
678                        }
679                       
680                        //
681                        // Multiply by the inverse of the diagonal block.
682                        //
683                        akm1k = a[k-1,k];
684                        akm1 = a[k-1,k-1]/akm1k;
685                        ak = a[k,k]/akm1k;
686                        denom = akm1*ak-1;
687                        bkm1 = b[k-1]/akm1k;
688                        bk = b[k]/akm1k;
689                        b[k-1] = (ak*bkm1-bk)/denom;
690                        b[k] = (akm1*bk-bkm1)/denom;
691                        k = k-2;
692                    }
693                }
694               
695                //
696                // Next solve U'*X = B, overwriting B with X.
697                //
698                // K is the main loop index, increasing from 1 to N in steps of
699                // 1 or 2, depending on the size of the diagonal blocks.
700                //
701                k = 1;
702                while( k<=n )
703                {
704                    if( pivots[k]>0 )
705                    {
706                       
707                        //
708                        // 1 x 1 diagonal block
709                        //
710                        // Multiply by inv(U'(K)), where U(K) is the transformation
711                        // stored in column K of A.
712                        //
713                        km1 = k-1;
714                        v = 0.0;
715                        for(i_=1; i_<=km1;i_++)
716                        {
717                            v += b[i_]*a[i_,k];
718                        }
719                        b[k] = b[k]-v;
720                       
721                        //
722                        // Interchange rows K and IPIV(K).
723                        //
724                        kp = pivots[k];
725                        if( kp!=k )
726                        {
727                            v = b[k];
728                            b[k] = b[kp];
729                            b[kp] = v;
730                        }
731                        k = k+1;
732                    }
733                    else
734                    {
735                       
736                        //
737                        // 2 x 2 diagonal block
738                        //
739                        // Multiply by inv(U'(K+1)), where U(K+1) is the transformation
740                        // stored in columns K and K+1 of A.
741                        //
742                        km1 = k-1;
743                        kp1 = k+1;
744                        v = 0.0;
745                        for(i_=1; i_<=km1;i_++)
746                        {
747                            v += b[i_]*a[i_,k];
748                        }
749                        b[k] = b[k]-v;
750                        v = 0.0;
751                        for(i_=1; i_<=km1;i_++)
752                        {
753                            v += b[i_]*a[i_,kp1];
754                        }
755                        b[k+1] = b[k+1]-v;
756                       
757                        //
758                        // Interchange rows K and -IPIV(K).
759                        //
760                        kp = -pivots[k];
761                        if( kp!=k )
762                        {
763                            v = b[k];
764                            b[k] = b[kp];
765                            b[kp] = v;
766                        }
767                        k = k+2;
768                    }
769                }
770            }
771            else
772            {
773               
774                //
775                // Solve A*X = B, where A = L*D*L'.
776                //
777                // First solve L*D*X = B, overwriting B with X.
778                //
779                // K is the main loop index, increasing from 1 to N in steps of
780                // 1 or 2, depending on the size of the diagonal blocks.
781                //
782                k = 1;
783                while( k<=n )
784                {
785                    if( pivots[k]>0 )
786                    {
787                       
788                        //
789                        // 1 x 1 diagonal block
790                        //
791                        // Interchange rows K and IPIV(K).
792                        //
793                        kp = pivots[k];
794                        if( kp!=k )
795                        {
796                            v = b[k];
797                            b[k] = b[kp];
798                            b[kp] = v;
799                        }
800                       
801                        //
802                        // Multiply by inv(L(K)), where L(K) is the transformation
803                        // stored in column K of A.
804                        //
805                        if( k<n )
806                        {
807                            kp1 = k+1;
808                            v = b[k];
809                            for(i_=kp1; i_<=n;i_++)
810                            {
811                                b[i_] = b[i_] - v*a[i_,k];
812                            }
813                        }
814                       
815                        //
816                        // Multiply by the inverse of the diagonal block.
817                        //
818                        b[k] = b[k]/a[k,k];
819                        k = k+1;
820                    }
821                    else
822                    {
823                       
824                        //
825                        // 2 x 2 diagonal block
826                        //
827                        // Interchange rows K+1 and -IPIV(K).
828                        //
829                        kp = -pivots[k];
830                        if( kp!=k+1 )
831                        {
832                            v = b[k+1];
833                            b[k+1] = b[kp];
834                            b[kp] = v;
835                        }
836                       
837                        //
838                        // Multiply by inv(L(K)), where L(K) is the transformation
839                        // stored in columns K and K+1 of A.
840                        //
841                        if( k<n-1 )
842                        {
843                            kp1 = k+1;
844                            kp2 = k+2;
845                            v = b[k];
846                            for(i_=kp2; i_<=n;i_++)
847                            {
848                                b[i_] = b[i_] - v*a[i_,k];
849                            }
850                            v = b[k+1];
851                            for(i_=kp2; i_<=n;i_++)
852                            {
853                                b[i_] = b[i_] - v*a[i_,kp1];
854                            }
855                        }
856                       
857                        //
858                        // Multiply by the inverse of the diagonal block.
859                        //
860                        akm1k = a[k+1,k];
861                        akm1 = a[k,k]/akm1k;
862                        ak = a[k+1,k+1]/akm1k;
863                        denom = akm1*ak-1;
864                        bkm1 = b[k]/akm1k;
865                        bk = b[k+1]/akm1k;
866                        b[k] = (ak*bkm1-bk)/denom;
867                        b[k+1] = (akm1*bk-bkm1)/denom;
868                        k = k+2;
869                    }
870                }
871               
872                //
873                // Next solve L'*X = B, overwriting B with X.
874                //
875                // K is the main loop index, decreasing from N to 1 in steps of
876                // 1 or 2, depending on the size of the diagonal blocks.
877                //
878                k = n;
879                while( k>=1 )
880                {
881                    if( pivots[k]>0 )
882                    {
883                       
884                        //
885                        // 1 x 1 diagonal block
886                        //
887                        // Multiply by inv(L'(K)), where L(K) is the transformation
888                        // stored in column K of A.
889                        //
890                        if( k<n )
891                        {
892                            kp1 = k+1;
893                            v = 0.0;
894                            for(i_=kp1; i_<=n;i_++)
895                            {
896                                v += b[i_]*a[i_,k];
897                            }
898                            b[k] = b[k]-v;
899                        }
900                       
901                        //
902                        // Interchange rows K and IPIV(K).
903                        //
904                        kp = pivots[k];
905                        if( kp!=k )
906                        {
907                            v = b[k];
908                            b[k] = b[kp];
909                            b[kp] = v;
910                        }
911                        k = k-1;
912                    }
913                    else
914                    {
915                       
916                        //
917                        // 2 x 2 diagonal block
918                        //
919                        // Multiply by inv(L'(K-1)), where L(K-1) is the transformation
920                        // stored in columns K-1 and K of A.
921                        //
922                        if( k<n )
923                        {
924                            kp1 = k+1;
925                            km1 = k-1;
926                            v = 0.0;
927                            for(i_=kp1; i_<=n;i_++)
928                            {
929                                v += b[i_]*a[i_,k];
930                            }
931                            b[k] = b[k]-v;
932                            v = 0.0;
933                            for(i_=kp1; i_<=n;i_++)
934                            {
935                                v += b[i_]*a[i_,km1];
936                            }
937                            b[k-1] = b[k-1]-v;
938                        }
939                       
940                        //
941                        // Interchange rows K and -IPIV(K).
942                        //
943                        kp = -pivots[k];
944                        if( kp!=k )
945                        {
946                            v = b[k];
947                            b[k] = b[kp];
948                            b[kp] = v;
949                        }
950                        k = k-2;
951                    }
952                }
953            }
954            x = new double[n+1];
955            for(i_=1; i_<=n;i_++)
956            {
957                x[i_] = b[i_];
958            }
959            return result;
960        }
961
962
963        public static bool solvesymmetricsystem(double[,] a,
964            double[] b,
965            int n,
966            bool isupper,
967            ref double[] x)
968        {
969            bool result = new bool();
970            int[] pivots = new int[0];
971
972            a = (double[,])a.Clone();
973            b = (double[])b.Clone();
974
975            ldlt.ldltdecomposition(ref a, n, isupper, ref pivots);
976            result = solvesystemldlt(ref a, ref pivots, b, n, isupper, ref x);
977            return result;
978        }
979    }
980}
Note: See TracBrowser for help on using the repository browser.