Free cookie consent management tool by TermsFeed Policy Generator

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