Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/ALGLIB/hcholesky.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: 13.3 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 hcholesky
32    {
33        /*************************************************************************
34        Cholesky decomposition
35
36        The algorithm computes Cholesky decomposition  of  a  Hermitian  positive-
37        definite matrix.
38
39        The result of an algorithm is a representation of matrix A as A = U'*U  or
40        A = L*L' (here X' detones conj(X^T)).
41
42        Input parameters:
43            A       -   upper or lower triangle of a factorized matrix.
44                        array with elements [0..N-1, 0..N-1].
45            N       -   size of matrix A.
46            IsUpper -   if IsUpper=True, then A contains an upper triangle of
47                        a symmetric matrix, otherwise A contains a lower one.
48
49        Output parameters:
50            A       -   the result of factorization. If IsUpper=True, then
51                        the upper triangle contains matrix U, so that A = U'*U,
52                        and the elements below the main diagonal are not modified.
53                        Similarly, if IsUpper = False.
54
55        Result:
56            If the matrix is positive-definite, the function returns True.
57            Otherwise, the function returns False. This means that the
58            factorization could not be carried out.
59
60          -- LAPACK routine (version 3.0) --
61             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
62             Courant Institute, Argonne National Lab, and Rice University
63             February 29, 1992
64        *************************************************************************/
65        public static bool hmatrixcholesky(ref AP.Complex[,] a,
66            int n,
67            bool isupper)
68        {
69            bool result = new bool();
70            int j = 0;
71            double ajj = 0;
72            AP.Complex v = 0;
73            double r = 0;
74            AP.Complex[] t = new AP.Complex[0];
75            AP.Complex[] t2 = new AP.Complex[0];
76            AP.Complex[] t3 = new AP.Complex[0];
77            int i = 0;
78            AP.Complex[,] a1 = new AP.Complex[0,0];
79            int i_ = 0;
80
81            if( !isupper )
82            {
83                a1 = new AP.Complex[n+1, n+1];
84                for(i=1; i<=n; i++)
85                {
86                    for(j=1; j<=n; j++)
87                    {
88                        a1[i,j] = a[i-1,j-1];
89                    }
90                }
91                result = hermitiancholeskydecomposition(ref a1, n, isupper);
92                for(i=1; i<=n; i++)
93                {
94                    for(j=1; j<=n; j++)
95                    {
96                        a[i-1,j-1] = a1[i,j];
97                    }
98                }
99                return result;
100            }
101            t = new AP.Complex[n-1+1];
102            t2 = new AP.Complex[n-1+1];
103            t3 = new AP.Complex[n-1+1];
104            result = true;
105            if( n<0 )
106            {
107                result = false;
108                return result;
109            }
110           
111            //
112            // Quick return if possible
113            //
114            if( n==0 )
115            {
116                return result;
117            }
118            if( isupper )
119            {
120               
121                //
122                // Compute the Cholesky factorization A = U'*U.
123                //
124                for(j=0; j<=n-1; j++)
125                {
126                   
127                    //
128                    // Compute U(J,J) and test for non-positive-definiteness.
129                    //
130                    v = 0.0;
131                    for(i_=0; i_<=j-1;i_++)
132                    {
133                        v += AP.Math.Conj(a[i_,j])*a[i_,j];
134                    }
135                    ajj = (a[j,j]-v).x;
136                    if( (double)(ajj)<=(double)(0) )
137                    {
138                        a[j,j] = ajj;
139                        result = false;
140                        return result;
141                    }
142                    ajj = Math.Sqrt(ajj);
143                    a[j,j] = ajj;
144                   
145                    //
146                    // Compute elements J+1:N-1 of row J.
147                    //
148                    if( j<n-1 )
149                    {
150                        for(i_=0; i_<=j-1;i_++)
151                        {
152                            t2[i_] = AP.Math.Conj(a[i_,j]);
153                        }
154                        for(i_=j+1; i_<=n-1;i_++)
155                        {
156                            t3[i_] = a[j,i_];
157                        }
158                        cblas.complexmatrixvectormultiply(ref a, 0, j-1, j+1, n-1, true, false, ref t2, 0, j-1, -1.0, ref t3, j+1, n-1, 1.0, ref t);
159                        for(i_=j+1; i_<=n-1;i_++)
160                        {
161                            a[j,i_] = t3[i_];
162                        }
163                        r = 1/ajj;
164                        for(i_=j+1; i_<=n-1;i_++)
165                        {
166                            a[j,i_] = r*a[j,i_];
167                        }
168                    }
169                }
170            }
171            else
172            {
173               
174                //
175                // Compute the Cholesky factorization A = L*L'.
176                //
177                for(j=0; j<=n-1; j++)
178                {
179                   
180                    //
181                    // Compute L(J+1,J+1) and test for non-positive-definiteness.
182                    //
183                    v = 0.0;
184                    for(i_=0; i_<=j-1;i_++)
185                    {
186                        v += AP.Math.Conj(a[j,i_])*a[j,i_];
187                    }
188                    ajj = (a[j,j]-v).x;
189                    if( (double)(ajj)<=(double)(0) )
190                    {
191                        a[j,j] = ajj;
192                        result = false;
193                        return result;
194                    }
195                    ajj = Math.Sqrt(ajj);
196                    a[j,j] = ajj;
197                   
198                    //
199                    // Compute elements J+1:N of column J.
200                    //
201                    if( j<n-1 )
202                    {
203                        for(i_=0; i_<=j-1;i_++)
204                        {
205                            t2[i_] = AP.Math.Conj(a[j,i_]);
206                        }
207                        for(i_=j+1; i_<=n-1;i_++)
208                        {
209                            t3[i_] = a[i_,j];
210                        }
211                        cblas.complexmatrixvectormultiply(ref a, j+1, n-1, 0, j-1, false, false, ref t2, 0, j-1, -1.0, ref t3, j+1, n-1, 1.0, ref t);
212                        for(i_=j+1; i_<=n-1;i_++)
213                        {
214                            a[i_,j] = t3[i_];
215                        }
216                        r = 1/ajj;
217                        for(i_=j+1; i_<=n-1;i_++)
218                        {
219                            a[i_,j] = r*a[i_,j];
220                        }
221                    }
222                }
223            }
224            return result;
225        }
226
227
228        public static bool hermitiancholeskydecomposition(ref AP.Complex[,] a,
229            int n,
230            bool isupper)
231        {
232            bool result = new bool();
233            int j = 0;
234            double ajj = 0;
235            AP.Complex v = 0;
236            double r = 0;
237            AP.Complex[] t = new AP.Complex[0];
238            AP.Complex[] t2 = new AP.Complex[0];
239            AP.Complex[] t3 = new AP.Complex[0];
240            int i_ = 0;
241
242            t = new AP.Complex[n+1];
243            t2 = new AP.Complex[n+1];
244            t3 = new AP.Complex[n+1];
245            result = true;
246            if( n<0 )
247            {
248                result = false;
249                return result;
250            }
251           
252            //
253            // Quick return if possible
254            //
255            if( n==0 )
256            {
257                return result;
258            }
259            if( isupper )
260            {
261               
262                //
263                // Compute the Cholesky factorization A = U'*U.
264                //
265                for(j=1; j<=n; j++)
266                {
267                   
268                    //
269                    // Compute U(J,J) and test for non-positive-definiteness.
270                    //
271                    v = 0.0;
272                    for(i_=1; i_<=j-1;i_++)
273                    {
274                        v += AP.Math.Conj(a[i_,j])*a[i_,j];
275                    }
276                    ajj = (a[j,j]-v).x;
277                    if( (double)(ajj)<=(double)(0) )
278                    {
279                        a[j,j] = ajj;
280                        result = false;
281                        return result;
282                    }
283                    ajj = Math.Sqrt(ajj);
284                    a[j,j] = ajj;
285                   
286                    //
287                    // Compute elements J+1:N of row J.
288                    //
289                    if( j<n )
290                    {
291                        for(i_=1; i_<=j-1;i_++)
292                        {
293                            a[i_,j] = AP.Math.Conj(a[i_,j]);
294                        }
295                        for(i_=1; i_<=j-1;i_++)
296                        {
297                            t2[i_] = a[i_,j];
298                        }
299                        for(i_=j+1; i_<=n;i_++)
300                        {
301                            t3[i_] = a[j,i_];
302                        }
303                        cblas.complexmatrixvectormultiply(ref a, 1, j-1, j+1, n, true, false, ref t2, 1, j-1, -1.0, ref t3, j+1, n, 1.0, ref t);
304                        for(i_=j+1; i_<=n;i_++)
305                        {
306                            a[j,i_] = t3[i_];
307                        }
308                        for(i_=1; i_<=j-1;i_++)
309                        {
310                            a[i_,j] = AP.Math.Conj(a[i_,j]);
311                        }
312                        r = 1/ajj;
313                        for(i_=j+1; i_<=n;i_++)
314                        {
315                            a[j,i_] = r*a[j,i_];
316                        }
317                    }
318                }
319            }
320            else
321            {
322               
323                //
324                // Compute the Cholesky factorization A = L*L'.
325                //
326                for(j=1; j<=n; j++)
327                {
328                   
329                    //
330                    // Compute L(J,J) and test for non-positive-definiteness.
331                    //
332                    v = 0.0;
333                    for(i_=1; i_<=j-1;i_++)
334                    {
335                        v += AP.Math.Conj(a[j,i_])*a[j,i_];
336                    }
337                    ajj = (a[j,j]-v).x;
338                    if( (double)(ajj)<=(double)(0) )
339                    {
340                        a[j,j] = ajj;
341                        result = false;
342                        return result;
343                    }
344                    ajj = Math.Sqrt(ajj);
345                    a[j,j] = ajj;
346                   
347                    //
348                    // Compute elements J+1:N of column J.
349                    //
350                    if( j<n )
351                    {
352                        for(i_=1; i_<=j-1;i_++)
353                        {
354                            a[j,i_] = AP.Math.Conj(a[j,i_]);
355                        }
356                        for(i_=1; i_<=j-1;i_++)
357                        {
358                            t2[i_] = a[j,i_];
359                        }
360                        for(i_=j+1; i_<=n;i_++)
361                        {
362                            t3[i_] = a[i_,j];
363                        }
364                        cblas.complexmatrixvectormultiply(ref a, j+1, n, 1, j-1, false, false, ref t2, 1, j-1, -1.0, ref t3, j+1, n, 1.0, ref t);
365                        for(i_=j+1; i_<=n;i_++)
366                        {
367                            a[i_,j] = t3[i_];
368                        }
369                        for(i_=1; i_<=j-1;i_++)
370                        {
371                            a[j,i_] = AP.Math.Conj(a[j,i_]);
372                        }
373                        r = 1/ajj;
374                        for(i_=j+1; i_<=n;i_++)
375                        {
376                            a[i_,j] = r*a[i_,j];
377                        }
378                    }
379                }
380            }
381            return result;
382        }
383    }
384}
Note: See TracBrowser for help on using the repository browser.