Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/trinverse.cs @ 2645

Last change on this file since 2645 was 2645, checked in by mkommend, 14 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 12.5 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 trinverse
32    {
33        /*************************************************************************
34        Triangular matrix inversion
35
36        The subroutine inverts the following types of matrices:
37            * upper triangular
38            * upper triangular with unit diagonal
39            * lower triangular
40            * lower triangular with unit diagonal
41
42        In case of an upper (lower) triangular matrix,  the  inverse  matrix  will
43        also be upper (lower) triangular, and after the end of the algorithm,  the
44        inverse matrix replaces the source matrix. The elements  below (above) the
45        main diagonal are not changed by the algorithm.
46
47        If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
48        diagonal, and the diagonal elements are not passed to the algorithm.
49
50        Input parameters:
51            A       -   matrix.
52                        Array whose indexes range within [0..N-1, 0..N-1].
53            N       -   size of matrix A.
54            IsUpper -   True, if the matrix is upper triangular.
55            IsunitTriangular
56                    -   True, if the matrix has a unit diagonal.
57
58        Output parameters:
59            A       -   inverse matrix (if the problem is not degenerate).
60
61        Result:
62            True, if the matrix is not singular.
63            False, if the matrix is singular.
64
65          -- LAPACK routine (version 3.0) --
66             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
67             Courant Institute, Argonne National Lab, and Rice University
68             February 29, 1992
69        *************************************************************************/
70        public static bool rmatrixtrinverse(ref double[,] a,
71            int n,
72            bool isupper,
73            bool isunittriangular)
74        {
75            bool result = new bool();
76            bool nounit = new bool();
77            int i = 0;
78            int j = 0;
79            double v = 0;
80            double ajj = 0;
81            double[] t = new double[0];
82            int i_ = 0;
83
84            result = true;
85            t = new double[n-1+1];
86           
87            //
88            // Test the input parameters.
89            //
90            nounit = !isunittriangular;
91            if( isupper )
92            {
93               
94                //
95                // Compute inverse of upper triangular matrix.
96                //
97                for(j=0; j<=n-1; j++)
98                {
99                    if( nounit )
100                    {
101                        if( (double)(a[j,j])==(double)(0) )
102                        {
103                            result = false;
104                            return result;
105                        }
106                        a[j,j] = 1/a[j,j];
107                        ajj = -a[j,j];
108                    }
109                    else
110                    {
111                        ajj = -1;
112                    }
113                   
114                    //
115                    // Compute elements 1:j-1 of j-th column.
116                    //
117                    if( j>0 )
118                    {
119                        for(i_=0; i_<=j-1;i_++)
120                        {
121                            t[i_] = a[i_,j];
122                        }
123                        for(i=0; i<=j-1; i++)
124                        {
125                            if( i<j-1 )
126                            {
127                                v = 0.0;
128                                for(i_=i+1; i_<=j-1;i_++)
129                                {
130                                    v += a[i,i_]*t[i_];
131                                }
132                            }
133                            else
134                            {
135                                v = 0;
136                            }
137                            if( nounit )
138                            {
139                                a[i,j] = v+a[i,i]*t[i];
140                            }
141                            else
142                            {
143                                a[i,j] = v+t[i];
144                            }
145                        }
146                        for(i_=0; i_<=j-1;i_++)
147                        {
148                            a[i_,j] = ajj*a[i_,j];
149                        }
150                    }
151                }
152            }
153            else
154            {
155               
156                //
157                // Compute inverse of lower triangular matrix.
158                //
159                for(j=n-1; j>=0; j--)
160                {
161                    if( nounit )
162                    {
163                        if( (double)(a[j,j])==(double)(0) )
164                        {
165                            result = false;
166                            return result;
167                        }
168                        a[j,j] = 1/a[j,j];
169                        ajj = -a[j,j];
170                    }
171                    else
172                    {
173                        ajj = -1;
174                    }
175                    if( j<n-1 )
176                    {
177                       
178                        //
179                        // Compute elements j+1:n of j-th column.
180                        //
181                        for(i_=j+1; i_<=n-1;i_++)
182                        {
183                            t[i_] = a[i_,j];
184                        }
185                        for(i=j+1; i<=n-1; i++)
186                        {
187                            if( i>j+1 )
188                            {
189                                v = 0.0;
190                                for(i_=j+1; i_<=i-1;i_++)
191                                {
192                                    v += a[i,i_]*t[i_];
193                                }
194                            }
195                            else
196                            {
197                                v = 0;
198                            }
199                            if( nounit )
200                            {
201                                a[i,j] = v+a[i,i]*t[i];
202                            }
203                            else
204                            {
205                                a[i,j] = v+t[i];
206                            }
207                        }
208                        for(i_=j+1; i_<=n-1;i_++)
209                        {
210                            a[i_,j] = ajj*a[i_,j];
211                        }
212                    }
213                }
214            }
215            return result;
216        }
217
218
219        public static bool invtriangular(ref double[,] a,
220            int n,
221            bool isupper,
222            bool isunittriangular)
223        {
224            bool result = new bool();
225            bool nounit = new bool();
226            int i = 0;
227            int j = 0;
228            int nmj = 0;
229            int jm1 = 0;
230            int jp1 = 0;
231            double v = 0;
232            double ajj = 0;
233            double[] t = new double[0];
234            int i_ = 0;
235
236            result = true;
237            t = new double[n+1];
238           
239            //
240            // Test the input parameters.
241            //
242            nounit = !isunittriangular;
243            if( isupper )
244            {
245               
246                //
247                // Compute inverse of upper triangular matrix.
248                //
249                for(j=1; j<=n; j++)
250                {
251                    if( nounit )
252                    {
253                        if( (double)(a[j,j])==(double)(0) )
254                        {
255                            result = false;
256                            return result;
257                        }
258                        a[j,j] = 1/a[j,j];
259                        ajj = -a[j,j];
260                    }
261                    else
262                    {
263                        ajj = -1;
264                    }
265                   
266                    //
267                    // Compute elements 1:j-1 of j-th column.
268                    //
269                    if( j>1 )
270                    {
271                        jm1 = j-1;
272                        for(i_=1; i_<=jm1;i_++)
273                        {
274                            t[i_] = a[i_,j];
275                        }
276                        for(i=1; i<=j-1; i++)
277                        {
278                            if( i<j-1 )
279                            {
280                                v = 0.0;
281                                for(i_=i+1; i_<=jm1;i_++)
282                                {
283                                    v += a[i,i_]*t[i_];
284                                }
285                            }
286                            else
287                            {
288                                v = 0;
289                            }
290                            if( nounit )
291                            {
292                                a[i,j] = v+a[i,i]*t[i];
293                            }
294                            else
295                            {
296                                a[i,j] = v+t[i];
297                            }
298                        }
299                        for(i_=1; i_<=jm1;i_++)
300                        {
301                            a[i_,j] = ajj*a[i_,j];
302                        }
303                    }
304                }
305            }
306            else
307            {
308               
309                //
310                // Compute inverse of lower triangular matrix.
311                //
312                for(j=n; j>=1; j--)
313                {
314                    if( nounit )
315                    {
316                        if( (double)(a[j,j])==(double)(0) )
317                        {
318                            result = false;
319                            return result;
320                        }
321                        a[j,j] = 1/a[j,j];
322                        ajj = -a[j,j];
323                    }
324                    else
325                    {
326                        ajj = -1;
327                    }
328                    if( j<n )
329                    {
330                       
331                        //
332                        // Compute elements j+1:n of j-th column.
333                        //
334                        nmj = n-j;
335                        jp1 = j+1;
336                        for(i_=jp1; i_<=n;i_++)
337                        {
338                            t[i_] = a[i_,j];
339                        }
340                        for(i=j+1; i<=n; i++)
341                        {
342                            if( i>j+1 )
343                            {
344                                v = 0.0;
345                                for(i_=jp1; i_<=i-1;i_++)
346                                {
347                                    v += a[i,i_]*t[i_];
348                                }
349                            }
350                            else
351                            {
352                                v = 0;
353                            }
354                            if( nounit )
355                            {
356                                a[i,j] = v+a[i,i]*t[i];
357                            }
358                            else
359                            {
360                                a[i,j] = v+t[i];
361                            }
362                        }
363                        for(i_=jp1; i_<=n;i_++)
364                        {
365                            a[i_,j] = ajj*a[i_,j];
366                        }
367                    }
368                }
369            }
370            return result;
371        }
372    }
373}
Note: See TracBrowser for help on using the repository browser.