Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/trfac.cs @ 4065

Last change on this file since 4065 was 3839, checked in by mkommend, 15 years ago

implemented first version of LR (ticket #1012)
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+1<=m-1 )
1183                    {
1184                        s = 1/a[offs+j,offs+j];
1185                        for(i_=offs+j+1; i_<=offs+m-1;i_++)
1186                        {
1187                            a[i_,offs+j] = s*a[i_,offs+j];
1188                        }
1189                    }
1190                }
1191                if( j<Math.Min(m, n)-1 )
1192                {
1193                   
1194                    //
1195                    //Update trailing submatrix.
1196                    //
1197                    i1_ = (offs+j+1) - (0);
1198                    for(i_=0; i_<=m-j-2;i_++)
1199                    {
1200                        tmp[i_] = a[i_+i1_,offs+j];
1201                    }
1202                    i1_ = (offs+j+1) - (m);
1203                    for(i_=m; i_<=m+n-j-2;i_++)
1204                    {
1205                        tmp[i_] = -a[offs+j,i_+i1_];
1206                    }
1207                    ablas.cmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1208                }
1209            }
1210        }
1211
1212
1213        /*************************************************************************
1214        Real PLU kernel
1215
1216          -- LAPACK routine (version 3.0) --
1217             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1218             Courant Institute, Argonne National Lab, and Rice University
1219             June 30, 1992
1220        *************************************************************************/
1221        private static void rmatrixplu2(ref double[,] a,
1222            int offs,
1223            int m,
1224            int n,
1225            ref int[] pivots,
1226            ref double[] tmp)
1227        {
1228            int i = 0;
1229            int j = 0;
1230            int jp = 0;
1231            double s = 0;
1232            int i_ = 0;
1233            int i1_ = 0;
1234
1235           
1236            //
1237            // Quick return if possible
1238            //
1239            if( m==0 | n==0 )
1240            {
1241                return;
1242            }
1243            for(j=0; j<=Math.Min(m-1, n-1); j++)
1244            {
1245               
1246                //
1247                // Find pivot and test for singularity.
1248                //
1249                jp = j;
1250                for(i=j+1; i<=m-1; i++)
1251                {
1252                    if( (double)(Math.Abs(a[offs+i,offs+j]))>(double)(Math.Abs(a[offs+jp,offs+j])) )
1253                    {
1254                        jp = i;
1255                    }
1256                }
1257                pivots[offs+j] = offs+jp;
1258                if( (double)(a[offs+jp,offs+j])!=(double)(0) )
1259                {
1260                   
1261                    //
1262                    //Apply the interchange to rows
1263                    //
1264                    if( jp!=j )
1265                    {
1266                        for(i=0; i<=n-1; i++)
1267                        {
1268                            s = a[offs+j,offs+i];
1269                            a[offs+j,offs+i] = a[offs+jp,offs+i];
1270                            a[offs+jp,offs+i] = s;
1271                        }
1272                    }
1273                   
1274                    //
1275                    //Compute elements J+1:M of J-th column.
1276                    //
1277                    if( j+1<=m-1 )
1278                    {
1279                        s = 1/a[offs+j,offs+j];
1280                        for(i_=offs+j+1; i_<=offs+m-1;i_++)
1281                        {
1282                            a[i_,offs+j] = s*a[i_,offs+j];
1283                        }
1284                    }
1285                }
1286                if( j<Math.Min(m, n)-1 )
1287                {
1288                   
1289                    //
1290                    //Update trailing submatrix.
1291                    //
1292                    i1_ = (offs+j+1) - (0);
1293                    for(i_=0; i_<=m-j-2;i_++)
1294                    {
1295                        tmp[i_] = a[i_+i1_,offs+j];
1296                    }
1297                    i1_ = (offs+j+1) - (m);
1298                    for(i_=m; i_<=m+n-j-2;i_++)
1299                    {
1300                        tmp[i_] = -a[offs+j,i_+i1_];
1301                    }
1302                    ablas.rmatrixrank1(m-j-1, n-j-1, ref a, offs+j+1, offs+j+1, ref tmp, 0, ref tmp, m);
1303                }
1304            }
1305        }
1306
1307
1308        /*************************************************************************
1309        Recursive computational subroutine for HPDMatrixCholesky
1310
1311          -- ALGLIB routine --
1312             15.12.2009
1313             Bochkanov Sergey
1314        *************************************************************************/
1315        private static bool hpdmatrixcholeskyrec(ref AP.Complex[,] a,
1316            int offs,
1317            int n,
1318            bool isupper,
1319            ref AP.Complex[] tmp)
1320        {
1321            bool result = new bool();
1322            int n1 = 0;
1323            int n2 = 0;
1324
1325           
1326            //
1327            // check N
1328            //
1329            if( n<1 )
1330            {
1331                result = false;
1332                return result;
1333            }
1334           
1335            //
1336            // special cases
1337            //
1338            if( n==1 )
1339            {
1340                if( (double)(a[offs,offs].x)>(double)(0) )
1341                {
1342                    a[offs,offs] = Math.Sqrt(a[offs,offs].x);
1343                    result = true;
1344                }
1345                else
1346                {
1347                    result = false;
1348                }
1349                return result;
1350            }
1351            if( n<=ablas.ablascomplexblocksize(ref a) )
1352            {
1353                result = hpdmatrixcholesky2(ref a, offs, n, isupper, ref tmp);
1354                return result;
1355            }
1356           
1357            //
1358            // general case: split task in cache-oblivious manner
1359            //
1360            result = true;
1361            ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1362            result = hpdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp);
1363            if( !result )
1364            {
1365                return result;
1366            }
1367            if( n2>0 )
1368            {
1369                if( isupper )
1370                {
1371                    ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 2, ref a, offs, offs+n1);
1372                    ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 2, +1.0, ref a, offs+n1, offs+n1, isupper);
1373                }
1374                else
1375                {
1376                    ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 2, ref a, offs+n1, offs);
1377                    ablas.cmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper);
1378                }
1379                result = hpdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp);
1380                if( !result )
1381                {
1382                    return result;
1383                }
1384            }
1385            return result;
1386        }
1387
1388
1389        /*************************************************************************
1390        Recursive computational subroutine for SPDMatrixCholesky
1391
1392          -- ALGLIB routine --
1393             15.12.2009
1394             Bochkanov Sergey
1395        *************************************************************************/
1396        private static bool spdmatrixcholeskyrec(ref double[,] a,
1397            int offs,
1398            int n,
1399            bool isupper,
1400            ref double[] tmp)
1401        {
1402            bool result = new bool();
1403            int n1 = 0;
1404            int n2 = 0;
1405
1406           
1407            //
1408            // check N
1409            //
1410            if( n<1 )
1411            {
1412                result = false;
1413                return result;
1414            }
1415           
1416            //
1417            // special cases
1418            //
1419            if( n==1 )
1420            {
1421                if( (double)(a[offs,offs])>(double)(0) )
1422                {
1423                    a[offs,offs] = Math.Sqrt(a[offs,offs]);
1424                    result = true;
1425                }
1426                else
1427                {
1428                    result = false;
1429                }
1430                return result;
1431            }
1432            if( n<=ablas.ablasblocksize(ref a) )
1433            {
1434                result = spdmatrixcholesky2(ref a, offs, n, isupper, ref tmp);
1435                return result;
1436            }
1437           
1438            //
1439            // general case: split task in cache-oblivious manner
1440            //
1441            result = true;
1442            ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1443            result = spdmatrixcholeskyrec(ref a, offs, n1, isupper, ref tmp);
1444            if( !result )
1445            {
1446                return result;
1447            }
1448            if( n2>0 )
1449            {
1450                if( isupper )
1451                {
1452                    ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 1, ref a, offs, offs+n1);
1453                    ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs, offs+n1, 1, +1.0, ref a, offs+n1, offs+n1, isupper);
1454                }
1455                else
1456                {
1457                    ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 1, ref a, offs+n1, offs);
1458                    ablas.rmatrixsyrk(n2, n1, -1.0, ref a, offs+n1, offs, 0, +1.0, ref a, offs+n1, offs+n1, isupper);
1459                }
1460                result = spdmatrixcholeskyrec(ref a, offs+n1, n2, isupper, ref tmp);
1461                if( !result )
1462                {
1463                    return result;
1464                }
1465            }
1466            return result;
1467        }
1468
1469
1470        /*************************************************************************
1471        Level-2 Hermitian Cholesky subroutine.
1472
1473          -- LAPACK routine (version 3.0) --
1474             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1475             Courant Institute, Argonne National Lab, and Rice University
1476             February 29, 1992
1477        *************************************************************************/
1478        private static bool hpdmatrixcholesky2(ref AP.Complex[,] aaa,
1479            int offs,
1480            int n,
1481            bool isupper,
1482            ref AP.Complex[] tmp)
1483        {
1484            bool result = new bool();
1485            int i = 0;
1486            int j = 0;
1487            int k = 0;
1488            int j1 = 0;
1489            int j2 = 0;
1490            double ajj = 0;
1491            AP.Complex v = 0;
1492            double r = 0;
1493            int i_ = 0;
1494            int i1_ = 0;
1495
1496            result = true;
1497            if( n<0 )
1498            {
1499                result = false;
1500                return result;
1501            }
1502           
1503            //
1504            // Quick return if possible
1505            //
1506            if( n==0 )
1507            {
1508                return result;
1509            }
1510            if( isupper )
1511            {
1512               
1513                //
1514                // Compute the Cholesky factorization A = U'*U.
1515                //
1516                for(j=0; j<=n-1; j++)
1517                {
1518                   
1519                    //
1520                    // Compute U(J,J) and test for non-positive-definiteness.
1521                    //
1522                    v = 0.0;
1523                    for(i_=offs; i_<=offs+j-1;i_++)
1524                    {
1525                        v += AP.Math.Conj(aaa[i_,offs+j])*aaa[i_,offs+j];
1526                    }
1527                    ajj = (aaa[offs+j,offs+j]-v).x;
1528                    if( (double)(ajj)<=(double)(0) )
1529                    {
1530                        aaa[offs+j,offs+j] = ajj;
1531                        result = false;
1532                        return result;
1533                    }
1534                    ajj = Math.Sqrt(ajj);
1535                    aaa[offs+j,offs+j] = ajj;
1536                   
1537                    //
1538                    // Compute elements J+1:N-1 of row J.
1539                    //
1540                    if( j<n-1 )
1541                    {
1542                        if( j>0 )
1543                        {
1544                            i1_ = (offs) - (0);
1545                            for(i_=0; i_<=j-1;i_++)
1546                            {
1547                                tmp[i_] = -AP.Math.Conj(aaa[i_+i1_,offs+j]);
1548                            }
1549                            ablas.cmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n);
1550                            i1_ = (n) - (offs+j+1);
1551                            for(i_=offs+j+1; i_<=offs+n-1;i_++)
1552                            {
1553                                aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_];
1554                            }
1555                        }
1556                        r = 1/ajj;
1557                        for(i_=offs+j+1; i_<=offs+n-1;i_++)
1558                        {
1559                            aaa[offs+j,i_] = r*aaa[offs+j,i_];
1560                        }
1561                    }
1562                }
1563            }
1564            else
1565            {
1566               
1567                //
1568                // Compute the Cholesky factorization A = L*L'.
1569                //
1570                for(j=0; j<=n-1; j++)
1571                {
1572                   
1573                    //
1574                    // Compute L(J+1,J+1) and test for non-positive-definiteness.
1575                    //
1576                    v = 0.0;
1577                    for(i_=offs; i_<=offs+j-1;i_++)
1578                    {
1579                        v += AP.Math.Conj(aaa[offs+j,i_])*aaa[offs+j,i_];
1580                    }
1581                    ajj = (aaa[offs+j,offs+j]-v).x;
1582                    if( (double)(ajj)<=(double)(0) )
1583                    {
1584                        aaa[offs+j,offs+j] = ajj;
1585                        result = false;
1586                        return result;
1587                    }
1588                    ajj = Math.Sqrt(ajj);
1589                    aaa[offs+j,offs+j] = ajj;
1590                   
1591                    //
1592                    // Compute elements J+1:N of column J.
1593                    //
1594                    if( j<n-1 )
1595                    {
1596                        if( j>0 )
1597                        {
1598                            i1_ = (offs) - (0);
1599                            for(i_=0; i_<=j-1;i_++)
1600                            {
1601                                tmp[i_] = AP.Math.Conj(aaa[offs+j,i_+i1_]);
1602                            }
1603                            ablas.cmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n);
1604                            for(i=0; i<=n-j-2; i++)
1605                            {
1606                                aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj;
1607                            }
1608                        }
1609                        else
1610                        {
1611                            for(i=0; i<=n-j-2; i++)
1612                            {
1613                                aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj;
1614                            }
1615                        }
1616                    }
1617                }
1618            }
1619            return result;
1620        }
1621
1622
1623        /*************************************************************************
1624        Level-2 Cholesky subroutine
1625
1626          -- LAPACK routine (version 3.0) --
1627             Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1628             Courant Institute, Argonne National Lab, and Rice University
1629             February 29, 1992
1630        *************************************************************************/
1631        private static bool spdmatrixcholesky2(ref double[,] aaa,
1632            int offs,
1633            int n,
1634            bool isupper,
1635            ref double[] tmp)
1636        {
1637            bool result = new bool();
1638            int i = 0;
1639            int j = 0;
1640            int k = 0;
1641            int j1 = 0;
1642            int j2 = 0;
1643            double ajj = 0;
1644            double v = 0;
1645            double r = 0;
1646            int i_ = 0;
1647            int i1_ = 0;
1648
1649            result = true;
1650            if( n<0 )
1651            {
1652                result = false;
1653                return result;
1654            }
1655           
1656            //
1657            // Quick return if possible
1658            //
1659            if( n==0 )
1660            {
1661                return result;
1662            }
1663            if( isupper )
1664            {
1665               
1666                //
1667                // Compute the Cholesky factorization A = U'*U.
1668                //
1669                for(j=0; j<=n-1; j++)
1670                {
1671                   
1672                    //
1673                    // Compute U(J,J) and test for non-positive-definiteness.
1674                    //
1675                    v = 0.0;
1676                    for(i_=offs; i_<=offs+j-1;i_++)
1677                    {
1678                        v += aaa[i_,offs+j]*aaa[i_,offs+j];
1679                    }
1680                    ajj = aaa[offs+j,offs+j]-v;
1681                    if( (double)(ajj)<=(double)(0) )
1682                    {
1683                        aaa[offs+j,offs+j] = ajj;
1684                        result = false;
1685                        return result;
1686                    }
1687                    ajj = Math.Sqrt(ajj);
1688                    aaa[offs+j,offs+j] = ajj;
1689                   
1690                    //
1691                    // Compute elements J+1:N-1 of row J.
1692                    //
1693                    if( j<n-1 )
1694                    {
1695                        if( j>0 )
1696                        {
1697                            i1_ = (offs) - (0);
1698                            for(i_=0; i_<=j-1;i_++)
1699                            {
1700                                tmp[i_] = -aaa[i_+i1_,offs+j];
1701                            }
1702                            ablas.rmatrixmv(n-j-1, j, ref aaa, offs, offs+j+1, 1, ref tmp, 0, ref tmp, n);
1703                            i1_ = (n) - (offs+j+1);
1704                            for(i_=offs+j+1; i_<=offs+n-1;i_++)
1705                            {
1706                                aaa[offs+j,i_] = aaa[offs+j,i_] + tmp[i_+i1_];
1707                            }
1708                        }
1709                        r = 1/ajj;
1710                        for(i_=offs+j+1; i_<=offs+n-1;i_++)
1711                        {
1712                            aaa[offs+j,i_] = r*aaa[offs+j,i_];
1713                        }
1714                    }
1715                }
1716            }
1717            else
1718            {
1719               
1720                //
1721                // Compute the Cholesky factorization A = L*L'.
1722                //
1723                for(j=0; j<=n-1; j++)
1724                {
1725                   
1726                    //
1727                    // Compute L(J+1,J+1) and test for non-positive-definiteness.
1728                    //
1729                    v = 0.0;
1730                    for(i_=offs; i_<=offs+j-1;i_++)
1731                    {
1732                        v += aaa[offs+j,i_]*aaa[offs+j,i_];
1733                    }
1734                    ajj = aaa[offs+j,offs+j]-v;
1735                    if( (double)(ajj)<=(double)(0) )
1736                    {
1737                        aaa[offs+j,offs+j] = ajj;
1738                        result = false;
1739                        return result;
1740                    }
1741                    ajj = Math.Sqrt(ajj);
1742                    aaa[offs+j,offs+j] = ajj;
1743                   
1744                    //
1745                    // Compute elements J+1:N of column J.
1746                    //
1747                    if( j<n-1 )
1748                    {
1749                        if( j>0 )
1750                        {
1751                            i1_ = (offs) - (0);
1752                            for(i_=0; i_<=j-1;i_++)
1753                            {
1754                                tmp[i_] = aaa[offs+j,i_+i1_];
1755                            }
1756                            ablas.rmatrixmv(n-j-1, j, ref aaa, offs+j+1, offs, 0, ref tmp, 0, ref tmp, n);
1757                            for(i=0; i<=n-j-2; i++)
1758                            {
1759                                aaa[offs+j+1+i,offs+j] = (aaa[offs+j+1+i,offs+j]-tmp[n+i])/ajj;
1760                            }
1761                        }
1762                        else
1763                        {
1764                            for(i=0; i<=n-j-2; i++)
1765                            {
1766                                aaa[offs+j+1+i,offs+j] = aaa[offs+j+1+i,offs+j]/ajj;
1767                            }
1768                        }
1769                    }
1770                }
1771            }
1772            return result;
1773        }
1774    }
1775}
Note: See TracBrowser for help on using the repository browser.