Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/trfac.cs @ 3021

Last change on this file since 3021 was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 60.2 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 trfac
32    {
33        /*************************************************************************
34        LU decomposition of a general real matrix with row pivoting
35
36        A is represented as A = P*L*U, where:
37        * L is lower unitriangular matrix
38        * U is upper triangular matrix
39        * P = P0*P1*...*PK, K=min(M,N)-1,
40          Pi - permutation matrix for I and Pivots[I]
41
42        This is cache-oblivous implementation of LU decomposition.
43        It is optimized for square matrices. As for rectangular matrices:
44        * best case - M>>N
45        * worst case - N>>M, small M, large N, matrix does not fit in CPU cache
46
47        INPUT PARAMETERS:
48            A       -   array[0..M-1, 0..N-1].
49            M       -   number of rows in matrix A.
50            N       -   number of columns in matrix A.
51
52
53        OUTPUT PARAMETERS:
54            A       -   matrices L and U in compact form:
55                        * L is stored under main diagonal
56                        * U is stored on and above main diagonal
57            Pivots  -   permutation matrix in compact form.
58                        array[0..Min(M-1,N-1)].
59
60          -- ALGLIB routine --
61             10.01.2010
62             Bochkanov Sergey
63        *************************************************************************/
64        public static void rmatrixlu(ref double[,] a,
65            int m,
66            int n,
67            ref int[] pivots)
68        {
69            System.Diagnostics.Debug.Assert(m>0, "RMatrixLU: incorrect M!");
70            System.Diagnostics.Debug.Assert(n>0, "RMatrixLU: incorrect N!");
71            rmatrixplu(ref a, m, n, ref pivots);
72        }
73
74
75        /*************************************************************************
76        LU decomposition of a general complex matrix with row pivoting
77
78        A is represented as A = P*L*U, where:
79        * L is lower unitriangular matrix
80        * U is upper triangular matrix
81        * P = P0*P1*...*PK, K=min(M,N)-1,
82          Pi - permutation matrix for I and Pivots[I]
83
84        This is cache-oblivous implementation of LU decomposition. It is optimized
85        for square matrices. As for rectangular matrices:
86        * best case - M>>N
87        * worst case - N>>M, small M, large N, matrix does not fit in CPU cache
88
89        INPUT PARAMETERS:
90            A       -   array[0..M-1, 0..N-1].
91            M       -   number of rows in matrix A.
92            N       -   number of columns in matrix A.
93
94
95        OUTPUT PARAMETERS:
96            A       -   matrices L and U in compact form:
97                        * L is stored under main diagonal
98                        * U is stored on and above main diagonal
99            Pivots  -   permutation matrix in compact form.
100                        array[0..Min(M-1,N-1)].
101
102          -- ALGLIB routine --
103             10.01.2010
104             Bochkanov Sergey
105        *************************************************************************/
106        public static void cmatrixlu(ref AP.Complex[,] a,
107            int m,
108            int n,
109            ref int[] pivots)
110        {
111            System.Diagnostics.Debug.Assert(m>0, "CMatrixLU: incorrect M!");
112            System.Diagnostics.Debug.Assert(n>0, "CMatrixLU: incorrect N!");
113            cmatrixplu(ref a, m, n, ref pivots);
114        }
115
116
117        /*************************************************************************
118        Cache-oblivious Cholesky decomposition
119
120        The algorithm computes Cholesky decomposition  of  a  Hermitian  positive-
121        definite matrix. The result of an algorithm is a representation  of  A  as
122        A=U'*U  or A=L*L' (here X' detones conj(X^T)).
123
124        INPUT PARAMETERS:
125            A       -   upper or lower triangle of a factorized matrix.
126                        array with elements [0..N-1, 0..N-1].
127            N       -   size of matrix A.
128            IsUpper -   if IsUpper=True, then A contains an upper triangle of
129                        a symmetric matrix, otherwise A contains a lower one.
130
131        OUTPUT PARAMETERS:
132            A       -   the result of factorization. If IsUpper=True, then
133                        the upper triangle contains matrix U, so that A = U'*U,
134                        and the elements below the main diagonal are not modified.
135                        Similarly, if IsUpper = False.
136
137        RESULT:
138            If  the  matrix  is  positive-definite,  the  function  returns  True.
139            Otherwise, the function returns False. Contents of A is not determined
140            in such case.
141
142          -- ALGLIB routine --
143             15.12.2009
144             Bochkanov Sergey
145        *************************************************************************/
146        public static bool hpdmatrixcholesky(ref AP.Complex[,] a,
147            int n,
148            bool isupper)
149        {
150            bool result = new bool();
151            AP.Complex[] tmp = new AP.Complex[0];
152
153            if( n<1 )
154            {
155                result = false;
156                return result;
157            }
158            tmp = new AP.Complex[2*n];
159            result = hpdmatrixcholeskyrec(ref a, 0, n, isupper, ref tmp);
160            return result;
161        }
162
163
164        /*************************************************************************
165        Cache-oblivious Cholesky decomposition
166
167        The algorithm computes Cholesky decomposition  of  a  symmetric  positive-
168        definite matrix. The result of an algorithm is a representation  of  A  as
169        A=U^T*U  or A=L*L^T
170
171        INPUT PARAMETERS:
172            A       -   upper or lower triangle of a factorized matrix.
173                        array with elements [0..N-1, 0..N-1].
174            N       -   size of matrix A.
175            IsUpper -   if IsUpper=True, then A contains an upper triangle of
176                        a symmetric matrix, otherwise A contains a lower one.
177
178        OUTPUT PARAMETERS:
179            A       -   the result of factorization. If IsUpper=True, then
180                        the upper triangle contains matrix U, so that A = U^T*U,
181                        and the elements below the main diagonal are not modified.
182                        Similarly, if IsUpper = False.
183
184        RESULT:
185            If  the  matrix  is  positive-definite,  the  function  returns  True.
186            Otherwise, the function returns False. Contents of A is not determined
187            in such case.
188
189          -- ALGLIB routine --
190             15.12.2009
191             Bochkanov Sergey
192        *************************************************************************/
193        public static bool spdmatrixcholesky(ref double[,] a,
194            int n,
195            bool isupper)
196        {
197            bool result = new bool();
198            double[] tmp = new double[0];
199
200            if( n<1 )
201            {
202                result = false;
203                return result;
204            }
205            tmp = new double[2*n];
206            result = spdmatrixcholeskyrec(ref a, 0, n, isupper, ref tmp);
207            return result;
208        }
209
210
211        public static void rmatrixlup(ref double[,] a,
212            int m,
213            int n,
214            ref int[] pivots)
215        {
216            double[] tmp = new double[0];
217            int i = 0;
218            int j = 0;
219            double mx = 0;
220            double v = 0;
221            int i_ = 0;
222
223           
224            //
225            // Internal LU decomposition subroutine.
226            // Never call it directly.
227            //
228            System.Diagnostics.Debug.Assert(m>0, "RMatrixLUP: incorrect M!");
229            System.Diagnostics.Debug.Assert(n>0, "RMatrixLUP: incorrect N!");
230           
231            //
232            // Scale matrix to avoid overflows,
233            // decompose it, then scale back.
234            //
235            mx = 0;
236            for(i=0; i<=m-1; i++)
237            {
238                for(j=0; j<=n-1; j++)
239                {
240                    mx = Math.Max(mx, Math.Abs(a[i,j]));
241                }
242            }
243            if( (double)(mx)!=(double)(0) )
244            {
245                v = 1/mx;
246                for(i=0; i<=m-1; i++)
247                {
248                    for(i_=0; i_<=n-1;i_++)
249                    {
250                        a[i,i_] = v*a[i,i_];
251                    }
252                }
253            }
254            pivots = new int[Math.Min(m, n)];
255            tmp = new double[2*Math.Max(m, n)];
256            rmatrixluprec(ref a, 0, m, n, ref pivots, ref tmp);
257            if( (double)(mx)!=(double)(0) )
258            {
259                v = mx;
260                for(i=0; i<=m-1; i++)
261                {
262                    for(i_=0; i_<=Math.Min(i, n-1);i_++)
263                    {
264                        a[i,i_] = v*a[i,i_];
265                    }
266                }
267            }
268        }
269
270
271        public static void cmatrixlup(ref AP.Complex[,] a,
272            int m,
273            int n,
274            ref int[] pivots)
275        {
276            AP.Complex[] tmp = new AP.Complex[0];
277            int i = 0;
278            int j = 0;
279            double mx = 0;
280            double v = 0;
281            int i_ = 0;
282
283           
284            //
285            // Internal LU decomposition subroutine.
286            // Never call it directly.
287            //
288            System.Diagnostics.Debug.Assert(m>0, "CMatrixLUP: incorrect M!");
289            System.Diagnostics.Debug.Assert(n>0, "CMatrixLUP: incorrect N!");
290           
291            //
292            // Scale matrix to avoid overflows,
293            // decompose it, then scale back.
294            //
295            mx = 0;
296            for(i=0; i<=m-1; i++)
297            {
298                for(j=0; j<=n-1; j++)
299                {
300                    mx = Math.Max(mx, AP.Math.AbsComplex(a[i,j]));
301                }
302            }
303            if( (double)(mx)!=(double)(0) )
304            {
305                v = 1/mx;
306                for(i=0; i<=m-1; i++)
307                {
308                    for(i_=0; i_<=n-1;i_++)
309                    {
310                        a[i,i_] = v*a[i,i_];
311                    }
312                }
313            }
314            pivots = new int[Math.Min(m, n)];
315            tmp = new AP.Complex[2*Math.Max(m, n)];
316            cmatrixluprec(ref a, 0, m, n, ref pivots, ref tmp);
317            if( (double)(mx)!=(double)(0) )
318            {
319                v = mx;
320                for(i=0; i<=m-1; i++)
321                {
322                    for(i_=0; i_<=Math.Min(i, n-1);i_++)
323                    {
324                        a[i,i_] = v*a[i,i_];
325                    }
326                }
327            }
328        }
329
330
331        public static void rmatrixplu(ref double[,] a,
332            int m,
333            int n,
334            ref int[] pivots)
335        {
336            double[] tmp = new double[0];
337            int i = 0;
338            int j = 0;
339            double mx = 0;
340            double v = 0;
341            int i_ = 0;
342
343           
344            //
345            // Internal LU decomposition subroutine.
346            // Never call it directly.
347            //
348            System.Diagnostics.Debug.Assert(m>0, "RMatrixPLU: incorrect M!");
349            System.Diagnostics.Debug.Assert(n>0, "RMatrixPLU: incorrect N!");
350            tmp = new double[2*Math.Max(m, n)];
351            pivots = new int[Math.Min(m, n)];
352           
353            //
354            // Scale matrix to avoid overflows,
355            // decompose it, then scale back.
356            //
357            mx = 0;
358            for(i=0; i<=m-1; i++)
359            {
360                for(j=0; j<=n-1; j++)
361                {
362                    mx = Math.Max(mx, Math.Abs(a[i,j]));
363                }
364            }
365            if( (double)(mx)!=(double)(0) )
366            {
367                v = 1/mx;
368                for(i=0; i<=m-1; i++)
369                {
370                    for(i_=0; i_<=n-1;i_++)
371                    {
372                        a[i,i_] = v*a[i,i_];
373                    }
374                }
375            }
376            rmatrixplurec(ref a, 0, m, n, ref pivots, ref tmp);
377            if( (double)(mx)!=(double)(0) )
378            {
379                v = mx;
380                for(i=0; i<=Math.Min(m, n)-1; i++)
381                {
382                    for(i_=i; i_<=n-1;i_++)
383                    {
384                        a[i,i_] = v*a[i,i_];
385                    }
386                }
387            }
388        }
389
390
391        public static void cmatrixplu(ref AP.Complex[,] a,
392            int m,
393            int n,
394            ref int[] pivots)
395        {
396            AP.Complex[] tmp = new AP.Complex[0];
397            int i = 0;
398            int j = 0;
399            double mx = 0;
400            AP.Complex v = 0;
401            int i_ = 0;
402
403           
404            //
405            // Internal LU decomposition subroutine.
406            // Never call it directly.
407            //
408            System.Diagnostics.Debug.Assert(m>0, "CMatrixPLU: incorrect M!");
409            System.Diagnostics.Debug.Assert(n>0, "CMatrixPLU: incorrect N!");
410            tmp = new AP.Complex[2*Math.Max(m, n)];
411            pivots = new int[Math.Min(m, n)];
412           
413            //
414            // Scale matrix to avoid overflows,
415            // decompose it, then scale back.
416            //
417            mx = 0;
418            for(i=0; i<=m-1; i++)
419            {
420                for(j=0; j<=n-1; j++)
421                {
422                    mx = Math.Max(mx, AP.Math.AbsComplex(a[i,j]));
423                }
424            }
425            if( (double)(mx)!=(double)(0) )
426            {
427                v = 1/mx;
428                for(i=0; i<=m-1; i++)
429                {
430                    for(i_=0; i_<=n-1;i_++)
431                    {
432                        a[i,i_] = v*a[i,i_];
433                    }
434                }
435            }
436            cmatrixplurec(ref a, 0, m, n, ref pivots, ref tmp);
437            if( (double)(mx)!=(double)(0) )
438            {
439                v = mx;
440                for(i=0; i<=Math.Min(m, n)-1; i++)
441                {
442                    for(i_=i; i_<=n-1;i_++)
443                    {
444                        a[i,i_] = v*a[i,i_];
445                    }
446                }
447            }
448        }
449
450
451        /*************************************************************************
452        Recurrent complex LU subroutine.
453        Never call it directly.
454
455          -- ALGLIB routine --
456             04.01.2010
457             Bochkanov Sergey
458        *************************************************************************/
459        private static void cmatrixluprec(ref AP.Complex[,] a,
460            int offs,
461            int m,
462            int n,
463            ref int[] pivots,
464            ref AP.Complex[] tmp)
465        {
466            int i = 0;
467            int m1 = 0;
468            int m2 = 0;
469            int i_ = 0;
470            int i1_ = 0;
471
472           
473            //
474            // Kernel case
475            //
476            if( Math.Min(m, n)<=ablas.ablascomplexblocksize(ref a) )
477            {
478                cmatrixlup2(ref a, offs, m, n, ref pivots, ref tmp);
479                return;
480            }
481           
482            //
483            // Preliminary step, make N>=M
484            //
485            //     ( A1 )
486            // A = (    ), where A1 is square
487            //     ( A2 )
488            //
489            // Factorize A1, update A2
490            //
491            if( m>n )
492            {
493                cmatrixluprec(ref a, offs, n, n, ref pivots, ref tmp);
494                for(i=0; i<=n-1; i++)
495                {
496                    i1_ = (offs+n) - (0);
497                    for(i_=0; i_<=m-n-1;i_++)
498                    {
499                        tmp[i_] = a[i_+i1_,offs+i];
500                    }
501                    for(i_=offs+n; i_<=offs+m-1;i_++)
502                    {
503                        a[i_,offs+i] = a[i_,pivots[offs+i]];
504                    }
505                    i1_ = (0) - (offs+n);
506                    for(i_=offs+n; i_<=offs+m-1;i_++)
507                    {
508                        a[i_,pivots[offs+i]] = tmp[i_+i1_];
509                    }
510                }
511                ablas.cmatrixrighttrsm(m-n, n, ref a, offs, offs, true, true, 0, ref a, offs+n, offs);
512                return;
513            }
514           
515            //
516            // Non-kernel case
517            //
518            ablas.ablascomplexsplitlength(ref a, m, ref m1, ref m2);
519            cmatrixluprec(ref a, offs, m1, n, ref pivots, ref tmp);
520            if( m2>0 )
521            {
522                for(i=0; i<=m1-1; i++)
523                {
524                    if( offs+i!=pivots[offs+i] )
525                    {
526                        i1_ = (offs+m1) - (0);
527                        for(i_=0; i_<=m2-1;i_++)
528                        {
529                            tmp[i_] = a[i_+i1_,offs+i];
530                        }
531                        for(i_=offs+m1; i_<=offs+m-1;i_++)
532                        {
533                            a[i_,offs+i] = a[i_,pivots[offs+i]];
534                        }
535                        i1_ = (0) - (offs+m1);
536                        for(i_=offs+m1; i_<=offs+m-1;i_++)
537                        {
538                            a[i_,pivots[offs+i]] = tmp[i_+i1_];
539                        }
540                    }
541                }
542                ablas.cmatrixrighttrsm(m2, m1, ref a, offs, offs, true, true, 0, ref a, offs+m1, offs);
543                ablas.cmatrixgemm(m-m1, n-m1, m1, -1.0, ref a, offs+m1, offs, 0, ref a, offs, offs+m1, 0, +1.0, ref a, offs+m1, offs+m1);
544                cmatrixluprec(ref a, offs+m1, m-m1, n-m1, ref pivots, ref tmp);
545                for(i=0; i<=m2-1; i++)
546                {
547                    if( offs+m1+i!=pivots[offs+m1+i] )
548                    {
549                        i1_ = (offs) - (0);
550                        for(i_=0; i_<=m1-1;i_++)
551                        {
552                            tmp[i_] = a[i_+i1_,offs+m1+i];
553                        }
554                        for(i_=offs; i_<=offs+m1-1;i_++)
555                        {
556                            a[i_,offs+m1+i] = a[i_,pivots[offs+m1+i]];
557                        }
558                        i1_ = (0) - (offs);
559                        for(i_=offs; i_<=offs+m1-1;i_++)
560                        {
561                            a[i_,pivots[offs+m1+i]] = tmp[i_+i1_];
562                        }
563                    }
564                }
565            }
566        }
567
568
569        /*************************************************************************
570        Recurrent real LU subroutine.
571        Never call it directly.
572
573          -- ALGLIB routine --
574             04.01.2010
575             Bochkanov Sergey
576        *************************************************************************/
577        private static void rmatrixluprec(ref double[,] a,
578            int offs,
579            int m,
580            int n,
581            ref int[] pivots,
582            ref double[] tmp)
583        {
584            int i = 0;
585            int m1 = 0;
586            int m2 = 0;
587            int i_ = 0;
588            int i1_ = 0;
589
590           
591            //
592            // Kernel case
593            //
594            if( Math.Min(m, n)<=ablas.ablasblocksize(ref a) )
595            {
596                rmatrixlup2(ref a, offs, m, n, ref pivots, ref tmp);
597                return;
598            }
599           
600            //
601            // Preliminary step, make N>=M
602            //
603            //     ( A1 )
604            // A = (    ), where A1 is square
605            //     ( A2 )
606            //
607            // Factorize A1, update A2
608            //
609            if( m>n )
610            {
611                rmatrixluprec(ref a, offs, n, n, ref pivots, ref tmp);
612                for(i=0; i<=n-1; i++)
613                {
614                    if( offs+i!=pivots[offs+i] )
615                    {
616                        i1_ = (offs+n) - (0);
617                        for(i_=0; i_<=m-n-1;i_++)
618                        {
619                            tmp[i_] = a[i_+i1_,offs+i];
620                        }
621                        for(i_=offs+n; i_<=offs+m-1;i_++)
622                        {
623                            a[i_,offs+i] = a[i_,pivots[offs+i]];
624                        }
625                        i1_ = (0) - (offs+n);
626                        for(i_=offs+n; i_<=offs+m-1;i_++)
627                        {
628                            a[i_,pivots[offs+i]] = tmp[i_+i1_];
629                        }
630                    }
631                }
632                ablas.rmatrixrighttrsm(m-n, n, ref a, offs, offs, true, true, 0, ref a, offs+n, offs);
633                return;
634            }
635           
636            //
637            // Non-kernel case
638            //
639            ablas.ablassplitlength(ref a, m, ref m1, ref m2);
640            rmatrixluprec(ref a, offs, m1, n, ref pivots, ref tmp);
641            if( m2>0 )
642            {
643                for(i=0; i<=m1-1; i++)
644                {
645                    if( offs+i!=pivots[offs+i] )
646                    {
647                        i1_ = (offs+m1) - (0);
648                        for(i_=0; i_<=m2-1;i_++)
649                        {
650                            tmp[i_] = a[i_+i1_,offs+i];
651                        }
652                        for(i_=offs+m1; i_<=offs+m-1;i_++)
653                        {
654                            a[i_,offs+i] = a[i_,pivots[offs+i]];
655                        }
656                        i1_ = (0) - (offs+m1);
657                        for(i_=offs+m1; i_<=offs+m-1;i_++)
658                        {
659                            a[i_,pivots[offs+i]] = tmp[i_+i1_];
660                        }
661                    }
662                }
663                ablas.rmatrixrighttrsm(m2, m1, ref a, offs, offs, true, true, 0, ref a, offs+m1, offs);
664                ablas.rmatrixgemm(m-m1, n-m1, m1, -1.0, ref a, offs+m1, offs, 0, ref a, offs, offs+m1, 0, +1.0, ref a, offs+m1, offs+m1);
665                rmatrixluprec(ref a, offs+m1, m-m1, n-m1, ref pivots, ref tmp);
666                for(i=0; i<=m2-1; i++)
667                {
668                    if( offs+m1+i!=pivots[offs+m1+i] )
669                    {
670                        i1_ = (offs) - (0);
671                        for(i_=0; i_<=m1-1;i_++)
672                        {
673                            tmp[i_] = a[i_+i1_,offs+m1+i];
674                        }
675                        for(i_=offs; i_<=offs+m1-1;i_++)
676                        {
677                            a[i_,offs+m1+i] = a[i_,pivots[offs+m1+i]];
678                        }
679                        i1_ = (0) - (offs);
680                        for(i_=offs; i_<=offs+m1-1;i_++)
681                        {
682                            a[i_,pivots[offs+m1+i]] = tmp[i_+i1_];
683                        }
684                    }
685                }
686            }
687        }
688
689
690        /*************************************************************************
691        Recurrent complex LU subroutine.
692        Never call it directly.
693
694          -- ALGLIB routine --
695             04.01.2010
696             Bochkanov Sergey
697        *************************************************************************/
698        private static void cmatrixplurec(ref AP.Complex[,] a,
699            int offs,
700            int m,
701            int n,
702            ref int[] pivots,
703            ref AP.Complex[] tmp)
704        {
705            int i = 0;
706            int n1 = 0;
707            int n2 = 0;
708            int i_ = 0;
709            int i1_ = 0;
710
711           
712            //
713            // Kernel case
714            //
715            if( Math.Min(m, n)<=ablas.ablascomplexblocksize(ref a) )
716            {
717                cmatrixplu2(ref a, offs, m, n, ref pivots, ref tmp);
718                return;
719            }
720           
721            //
722            // Preliminary step, make M>=N.
723            //
724            // A = (A1 A2), where A1 is square
725            // Factorize A1, update A2
726            //
727            if( n>m )
728            {
729                cmatrixplurec(ref a, offs, m, m, ref pivots, ref tmp);
730                for(i=0; i<=m-1; i++)
731                {
732                    i1_ = (offs+m) - (0);
733                    for(i_=0; i_<=n-m-1;i_++)
734                    {
735                        tmp[i_] = a[offs+i,i_+i1_];
736                    }
737                    for(i_=offs+m; i_<=offs+n-1;i_++)
738                    {
739                        a[offs+i,i_] = a[pivots[offs+i],i_];
740                    }
741                    i1_ = (0) - (offs+m);
742                    for(i_=offs+m; i_<=offs+n-1;i_++)
743                    {
744                        a[pivots[offs+i],i_] = tmp[i_+i1_];
745                    }
746                }
747                ablas.cmatrixlefttrsm(m, n-m, ref a, offs, offs, false, true, 0, ref a, offs, offs+m);
748                return;
749            }
750           
751            //
752            // Non-kernel case
753            //
754            ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
755            cmatrixplurec(ref a, offs, m, n1, ref pivots, ref tmp);
756            if( n2>0 )
757            {
758                for(i=0; i<=n1-1; i++)
759                {
760                    if( offs+i!=pivots[offs+i] )
761                    {
762                        i1_ = (offs+n1) - (0);
763                        for(i_=0; i_<=n2-1;i_++)
764                        {
765                            tmp[i_] = a[offs+i,i_+i1_];
766                        }
767                        for(i_=offs+n1; i_<=offs+n-1;i_++)
768                        {
769                            a[offs+i,i_] = a[pivots[offs+i],i_];
770                        }
771                        i1_ = (0) - (offs+n1);
772                        for(i_=offs+n1; i_<=offs+n-1;i_++)
773                        {
774                            a[pivots[offs+i],i_] = tmp[i_+i1_];
775                        }
776                    }
777                }
778                ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, false, true, 0, ref a, offs, offs+n1);
779                ablas.cmatrixgemm(m-n1, n-n1, n1, -1.0, ref a, offs+n1, offs, 0, ref a, offs, offs+n1, 0, +1.0, ref a, offs+n1, offs+n1);
780                cmatrixplurec(ref a, offs+n1, m-n1, n-n1, ref pivots, ref tmp);
781                for(i=0; i<=n2-1; i++)
782                {
783                    if( offs+n1+i!=pivots[offs+n1+i] )
784                    {
785                        i1_ = (offs) - (0);
786                        for(i_=0; i_<=n1-1;i_++)
787                        {
788                            tmp[i_] = a[offs+n1+i,i_+i1_];
789                        }
790                        for(i_=offs; i_<=offs+n1-1;i_++)
791                        {
792                            a[offs+n1+i,i_] = a[pivots[offs+n1+i],i_];
793                        }
794                        i1_ = (0) - (offs);
795                        for(i_=offs; i_<=offs+n1-1;i_++)
796                        {
797                            a[pivots[offs+n1+i],i_] = tmp[i_+i1_];
798                        }
799                    }
800                }
801            }
802        }
803
804
805        /*************************************************************************
806        Recurrent real LU subroutine.
807        Never call it directly.
808
809          -- ALGLIB routine --
810             04.01.2010
811             Bochkanov Sergey
812        *************************************************************************/
813        private static void rmatrixplurec(ref double[,] a,
814            int offs,
815            int m,
816            int n,
817            ref int[] pivots,
818            ref double[] tmp)
819        {
820            int i = 0;
821            int n1 = 0;
822            int n2 = 0;
823            int i_ = 0;
824            int i1_ = 0;
825
826           
827            //
828            // Kernel case
829            //
830            if( Math.Min(m, n)<=ablas.ablasblocksize(ref a) )
831            {
832                rmatrixplu2(ref a, offs, m, n, ref pivots, ref tmp);
833                return;
834            }
835           
836            //
837            // Preliminary step, make M>=N.
838            //
839            // A = (A1 A2), where A1 is square
840            // Factorize A1, update A2
841            //
842            if( n>m )
843            {
844                rmatrixplurec(ref a, offs, m, m, ref pivots, ref tmp);
845                for(i=0; i<=m-1; i++)
846                {
847                    i1_ = (offs+m) - (0);
848                    for(i_=0; i_<=n-m-1;i_++)
849                    {
850                        tmp[i_] = a[offs+i,i_+i1_];
851                    }
852                    for(i_=offs+m; i_<=offs+n-1;i_++)
853                    {
854                        a[offs+i,i_] = a[pivots[offs+i],i_];
855                    }
856                    i1_ = (0) - (offs+m);
857                    for(i_=offs+m; i_<=offs+n-1;i_++)
858                    {
859                        a[pivots[offs+i],i_] = tmp[i_+i1_];
860                    }
861                }
862                ablas.rmatrixlefttrsm(m, n-m, ref a, offs, offs, false, true, 0, ref a, offs, offs+m);
863                return;
864            }
865           
866            //
867            // Non-kernel case
868            //
869            ablas.ablassplitlength(ref a, n, ref n1, ref n2);
870            rmatrixplurec(ref a, offs, m, n1, ref pivots, ref tmp);
871            if( n2>0 )
872            {
873                for(i=0; i<=n1-1; i++)
874                {
875                    if( offs+i!=pivots[offs+i] )
876                    {
877                        i1_ = (offs+n1) - (0);
878                        for(i_=0; i_<=n2-1;i_++)
879                        {
880                            tmp[i_] = a[offs+i,i_+i1_];
881                        }
882                        for(i_=offs+n1; i_<=offs+n-1;i_++)
883                        {
884                            a[offs+i,i_] = a[pivots[offs+i],i_];
885                        }
886                        i1_ = (0) - (offs+n1);
887                        for(i_=offs+n1; i_<=offs+n-1;i_++)
888                        {
889                            a[pivots[offs+i],i_] = tmp[i_+i1_];
890                        }
891                    }
892                }
893                ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, false, true, 0, ref a, offs, offs+n1);
894                ablas.rmatrixgemm(m-n1, n-n1, n1, -1.0, ref a, offs+n1, offs, 0, ref a, offs, offs+n1, 0, +1.0, ref a, offs+n1, offs+n1);
895                rmatrixplurec(ref a, offs+n1, m-n1, n-n1, ref pivots, ref tmp);
896                for(i=0; i<=n2-1; i++)
897                {
898                    if( offs+n1+i!=pivots[offs+n1+i] )
899                    {
900                        i1_ = (offs) - (0);
901                        for(i_=0; i_<=n1-1;i_++)
902                        {
903                            tmp[i_] = a[offs+n1+i,i_+i1_];
904                        }
905                        for(i_=offs; i_<=offs+n1-1;i_++)
906                        {
907                            a[offs+n1+i,i_] = a[pivots[offs+n1+i],i_];
908                        }
909                        i1_ = (0) - (offs);
910                        for(i_=offs; i_<=offs+n1-1;i_++)
911                        {
912                            a[pivots[offs+n1+i],i_] = tmp[i_+i1_];
913                        }
914                    }
915                }
916            }
917        }
918
919
920        /*************************************************************************
921        Complex LUP kernel
922
923          -- ALGLIB routine --
924             10.01.2010
925             Bochkanov Sergey
926        *************************************************************************/
927        private static void cmatrixlup2(ref AP.Complex[,] a,
928            int offs,
929            int m,
930            int n,
931            ref int[] pivots,
932            ref AP.Complex[] tmp)
933        {
934            int i = 0;
935            int j = 0;
936            int jp = 0;
937            AP.Complex s = 0;
938            int i_ = 0;
939            int i1_ = 0;
940
941           
942            //
943            // Quick return if possible
944            //
945            if( m==0 | n==0 )
946            {
947                return;
948            }
949           
950            //
951            // main cycle
952            //
953            for(j=0; j<=Math.Min(m-1, n-1); j++)
954            {
955               
956                //
957                // Find pivot, swap columns
958                //
959                jp = j;
960                for(i=j+1; i<=n-1; i++)
961                {
962                    if( (double)(AP.Math.AbsComplex(a[offs+j,offs+i]))>(double)(AP.Math.AbsComplex(a[offs+j,offs+jp])) )
963                    {
964                        jp = i;
965                    }
966                }
967                pivots[offs+j] = offs+jp;
968                if( jp!=j )
969                {
970                    i1_ = (offs) - (0);
971                    for(i_=0; i_<=m-1;i_++)
972                    {
973                        tmp[i_] = a[i_+i1_,offs+j];
974                    }
975                    for(i_=offs; i_<=offs+m-1;i_++)
976                    {
977                        a[i_,offs+j] = a[i_,offs+jp];
978                    }
979                    i1_ = (0) - (offs);
980                    for(i_=offs; i_<=offs+m-1;i_++)
981                    {
982                        a[i_,offs+jp] = tmp[i_+i1_];
983                    }
984                }
985               
986                //
987                // LU decomposition of 1x(N-J) matrix
988                //
989                if( a[offs+j,offs+j]!=0 & j+1<=n-1 )
990                {
991                    s = 1/a[offs+j,offs+j];
992                    for(i_=offs+j+1; i_<=offs+n-1;i_++)
993                    {
994                        a[offs+j,i_] = s*a[offs+j,i_];
995                    }
996                }
997               
998                //
999                // Update trailing (M-J-1)x(N-J-1) matrix
1000                //
1001                if( j<Math.Min(m-1, n-1) )
1002                {
1003                    i1_ = (offs+j+1) - (0);
1004                    for(i_=0; i_<=m-j-2;i_++)
1005                    {
1006                        tmp[i_] = a[i_+i1_,offs+j];
1007                    }
1008                    i1_ = (offs+j+1) - (m);
1009                    for(i_=m; i_<=m+n-j-2;i_++)
1010                    {
1011                        tmp[i_] = -a[offs+j,i_+i1_];
1012                    }
1013                    ablas.cmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1014                }
1015            }
1016        }
1017
1018
1019        /*************************************************************************
1020        Real LUP kernel
1021
1022          -- ALGLIB routine --
1023             10.01.2010
1024             Bochkanov Sergey
1025        *************************************************************************/
1026        private static void rmatrixlup2(ref double[,] a,
1027            int offs,
1028            int m,
1029            int n,
1030            ref int[] pivots,
1031            ref double[] tmp)
1032        {
1033            int i = 0;
1034            int j = 0;
1035            int jp = 0;
1036            double s = 0;
1037            int i_ = 0;
1038            int i1_ = 0;
1039
1040           
1041            //
1042            // Quick return if possible
1043            //
1044            if( m==0 | n==0 )
1045            {
1046                return;
1047            }
1048           
1049            //
1050            // main cycle
1051            //
1052            for(j=0; j<=Math.Min(m-1, n-1); j++)
1053            {
1054               
1055                //
1056                // Find pivot, swap columns
1057                //
1058                jp = j;
1059                for(i=j+1; i<=n-1; i++)
1060                {
1061                    if( (double)(Math.Abs(a[offs+j,offs+i]))>(double)(Math.Abs(a[offs+j,offs+jp])) )
1062                    {
1063                        jp = i;
1064                    }
1065                }
1066                pivots[offs+j] = offs+jp;
1067                if( jp!=j )
1068                {
1069                    i1_ = (offs) - (0);
1070                    for(i_=0; i_<=m-1;i_++)
1071                    {
1072                        tmp[i_] = a[i_+i1_,offs+j];
1073                    }
1074                    for(i_=offs; i_<=offs+m-1;i_++)
1075                    {
1076                        a[i_,offs+j] = a[i_,offs+jp];
1077                    }
1078                    i1_ = (0) - (offs);
1079                    for(i_=offs; i_<=offs+m-1;i_++)
1080                    {
1081                        a[i_,offs+jp] = tmp[i_+i1_];
1082                    }
1083                }
1084               
1085                //
1086                // LU decomposition of 1x(N-J) matrix
1087                //
1088                if( (double)(a[offs+j,offs+j])!=(double)(0) & j+1<=n-1 )
1089                {
1090                    s = 1/a[offs+j,offs+j];
1091                    for(i_=offs+j+1; i_<=offs+n-1;i_++)
1092                    {
1093                        a[offs+j,i_] = s*a[offs+j,i_];
1094                    }
1095                }
1096               
1097                //
1098                // Update trailing (M-J-1)x(N-J-1) matrix
1099                //
1100                if( j<Math.Min(m-1, n-1) )
1101                {
1102                    i1_ = (offs+j+1) - (0);
1103                    for(i_=0; i_<=m-j-2;i_++)
1104                    {
1105                        tmp[i_] = a[i_+i1_,offs+j];
1106                    }
1107                    i1_ = (offs+j+1) - (m);
1108                    for(i_=m; i_<=m+n-j-2;i_++)
1109                    {
1110                        tmp[i_] = -a[offs+j,i_+i1_];
1111                    }
1112                    ablas.rmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1113                }
1114            }
1115        }
1116
1117
1118        /*************************************************************************
1119        Complex PLU kernel
1120
1121          -- LAPACK routine (version 3.0) --
1122             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1123             Courant Institute, Argonne National Lab, and Rice University
1124             June 30, 1992
1125        *************************************************************************/
1126        private static void cmatrixplu2(ref AP.Complex[,] a,
1127            int offs,
1128            int m,
1129            int n,
1130            ref int[] pivots,
1131            ref AP.Complex[] tmp)
1132        {
1133            int i = 0;
1134            int j = 0;
1135            int jp = 0;
1136            AP.Complex s = 0;
1137            int i_ = 0;
1138            int i1_ = 0;
1139
1140           
1141            //
1142            // Quick return if possible
1143            //
1144            if( m==0 | n==0 )
1145            {
1146                return;
1147            }
1148            for(j=0; j<=Math.Min(m-1, n-1); j++)
1149            {
1150               
1151                //
1152                // Find pivot and test for singularity.
1153                //
1154                jp = j;
1155                for(i=j+1; i<=m-1; i++)
1156                {
1157                    if( (double)(AP.Math.AbsComplex(a[offs+i,offs+j]))>(double)(AP.Math.AbsComplex(a[offs+jp,offs+j])) )
1158                    {
1159                        jp = i;
1160                    }
1161                }
1162                pivots[offs+j] = offs+jp;
1163                if( a[offs+jp,offs+j]!=0 )
1164                {
1165                   
1166                    //
1167                    //Apply the interchange to rows
1168                    //
1169                    if( jp!=j )
1170                    {
1171                        for(i=0; i<=n-1; i++)
1172                        {
1173                            s = a[offs+j,offs+i];
1174                            a[offs+j,offs+i] = a[offs+jp,offs+i];
1175                            a[offs+jp,offs+i] = s;
1176                        }
1177                    }
1178                   
1179                    //
1180                    //Compute elements J+1:M of J-th column.
1181                    //
1182                    if( j<m )
1183                    {
1184                        jp = j+1;
1185                        s = 1/a[offs+j,offs+j];
1186                        for(i_=offs+jp; i_<=offs+m-1;i_++)
1187                        {
1188                            a[i_,offs+j] = s*a[i_,offs+j];
1189                        }
1190                    }
1191                }
1192                if( j<Math.Min(m, n)-1 )
1193                {
1194                   
1195                    //
1196                    //Update trailing submatrix.
1197                    //
1198                    i1_ = (offs+j+1) - (0);
1199                    for(i_=0; i_<=m-j-2;i_++)
1200                    {
1201                        tmp[i_] = a[i_+i1_,offs+j];
1202                    }
1203                    i1_ = (offs+j+1) - (m);
1204                    for(i_=m; i_<=m+n-j-2;i_++)
1205                    {
1206                        tmp[i_] = -a[offs+j,i_+i1_];
1207                    }
1208                    ablas.cmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1209                }
1210            }
1211        }
1212
1213
1214        /*************************************************************************
1215        Real PLU kernel
1216
1217          -- LAPACK routine (version 3.0) --
1218             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1219             Courant Institute, Argonne National Lab, and Rice University
1220             June 30, 1992
1221        *************************************************************************/
1222        private static void rmatrixplu2(ref double[,] a,
1223            int offs,
1224            int m,
1225            int n,
1226            ref int[] pivots,
1227            ref double[] tmp)
1228        {
1229            int i = 0;
1230            int j = 0;
1231            int jp = 0;
1232            double s = 0;
1233            int i_ = 0;
1234            int i1_ = 0;
1235
1236           
1237            //
1238            // Quick return if possible
1239            //
1240            if( m==0 | n==0 )
1241            {
1242                return;
1243            }
1244            for(j=0; j<=Math.Min(m-1, n-1); j++)
1245            {
1246               
1247                //
1248                // Find pivot and test for singularity.
1249                //
1250                jp = j;
1251                for(i=j+1; i<=m-1; i++)
1252                {
1253                    if( (double)(Math.Abs(a[offs+i,offs+j]))>(double)(Math.Abs(a[offs+jp,offs+j])) )
1254                    {
1255                        jp = i;
1256                    }
1257                }
1258                pivots[offs+j] = offs+jp;
1259                if( (double)(a[offs+jp,offs+j])!=(double)(0) )
1260                {
1261                   
1262                    //
1263                    //Apply the interchange to rows
1264                    //
1265                    if( jp!=j )
1266                    {
1267                        for(i=0; i<=n-1; i++)
1268                        {
1269                            s = a[offs+j,offs+i];
1270                            a[offs+j,offs+i] = a[offs+jp,offs+i];
1271                            a[offs+jp,offs+i] = s;
1272                        }
1273                    }
1274                   
1275                    //
1276                    //Compute elements J+1:M of J-th column.
1277                    //
1278                    if( j<m )
1279                    {
1280                        jp = j+1;
1281                        s = 1/a[offs+j,offs+j];
1282                        for(i_=offs+jp; i_<=offs+m-1;i_++)
1283                        {
1284                            a[i_,offs+j] = s*a[i_,offs+j];
1285                        }
1286                    }
1287                }
1288                if( j<Math.Min(m, n)-1 )
1289                {
1290                   
1291                    //
1292                    //Update trailing submatrix.
1293                    //
1294                    i1_ = (offs+j+1) - (0);
1295                    for(i_=0; i_<=m-j-2;i_++)
1296                    {
1297                        tmp[i_] = a[i_+i1_,offs+j];
1298                    }
1299                    i1_ = (offs+j+1) - (m);
1300                    for(i_=m; i_<=m+n-j-2;i_++)
1301                    {
1302                        tmp[i_] = -a[offs+j,i_+i1_];
1303                    }
1304                    ablas.rmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1305                }
1306            }
1307        }
1308
1309
1310        /*************************************************************************
1311        Recursive computational subroutine for HPDMatrixCholesky
1312
1313          -- ALGLIB routine --
1314             15.12.2009
1315             Bochkanov Sergey
1316        *************************************************************************/
1317        private static bool hpdmatrixcholeskyrec(ref AP.Complex[,] a,
1318            int offs,
1319            int n,
1320            bool isupper,
1321            ref AP.Complex[] tmp)
1322        {
1323            bool result = new bool();
1324            int n1 = 0;
1325            int n2 = 0;
1326
1327           
1328            //
1329            // check N
1330            //
1331            if( n<1 )
1332            {
1333                result = false;
1334                return result;
1335            }
1336           
1337            //
1338            // special cases
1339            //
1340            if( n==1 )
1341            {
1342                if( (double)(a[offs,offs].x)>(double)(0) )
1343                {
1344                    a[offs,offs] = Math.Sqrt(a[offs,offs].x);
1345                    result = true;
1346                }
1347                else
1348                {
1349                    result = false;
1350                }
1351                return result;
1352            }
1353            if( n<=ablas.ablascomplexblocksize(ref a) )
1354            {
1355                result = hpdmatrixcholesky2(ref a, offs, n, isupper, ref tmp);
1356                return result;
1357            }
1358           
1359            //
1360            // general case: split task in cache-oblivious manner
1361            //
1362            result = true;
1363            ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1364            result = hpdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp);
1365            if( !result )
1366            {
1367                return result;
1368            }
1369            if( n2>0 )
1370            {
1371                if( isupper )
1372                {
1373                    ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 2, ref a, offs, offs+n1);
1374                    ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 2, +1.0, ref a, offs+n1, offs+n1, isupper);
1375                }
1376                else
1377                {
1378                    ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 2, ref a, offs+n1, offs);
1379                    ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper);
1380                }
1381                result = hpdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp);
1382                if( !result )
1383                {
1384                    return result;
1385                }
1386            }
1387            return result;
1388        }
1389
1390
1391        /*************************************************************************
1392        Recursive computational subroutine for SPDMatrixCholesky
1393
1394          -- ALGLIB routine --
1395             15.12.2009
1396             Bochkanov Sergey
1397        *************************************************************************/
1398        private static bool spdmatrixcholeskyrec(ref double[,] a,
1399            int offs,
1400            int n,
1401            bool isupper,
1402            ref double[] tmp)
1403        {
1404            bool result = new bool();
1405            int n1 = 0;
1406            int n2 = 0;
1407
1408           
1409            //
1410            // check N
1411            //
1412            if( n<1 )
1413            {
1414                result = false;
1415                return result;
1416            }
1417           
1418            //
1419            // special cases
1420            //
1421            if( n==1 )
1422            {
1423                if( (double)(a[offs,offs])>(double)(0) )
1424                {
1425                    a[offs,offs] = Math.Sqrt(a[offs,offs]);
1426                    result = true;
1427                }
1428                else
1429                {
1430                    result = false;
1431                }
1432                return result;
1433            }
1434            if( n<=ablas.ablasblocksize(ref a) )
1435            {
1436                result = spdmatrixcholesky2(ref a, offs, n, isupper, ref tmp);
1437                return result;
1438            }
1439           
1440            //
1441            // general case: split task in cache-oblivious manner
1442            //
1443            result = true;
1444            ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1445            result = spdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp);
1446            if( !result )
1447            {
1448                return result;
1449            }
1450            if( n2>0 )
1451            {
1452                if( isupper )
1453                {
1454                    ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 1, ref a, offs, offs+n1);
1455                    ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 1, +1.0, ref a, offs+n1, offs+n1, isupper);
1456                }
1457                else
1458                {
1459                    ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 1, ref a, offs+n1, offs);
1460                    ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper);
1461                }
1462                result = spdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp);
1463                if( !result )
1464                {
1465                    return result;
1466                }
1467            }
1468            return result;
1469        }
1470
1471
1472        /*************************************************************************
1473        Level-2 Hermitian Cholesky subroutine.
1474
1475          -- LAPACK routine (version 3.0) --
1476             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1477             Courant Institute, Argonne National Lab, and Rice University
1478             February 29, 1992
1479        *************************************************************************/
1480        private static bool hpdmatrixcholesky2(ref AP.Complex[,] aaa,
1481            int offs,
1482            int n,
1483            bool isupper,
1484            ref AP.Complex[] tmp)
1485        {
1486            bool result = new bool();
1487            int i = 0;
1488            int j = 0;
1489            int k = 0;
1490            int j1 = 0;
1491            int j2 = 0;
1492            double ajj = 0;
1493            AP.Complex v = 0;
1494            double r = 0;
1495            int i_ = 0;
1496            int i1_ = 0;
1497
1498            result = true;
1499            if( n<0 )
1500            {
1501                result = false;
1502                return result;
1503            }
1504           
1505            //
1506            // Quick return if possible
1507            //
1508            if( n==0 )
1509            {
1510                return result;
1511            }
1512            if( isupper )
1513            {
1514               
1515                //
1516                // Compute the Cholesky factorization A = U'*U.
1517                //
1518                for(j=0; j<=n-1; j++)
1519                {
1520                   
1521                    //
1522                    // Compute U(J,J) and test for non-positive-definiteness.
1523                    //
1524                    v = 0.0;
1525                    for(i_=offs; i_<=offs+j-1;i_++)
1526                    {
1527                        v += AP.Math.Conj(aaa[i_,offs+j])*aaa[i_,offs+j];
1528                    }
1529                    ajj = (aaa[offs+j,offs+j]-v).x;
1530                    if( (double)(ajj)<=(double)(0) )
1531                    {
1532                        aaa[offs+j,offs+j] = ajj;
1533                        result = false;
1534                        return result;
1535                    }
1536                    ajj = Math.Sqrt(ajj);
1537                    aaa[offs+j,offs+j] = ajj;
1538                   
1539                    //
1540                    // Compute elements J+1:N-1 of row J.
1541                    //
1542                    if( j<n-1 )
1543                    {
1544                        if( j>0 )
1545                        {
1546                            i1_ = (offs) - (0);
1547                            for(i_=0; i_<=j-1;i_++)
1548                            {
1549                                tmp[i_] = -AP.Math.Conj(aaa[i_+i1_,offs+j]);
1550                            }
1551                            ablas.cmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n);
1552                            i1_ = (n) - (offs+j+1);
1553                            for(i_=offs+j+1; i_<=offs+n-1;i_++)
1554                            {
1555                                aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_];
1556                            }
1557                        }
1558                        r = 1/ajj;
1559                        for(i_=offs+j+1; i_<=offs+n-1;i_++)
1560                        {
1561                            aaa[offs+j,i_] = r*aaa[offs+j,i_];
1562                        }
1563                    }
1564                }
1565            }
1566            else
1567            {
1568               
1569                //
1570                // Compute the Cholesky factorization A = L*L'.
1571                //
1572                for(j=0; j<=n-1; j++)
1573                {
1574                   
1575                    //
1576                    // Compute L(J+1,J+1) and test for non-positive-definiteness.
1577                    //
1578                    v = 0.0;
1579                    for(i_=offs; i_<=offs+j-1;i_++)
1580                    {
1581                        v += AP.Math.Conj(aaa[offs+j,i_])*aaa[offs+j,i_];
1582                    }
1583                    ajj = (aaa[offs+j,offs+j]-v).x;
1584                    if( (double)(ajj)<=(double)(0) )
1585                    {
1586                        aaa[offs+j,offs+j] = ajj;
1587                        result = false;
1588                        return result;
1589                    }
1590                    ajj = Math.Sqrt(ajj);
1591                    aaa[offs+j,offs+j] = ajj;
1592                   
1593                    //
1594                    // Compute elements J+1:N of column J.
1595                    //
1596                    if( j<n-1 )
1597                    {
1598                        if( j>0 )
1599                        {
1600                            i1_ = (offs) - (0);
1601                            for(i_=0; i_<=j-1;i_++)
1602                            {
1603                                tmp[i_] = AP.Math.Conj(aaa[offs+j,i_+i1_]);
1604                            }
1605                            ablas.cmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n);
1606                            for(i=0; i<=n-j-2; i++)
1607                            {
1608                                aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj;
1609                            }
1610                        }
1611                        else
1612                        {
1613                            for(i=0; i<=n-j-2; i++)
1614                            {
1615                                aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj;
1616                            }
1617                        }
1618                    }
1619                }
1620            }
1621            return result;
1622        }
1623
1624
1625        /*************************************************************************
1626        Level-2 Cholesky subroutine
1627
1628          -- LAPACK routine (version 3.0) --
1629             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1630             Courant Institute, Argonne National Lab, and Rice University
1631             February 29, 1992
1632        *************************************************************************/
1633        private static bool spdmatrixcholesky2(ref double[,] aaa,
1634            int offs,
1635            int n,
1636            bool isupper,
1637            ref double[] tmp)
1638        {
1639            bool result = new bool();
1640            int i = 0;
1641            int j = 0;
1642            int k = 0;
1643            int j1 = 0;
1644            int j2 = 0;
1645            double ajj = 0;
1646            double v = 0;
1647            double r = 0;
1648            int i_ = 0;
1649            int i1_ = 0;
1650
1651            result = true;
1652            if( n<0 )
1653            {
1654                result = false;
1655                return result;
1656            }
1657           
1658            //
1659            // Quick return if possible
1660            //
1661            if( n==0 )
1662            {
1663                return result;
1664            }
1665            if( isupper )
1666            {
1667               
1668                //
1669                // Compute the Cholesky factorization A = U'*U.
1670                //
1671                for(j=0; j<=n-1; j++)
1672                {
1673                   
1674                    //
1675                    // Compute U(J,J) and test for non-positive-definiteness.
1676                    //
1677                    v = 0.0;
1678                    for(i_=offs; i_<=offs+j-1;i_++)
1679                    {
1680                        v += aaa[i_,offs+j]*aaa[i_,offs+j];
1681                    }
1682                    ajj = aaa[offs+j,offs+j]-v;
1683                    if( (double)(ajj)<=(double)(0) )
1684                    {
1685                        aaa[offs+j,offs+j] = ajj;
1686                        result = false;
1687                        return result;
1688                    }
1689                    ajj = Math.Sqrt(ajj);
1690                    aaa[offs+j,offs+j] = ajj;
1691                   
1692                    //
1693                    // Compute elements J+1:N-1 of row J.
1694                    //
1695                    if( j<n-1 )
1696                    {
1697                        if( j>0 )
1698                        {
1699                            i1_ = (offs) - (0);
1700                            for(i_=0; i_<=j-1;i_++)
1701                            {
1702                                tmp[i_] = -aaa[i_+i1_,offs+j];
1703                            }
1704                            ablas.rmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n);
1705                            i1_ = (n) - (offs+j+1);
1706                            for(i_=offs+j+1; i_<=offs+n-1;i_++)
1707                            {
1708                                aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_];
1709                            }
1710                        }
1711                        r = 1/ajj;
1712                        for(i_=offs+j+1; i_<=offs+n-1;i_++)
1713                        {
1714                            aaa[offs+j,i_] = r*aaa[offs+j,i_];
1715                        }
1716                    }
1717                }
1718            }
1719            else
1720            {
1721               
1722                //
1723                // Compute the Cholesky factorization A = L*L'.
1724                //
1725                for(j=0; j<=n-1; j++)
1726                {
1727                   
1728                    //
1729                    // Compute L(J+1,J+1) and test for non-positive-definiteness.
1730                    //
1731                    v = 0.0;
1732                    for(i_=offs; i_<=offs+j-1;i_++)
1733                    {
1734                        v += aaa[offs+j,i_]*aaa[offs+j,i_];
1735                    }
1736                    ajj = aaa[offs+j,offs+j]-v;
1737                    if( (double)(ajj)<=(double)(0) )
1738                    {
1739                        aaa[offs+j,offs+j] = ajj;
1740                        result = false;
1741                        return result;
1742                    }
1743                    ajj = Math.Sqrt(ajj);
1744                    aaa[offs+j,offs+j] = ajj;
1745                   
1746                    //
1747                    // Compute elements J+1:N of column J.
1748                    //
1749                    if( j<n-1 )
1750                    {
1751                        if( j>0 )
1752                        {
1753                            i1_ = (offs) - (0);
1754                            for(i_=0; i_<=j-1;i_++)
1755                            {
1756                                tmp[i_] = aaa[offs+j,i_+i1_];
1757                            }
1758                            ablas.rmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n);
1759                            for(i=0; i<=n-j-2; i++)
1760                            {
1761                                aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj;
1762                            }
1763                        }
1764                        else
1765                        {
1766                            for(i=0; i<=n-j-2; i++)
1767                            {
1768                                aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj;
1769                            }
1770                        }
1771                    }
1772                }
1773            }
1774            return result;
1775        }
1776    }
1777}
Note: See TracBrowser for help on using the repository browser.