Free cookie consent management tool by TermsFeed Policy Generator

source: branches/ParameterBinding/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/matinv.cs @ 6451

Last change on this file since 6451 was 4068, checked in by swagner, 14 years ago

Sorted usings and removed unused usings in entire solution (#1094)
File size: 49.7 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
27
28namespace alglib {
29  public class matinv {
30    public struct matinvreport {
31      public double r1;
32      public double rinf;
33    };
34
35
36
37
38    /*************************************************************************
39    Inversion of a matrix given by its LU decomposition.
40
41    INPUT PARAMETERS:
42        A       -   LU decomposition of the matrix (output of RMatrixLU subroutine).
43        Pivots  -   table of permutations which were made during the LU decomposition
44                    (the output of RMatrixLU subroutine).
45        N       -   size of matrix A.
46
47    OUTPUT PARAMETERS:
48        Info    -   return code:
49                    * -3    A is singular, or VERY close to singular.
50                            it is filled by zeros in such cases.
51                    * -1    N<=0 was passed, or incorrect Pivots was passed
52                    *  1    task is solved (but matrix A may be ill-conditioned,
53                            check R1/RInf parameters for condition numbers).
54        Rep     -   solver report, see below for more info
55        A       -   inverse of matrix A.
56                    Array whose indexes range within [0..N-1, 0..N-1].
57
58    SOLVER REPORT
59
60    Subroutine sets following fields of the Rep structure:
61    * R1        reciprocal of condition number: 1/cond(A), 1-norm.
62    * RInf      reciprocal of condition number: 1/cond(A), inf-norm.
63
64      -- ALGLIB routine --
65         05.02.2010
66         Bochkanov Sergey
67    *************************************************************************/
68    public static void rmatrixluinverse(ref double[,] a,
69        ref int[] pivots,
70        int n,
71        ref int info,
72        ref matinvreport rep) {
73      double[] work = new double[0];
74      int i = 0;
75      int j = 0;
76      int k = 0;
77      double v = 0;
78
79      info = 1;
80
81      //
82      // Quick return if possible
83      //
84      if (n == 0) {
85        info = -1;
86        return;
87      }
88      for (i = 0; i <= n - 1; i++) {
89        if (pivots[i] > n - 1 | pivots[i] < i) {
90          info = -1;
91          return;
92        }
93      }
94
95      //
96      // calculate condition numbers
97      //
98      rep.r1 = rcond.rmatrixlurcond1(ref a, n);
99      rep.rinf = rcond.rmatrixlurcondinf(ref a, n);
100      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
101        for (i = 0; i <= n - 1; i++) {
102          for (j = 0; j <= n - 1; j++) {
103            a[i, j] = 0;
104          }
105        }
106        rep.r1 = 0;
107        rep.rinf = 0;
108        info = -3;
109        return;
110      }
111
112      //
113      // Call cache-oblivious code
114      //
115      work = new double[n];
116      rmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);
117
118      //
119      // apply permutations
120      //
121      for (i = 0; i <= n - 1; i++) {
122        for (j = n - 2; j >= 0; j--) {
123          k = pivots[j];
124          v = a[i, j];
125          a[i, j] = a[i, k];
126          a[i, k] = v;
127        }
128      }
129    }
130
131
132    /*************************************************************************
133    Inversion of a general matrix.
134
135    Input parameters:
136        A   -   matrix. Array whose indexes range within [0..N-1, 0..N-1].
137        N   -   size of matrix A.
138
139    Output parameters:
140        Info    -   return code, same as in RMatrixLUInverse
141        Rep     -   solver report, same as in RMatrixLUInverse
142        A       -   inverse of matrix A, same as in RMatrixLUInverse
143
144    Result:
145        True, if the matrix is not singular.
146        False, if the matrix is singular.
147
148      -- ALGLIB --
149         Copyright 2005 by Bochkanov Sergey
150    *************************************************************************/
151    public static void rmatrixinverse(ref double[,] a,
152        int n,
153        ref int info,
154        ref matinvreport rep) {
155      int[] pivots = new int[0];
156
157      trfac.rmatrixlu(ref a, n, n, ref pivots);
158      rmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
159    }
160
161
162    /*************************************************************************
163    Inversion of a matrix given by its LU decomposition.
164
165    INPUT PARAMETERS:
166        A       -   LU decomposition of the matrix (output of CMatrixLU subroutine).
167        Pivots  -   table of permutations which were made during the LU decomposition
168                    (the output of CMatrixLU subroutine).
169        N       -   size of matrix A.
170
171    OUTPUT PARAMETERS:
172        Info    -   return code, same as in RMatrixLUInverse
173        Rep     -   solver report, same as in RMatrixLUInverse
174        A       -   inverse of matrix A, same as in RMatrixLUInverse
175
176      -- ALGLIB routine --
177         05.02.2010
178         Bochkanov Sergey
179    *************************************************************************/
180    public static void cmatrixluinverse(ref AP.Complex[,] a,
181        ref int[] pivots,
182        int n,
183        ref int info,
184        ref matinvreport rep) {
185      AP.Complex[] work = new AP.Complex[0];
186      int i = 0;
187      int j = 0;
188      int k = 0;
189      AP.Complex v = 0;
190
191      info = 1;
192
193      //
194      // Quick return if possible
195      //
196      if (n == 0) {
197        info = -1;
198        return;
199      }
200      for (i = 0; i <= n - 1; i++) {
201        if (pivots[i] > n - 1 | pivots[i] < i) {
202          info = -1;
203          return;
204        }
205      }
206
207      //
208      // calculate condition numbers
209      //
210      rep.r1 = rcond.cmatrixlurcond1(ref a, n);
211      rep.rinf = rcond.cmatrixlurcondinf(ref a, n);
212      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
213        for (i = 0; i <= n - 1; i++) {
214          for (j = 0; j <= n - 1; j++) {
215            a[i, j] = 0;
216          }
217        }
218        rep.r1 = 0;
219        rep.rinf = 0;
220        info = -3;
221        return;
222      }
223
224      //
225      // Call cache-oblivious code
226      //
227      work = new AP.Complex[n];
228      cmatrixluinverserec(ref a, 0, n, ref work, ref info, ref rep);
229
230      //
231      // apply permutations
232      //
233      for (i = 0; i <= n - 1; i++) {
234        for (j = n - 2; j >= 0; j--) {
235          k = pivots[j];
236          v = a[i, j];
237          a[i, j] = a[i, k];
238          a[i, k] = v;
239        }
240      }
241    }
242
243
244    /*************************************************************************
245    Inversion of a general matrix.
246
247    Input parameters:
248        A   -   matrix, array[0..N-1,0..N-1].
249        N   -   size of A.
250
251    Output parameters:
252        Info    -   return code, same as in RMatrixLUInverse
253        Rep     -   solver report, same as in RMatrixLUInverse
254        A       -   inverse of matrix A, same as in RMatrixLUInverse
255
256      -- ALGLIB --
257         Copyright 2005 by Bochkanov Sergey
258    *************************************************************************/
259    public static void cmatrixinverse(ref AP.Complex[,] a,
260        int n,
261        ref int info,
262        ref matinvreport rep) {
263      int[] pivots = new int[0];
264
265      trfac.cmatrixlu(ref a, n, n, ref pivots);
266      cmatrixluinverse(ref a, ref pivots, n, ref info, ref rep);
267    }
268
269
270    /*************************************************************************
271    Inversion of a symmetric positive definite matrix which is given
272    by Cholesky decomposition.
273
274    Input parameters:
275        A       -   Cholesky decomposition of the matrix to be inverted:
276                    A=U’*U or A = L*L'.
277                    Output of  SPDMatrixCholesky subroutine.
278        N       -   size of matrix A.
279        IsUpper –   storage format.
280                    If IsUpper = True, then matrix A is given as A = U'*U
281                    (matrix contains upper triangle).
282                    Similarly, if IsUpper = False, then A = L*L'.
283
284    Output parameters:
285        Info    -   return code, same as in RMatrixLUInverse
286        Rep     -   solver report, same as in RMatrixLUInverse
287        A       -   inverse of matrix A, same as in RMatrixLUInverse
288
289      -- ALGLIB routine --
290         10.02.2010
291         Bochkanov Sergey
292    *************************************************************************/
293    public static void spdmatrixcholeskyinverse(ref double[,] a,
294        int n,
295        bool isupper,
296        ref int info,
297        ref matinvreport rep) {
298      int i = 0;
299      int j = 0;
300      int k = 0;
301      double v = 0;
302      double ajj = 0;
303      double aii = 0;
304      double[] tmp = new double[0];
305      int info2 = 0;
306      matinvreport rep2 = new matinvreport();
307
308      if (n < 1) {
309        info = -1;
310        return;
311      }
312      info = 1;
313
314      //
315      // calculate condition numbers
316      //
317      rep.r1 = rcond.spdmatrixcholeskyrcond(ref a, n, isupper);
318      rep.rinf = rep.r1;
319      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
320        if (isupper) {
321          for (i = 0; i <= n - 1; i++) {
322            for (j = i; j <= n - 1; j++) {
323              a[i, j] = 0;
324            }
325          }
326        } else {
327          for (i = 0; i <= n - 1; i++) {
328            for (j = 0; j <= i; j++) {
329              a[i, j] = 0;
330            }
331          }
332        }
333        rep.r1 = 0;
334        rep.rinf = 0;
335        info = -3;
336        return;
337      }
338
339      //
340      // Inverse
341      //
342      tmp = new double[n];
343      spdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
344    }
345
346
347    /*************************************************************************
348    Inversion of a symmetric positive definite matrix.
349
350    Given an upper or lower triangle of a symmetric positive definite matrix,
351    the algorithm generates matrix A^-1 and saves the upper or lower triangle
352    depending on the input.
353
354    Input parameters:
355        A       -   matrix to be inverted (upper or lower triangle).
356                    Array with elements [0..N-1,0..N-1].
357        N       -   size of matrix A.
358        IsUpper -   storage format.
359                    If IsUpper = True, then the upper triangle of matrix A is
360                    given, otherwise the lower triangle is given.
361
362    Output parameters:
363        Info    -   return code, same as in RMatrixLUInverse
364        Rep     -   solver report, same as in RMatrixLUInverse
365        A       -   inverse of matrix A, same as in RMatrixLUInverse
366
367      -- ALGLIB routine --
368         10.02.2010
369         Bochkanov Sergey
370    *************************************************************************/
371    public static void spdmatrixinverse(ref double[,] a,
372        int n,
373        bool isupper,
374        ref int info,
375        ref matinvreport rep) {
376      if (n < 1) {
377        info = -1;
378        return;
379      }
380      info = 1;
381      if (trfac.spdmatrixcholesky(ref a, n, isupper)) {
382        spdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
383      } else {
384        info = -3;
385      }
386    }
387
388
389    /*************************************************************************
390    Inversion of a Hermitian positive definite matrix which is given
391    by Cholesky decomposition.
392
393    Input parameters:
394        A       -   Cholesky decomposition of the matrix to be inverted:
395                    A=U’*U or A = L*L'.
396                    Output of  HPDMatrixCholesky subroutine.
397        N       -   size of matrix A.
398        IsUpper –   storage format.
399                    If IsUpper = True, then matrix A is given as A = U'*U
400                    (matrix contains upper triangle).
401                    Similarly, if IsUpper = False, then A = L*L'.
402
403    Output parameters:
404        Info    -   return code, same as in RMatrixLUInverse
405        Rep     -   solver report, same as in RMatrixLUInverse
406        A       -   inverse of matrix A, same as in RMatrixLUInverse
407
408      -- ALGLIB routine --
409         10.02.2010
410         Bochkanov Sergey
411    *************************************************************************/
412    public static void hpdmatrixcholeskyinverse(ref AP.Complex[,] a,
413        int n,
414        bool isupper,
415        ref int info,
416        ref matinvreport rep) {
417      int i = 0;
418      int j = 0;
419      int info2 = 0;
420      matinvreport rep2 = new matinvreport();
421      AP.Complex[] tmp = new AP.Complex[0];
422      AP.Complex v = 0;
423
424      if (n < 1) {
425        info = -1;
426        return;
427      }
428      info = 1;
429
430      //
431      // calculate condition numbers
432      //
433      rep.r1 = rcond.hpdmatrixcholeskyrcond(ref a, n, isupper);
434      rep.rinf = rep.r1;
435      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
436        if (isupper) {
437          for (i = 0; i <= n - 1; i++) {
438            for (j = i; j <= n - 1; j++) {
439              a[i, j] = 0;
440            }
441          }
442        } else {
443          for (i = 0; i <= n - 1; i++) {
444            for (j = 0; j <= i; j++) {
445              a[i, j] = 0;
446            }
447          }
448        }
449        rep.r1 = 0;
450        rep.rinf = 0;
451        info = -3;
452        return;
453      }
454
455      //
456      // Inverse
457      //
458      tmp = new AP.Complex[n];
459      hpdmatrixcholeskyinverserec(ref a, 0, n, isupper, ref tmp);
460    }
461
462
463    /*************************************************************************
464    Inversion of a Hermitian positive definite matrix.
465
466    Given an upper or lower triangle of a Hermitian positive definite matrix,
467    the algorithm generates matrix A^-1 and saves the upper or lower triangle
468    depending on the input.
469
470    Input parameters:
471        A       -   matrix to be inverted (upper or lower triangle).
472                    Array with elements [0..N-1,0..N-1].
473        N       -   size of matrix A.
474        IsUpper -   storage format.
475                    If IsUpper = True, then the upper triangle of matrix A is
476                    given, otherwise the lower triangle is given.
477
478    Output parameters:
479        Info    -   return code, same as in RMatrixLUInverse
480        Rep     -   solver report, same as in RMatrixLUInverse
481        A       -   inverse of matrix A, same as in RMatrixLUInverse
482
483      -- ALGLIB routine --
484         10.02.2010
485         Bochkanov Sergey
486    *************************************************************************/
487    public static void hpdmatrixinverse(ref AP.Complex[,] a,
488        int n,
489        bool isupper,
490        ref int info,
491        ref matinvreport rep) {
492      if (n < 1) {
493        info = -1;
494        return;
495      }
496      info = 1;
497      if (trfac.hpdmatrixcholesky(ref a, n, isupper)) {
498        hpdmatrixcholeskyinverse(ref a, n, isupper, ref info, ref rep);
499      } else {
500        info = -3;
501      }
502    }
503
504
505    /*************************************************************************
506    Triangular matrix inverse (real)
507
508    The subroutine inverts the following types of matrices:
509        * upper triangular
510        * upper triangular with unit diagonal
511        * lower triangular
512        * lower triangular with unit diagonal
513
514    In case of an upper (lower) triangular matrix,  the  inverse  matrix  will
515    also be upper (lower) triangular, and after the end of the algorithm,  the
516    inverse matrix replaces the source matrix. The elements  below (above) the
517    main diagonal are not changed by the algorithm.
518
519    If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
520    diagonal, and the diagonal elements are not passed to the algorithm.
521
522    Input parameters:
523        A       -   matrix, array[0..N-1, 0..N-1].
524        N       -   size of A.
525        IsUpper -   True, if the matrix is upper triangular.
526        IsUnit  -   True, if the matrix has a unit diagonal.
527
528    Output parameters:
529        Info    -   same as for RMatrixLUInverse
530        Rep     -   same as for RMatrixLUInverse
531        A       -   same as for RMatrixLUInverse.
532
533      -- ALGLIB --
534         Copyright 05.02.2010 by Bochkanov Sergey
535    *************************************************************************/
536    public static void rmatrixtrinverse(ref double[,] a,
537        int n,
538        bool isupper,
539        bool isunit,
540        ref int info,
541        ref matinvreport rep) {
542      int i = 0;
543      int j = 0;
544      double[] tmp = new double[0];
545
546      if (n < 1) {
547        info = -1;
548        return;
549      }
550      info = 1;
551
552      //
553      // calculate condition numbers
554      //
555      rep.r1 = rcond.rmatrixtrrcond1(ref a, n, isupper, isunit);
556      rep.rinf = rcond.rmatrixtrrcondinf(ref a, n, isupper, isunit);
557      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
558        for (i = 0; i <= n - 1; i++) {
559          for (j = 0; j <= n - 1; j++) {
560            a[i, j] = 0;
561          }
562        }
563        rep.r1 = 0;
564        rep.rinf = 0;
565        info = -3;
566        return;
567      }
568
569      //
570      // Invert
571      //
572      tmp = new double[n];
573      rmatrixtrinverserec(ref a, 0, n, isupper, isunit, ref tmp, ref info, ref rep);
574    }
575
576
577    /*************************************************************************
578    Triangular matrix inverse (complex)
579
580    The subroutine inverts the following types of matrices:
581        * upper triangular
582        * upper triangular with unit diagonal
583        * lower triangular
584        * lower triangular with unit diagonal
585
586    In case of an upper (lower) triangular matrix,  the  inverse  matrix  will
587    also be upper (lower) triangular, and after the end of the algorithm,  the
588    inverse matrix replaces the source matrix. The elements  below (above) the
589    main diagonal are not changed by the algorithm.
590
591    If  the matrix  has a unit diagonal, the inverse matrix also  has  a  unit
592    diagonal, and the diagonal elements are not passed to the algorithm.
593
594    Input parameters:
595        A       -   matrix, array[0..N-1, 0..N-1].
596        N       -   size of A.
597        IsUpper -   True, if the matrix is upper triangular.
598        IsUnit  -   True, if the matrix has a unit diagonal.
599
600    Output parameters:
601        Info    -   same as for RMatrixLUInverse
602        Rep     -   same as for RMatrixLUInverse
603        A       -   same as for RMatrixLUInverse.
604
605      -- ALGLIB --
606         Copyright 05.02.2010 by Bochkanov Sergey
607    *************************************************************************/
608    public static void cmatrixtrinverse(ref AP.Complex[,] a,
609        int n,
610        bool isupper,
611        bool isunit,
612        ref int info,
613        ref matinvreport rep) {
614      int i = 0;
615      int j = 0;
616      AP.Complex[] tmp = new AP.Complex[0];
617
618      if (n < 1) {
619        info = -1;
620        return;
621      }
622      info = 1;
623
624      //
625      // calculate condition numbers
626      //
627      rep.r1 = rcond.cmatrixtrrcond1(ref a, n, isupper, isunit);
628      rep.rinf = rcond.cmatrixtrrcondinf(ref a, n, isupper, isunit);
629      if ((double)(rep.r1) < (double)(rcond.rcondthreshold()) | (double)(rep.rinf) < (double)(rcond.rcondthreshold())) {
630        for (i = 0; i <= n - 1; i++) {
631          for (j = 0; j <= n - 1; j++) {
632            a[i, j] = 0;
633          }
634        }
635        rep.r1 = 0;
636        rep.rinf = 0;
637        info = -3;
638        return;
639      }
640
641      //
642      // Invert
643      //
644      tmp = new AP.Complex[n];
645      cmatrixtrinverserec(ref a, 0, n, isupper, isunit, ref tmp, ref info, ref rep);
646    }
647
648
649    /*************************************************************************
650    Triangular matrix inversion, recursive subroutine
651
652      -- ALGLIB --
653         05.02.2010, Bochkanov Sergey.
654         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
655         Courant Institute, Argonne National Lab, and Rice University
656         February 29, 1992.
657    *************************************************************************/
658    private static void rmatrixtrinverserec(ref double[,] a,
659        int offs,
660        int n,
661        bool isupper,
662        bool isunit,
663        ref double[] tmp,
664        ref int info,
665        ref matinvreport rep) {
666      int n1 = 0;
667      int n2 = 0;
668      int i = 0;
669      int j = 0;
670      double v = 0;
671      double ajj = 0;
672      int i_ = 0;
673      int i1_ = 0;
674
675      if (n < 1) {
676        info = -1;
677        return;
678      }
679
680      //
681      // Base case
682      //
683      if (n <= ablas.ablasblocksize(ref a)) {
684        if (isupper) {
685
686          //
687          // Compute inverse of upper triangular matrix.
688          //
689          for (j = 0; j <= n - 1; j++) {
690            if (!isunit) {
691              if ((double)(a[offs + j, offs + j]) == (double)(0)) {
692                info = -3;
693                return;
694              }
695              a[offs + j, offs + j] = 1 / a[offs + j, offs + j];
696              ajj = -a[offs + j, offs + j];
697            } else {
698              ajj = -1;
699            }
700
701            //
702            // Compute elements 1:j-1 of j-th column.
703            //
704            if (j > 0) {
705              i1_ = (offs + 0) - (0);
706              for (i_ = 0; i_ <= j - 1; i_++) {
707                tmp[i_] = a[i_ + i1_, offs + j];
708              }
709              for (i = 0; i <= j - 1; i++) {
710                if (i < j - 1) {
711                  i1_ = (i + 1) - (offs + i + 1);
712                  v = 0.0;
713                  for (i_ = offs + i + 1; i_ <= offs + j - 1; i_++) {
714                    v += a[offs + i, i_] * tmp[i_ + i1_];
715                  }
716                } else {
717                  v = 0;
718                }
719                if (!isunit) {
720                  a[offs + i, offs + j] = v + a[offs + i, offs + i] * tmp[i];
721                } else {
722                  a[offs + i, offs + j] = v + tmp[i];
723                }
724              }
725              for (i_ = offs + 0; i_ <= offs + j - 1; i_++) {
726                a[i_, offs + j] = ajj * a[i_, offs + j];
727              }
728            }
729          }
730        } else {
731
732          //
733          // Compute inverse of lower triangular matrix.
734          //
735          for (j = n - 1; j >= 0; j--) {
736            if (!isunit) {
737              if ((double)(a[offs + j, offs + j]) == (double)(0)) {
738                info = -3;
739                return;
740              }
741              a[offs + j, offs + j] = 1 / a[offs + j, offs + j];
742              ajj = -a[offs + j, offs + j];
743            } else {
744              ajj = -1;
745            }
746            if (j < n - 1) {
747
748              //
749              // Compute elements j+1:n of j-th column.
750              //
751              i1_ = (offs + j + 1) - (j + 1);
752              for (i_ = j + 1; i_ <= n - 1; i_++) {
753                tmp[i_] = a[i_ + i1_, offs + j];
754              }
755              for (i = j + 1; i <= n - 1; i++) {
756                if (i > j + 1) {
757                  i1_ = (j + 1) - (offs + j + 1);
758                  v = 0.0;
759                  for (i_ = offs + j + 1; i_ <= offs + i - 1; i_++) {
760                    v += a[offs + i, i_] * tmp[i_ + i1_];
761                  }
762                } else {
763                  v = 0;
764                }
765                if (!isunit) {
766                  a[offs + i, offs + j] = v + a[offs + i, offs + i] * tmp[i];
767                } else {
768                  a[offs + i, offs + j] = v + tmp[i];
769                }
770              }
771              for (i_ = offs + j + 1; i_ <= offs + n - 1; i_++) {
772                a[i_, offs + j] = ajj * a[i_, offs + j];
773              }
774            }
775          }
776        }
777        return;
778      }
779
780      //
781      // Recursive case
782      //
783      ablas.ablassplitlength(ref a, n, ref n1, ref n2);
784      if (n2 > 0) {
785        if (isupper) {
786          for (i = 0; i <= n1 - 1; i++) {
787            for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
788              a[offs + i, i_] = -1 * a[offs + i, i_];
789            }
790          }
791          ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, isunit, 0, ref a, offs, offs + n1);
792          ablas.rmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, isunit, 0, ref a, offs, offs + n1);
793        } else {
794          for (i = 0; i <= n2 - 1; i++) {
795            for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
796              a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
797            }
798          }
799          ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, isunit, 0, ref a, offs + n1, offs);
800          ablas.rmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, isunit, 0, ref a, offs + n1, offs);
801        }
802        rmatrixtrinverserec(ref a, offs + n1, n2, isupper, isunit, ref tmp, ref info, ref rep);
803      }
804      rmatrixtrinverserec(ref a, offs, n1, isupper, isunit, ref tmp, ref info, ref rep);
805    }
806
807
808    /*************************************************************************
809    Triangular matrix inversion, recursive subroutine
810
811      -- ALGLIB --
812         05.02.2010, Bochkanov Sergey.
813         Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
814         Courant Institute, Argonne National Lab, and Rice University
815         February 29, 1992.
816    *************************************************************************/
817    private static void cmatrixtrinverserec(ref AP.Complex[,] a,
818        int offs,
819        int n,
820        bool isupper,
821        bool isunit,
822        ref AP.Complex[] tmp,
823        ref int info,
824        ref matinvreport rep) {
825      int n1 = 0;
826      int n2 = 0;
827      int i = 0;
828      int j = 0;
829      AP.Complex v = 0;
830      AP.Complex ajj = 0;
831      int i_ = 0;
832      int i1_ = 0;
833
834      if (n < 1) {
835        info = -1;
836        return;
837      }
838
839      //
840      // Base case
841      //
842      if (n <= ablas.ablascomplexblocksize(ref a)) {
843        if (isupper) {
844
845          //
846          // Compute inverse of upper triangular matrix.
847          //
848          for (j = 0; j <= n - 1; j++) {
849            if (!isunit) {
850              if (a[offs + j, offs + j] == 0) {
851                info = -3;
852                return;
853              }
854              a[offs + j, offs + j] = 1 / a[offs + j, offs + j];
855              ajj = -a[offs + j, offs + j];
856            } else {
857              ajj = -1;
858            }
859
860            //
861            // Compute elements 1:j-1 of j-th column.
862            //
863            if (j > 0) {
864              i1_ = (offs + 0) - (0);
865              for (i_ = 0; i_ <= j - 1; i_++) {
866                tmp[i_] = a[i_ + i1_, offs + j];
867              }
868              for (i = 0; i <= j - 1; i++) {
869                if (i < j - 1) {
870                  i1_ = (i + 1) - (offs + i + 1);
871                  v = 0.0;
872                  for (i_ = offs + i + 1; i_ <= offs + j - 1; i_++) {
873                    v += a[offs + i, i_] * tmp[i_ + i1_];
874                  }
875                } else {
876                  v = 0;
877                }
878                if (!isunit) {
879                  a[offs + i, offs + j] = v + a[offs + i, offs + i] * tmp[i];
880                } else {
881                  a[offs + i, offs + j] = v + tmp[i];
882                }
883              }
884              for (i_ = offs + 0; i_ <= offs + j - 1; i_++) {
885                a[i_, offs + j] = ajj * a[i_, offs + j];
886              }
887            }
888          }
889        } else {
890
891          //
892          // Compute inverse of lower triangular matrix.
893          //
894          for (j = n - 1; j >= 0; j--) {
895            if (!isunit) {
896              if (a[offs + j, offs + j] == 0) {
897                info = -3;
898                return;
899              }
900              a[offs + j, offs + j] = 1 / a[offs + j, offs + j];
901              ajj = -a[offs + j, offs + j];
902            } else {
903              ajj = -1;
904            }
905            if (j < n - 1) {
906
907              //
908              // Compute elements j+1:n of j-th column.
909              //
910              i1_ = (offs + j + 1) - (j + 1);
911              for (i_ = j + 1; i_ <= n - 1; i_++) {
912                tmp[i_] = a[i_ + i1_, offs + j];
913              }
914              for (i = j + 1; i <= n - 1; i++) {
915                if (i > j + 1) {
916                  i1_ = (j + 1) - (offs + j + 1);
917                  v = 0.0;
918                  for (i_ = offs + j + 1; i_ <= offs + i - 1; i_++) {
919                    v += a[offs + i, i_] * tmp[i_ + i1_];
920                  }
921                } else {
922                  v = 0;
923                }
924                if (!isunit) {
925                  a[offs + i, offs + j] = v + a[offs + i, offs + i] * tmp[i];
926                } else {
927                  a[offs + i, offs + j] = v + tmp[i];
928                }
929              }
930              for (i_ = offs + j + 1; i_ <= offs + n - 1; i_++) {
931                a[i_, offs + j] = ajj * a[i_, offs + j];
932              }
933            }
934          }
935        }
936        return;
937      }
938
939      //
940      // Recursive case
941      //
942      ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
943      if (n2 > 0) {
944        if (isupper) {
945          for (i = 0; i <= n1 - 1; i++) {
946            for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
947              a[offs + i, i_] = -1 * a[offs + i, i_];
948            }
949          }
950          ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, isunit, 0, ref a, offs, offs + n1);
951          ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, isunit, 0, ref a, offs, offs + n1);
952        } else {
953          for (i = 0; i <= n2 - 1; i++) {
954            for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
955              a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
956            }
957          }
958          ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, isunit, 0, ref a, offs + n1, offs);
959          ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, isunit, 0, ref a, offs + n1, offs);
960        }
961        cmatrixtrinverserec(ref a, offs + n1, n2, isupper, isunit, ref tmp, ref info, ref rep);
962      }
963      cmatrixtrinverserec(ref a, offs, n1, isupper, isunit, ref tmp, ref info, ref rep);
964    }
965
966
967    private static void rmatrixluinverserec(ref double[,] a,
968        int offs,
969        int n,
970        ref double[] work,
971        ref int info,
972        ref matinvreport rep) {
973      int i = 0;
974      int iws = 0;
975      int j = 0;
976      int jb = 0;
977      int jj = 0;
978      int jp = 0;
979      int k = 0;
980      double v = 0;
981      int n1 = 0;
982      int n2 = 0;
983      int i_ = 0;
984      int i1_ = 0;
985
986      if (n < 1) {
987        info = -1;
988        return;
989      }
990
991      //
992      // Base case
993      //
994      if (n <= ablas.ablasblocksize(ref a)) {
995
996        //
997        // Form inv(U)
998        //
999        rmatrixtrinverserec(ref a, offs, n, true, false, ref work, ref info, ref rep);
1000        if (info <= 0) {
1001          return;
1002        }
1003
1004        //
1005        // Solve the equation inv(A)*L = inv(U) for inv(A).
1006        //
1007        for (j = n - 1; j >= 0; j--) {
1008
1009          //
1010          // Copy current column of L to WORK and replace with zeros.
1011          //
1012          for (i = j + 1; i <= n - 1; i++) {
1013            work[i] = a[offs + i, offs + j];
1014            a[offs + i, offs + j] = 0;
1015          }
1016
1017          //
1018          // Compute current column of inv(A).
1019          //
1020          if (j < n - 1) {
1021            for (i = 0; i <= n - 1; i++) {
1022              i1_ = (j + 1) - (offs + j + 1);
1023              v = 0.0;
1024              for (i_ = offs + j + 1; i_ <= offs + n - 1; i_++) {
1025                v += a[offs + i, i_] * work[i_ + i1_];
1026              }
1027              a[offs + i, offs + j] = a[offs + i, offs + j] - v;
1028            }
1029          }
1030        }
1031        return;
1032      }
1033
1034      //
1035      // Recursive code:
1036      //
1037      //         ( L1      )   ( U1  U12 )
1038      // A    =  (         ) * (         )
1039      //         ( L12  L2 )   (     U2  )
1040      //
1041      //         ( W   X )
1042      // A^-1 =  (       )
1043      //         ( Y   Z )
1044      //
1045      ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1046      System.Diagnostics.Debug.Assert(n2 > 0, "LUInverseRec: internal error!");
1047
1048      //
1049      // X := inv(U1)*U12*inv(U2)
1050      //
1051      ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, true, false, 0, ref a, offs, offs + n1);
1052      ablas.rmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, true, false, 0, ref a, offs, offs + n1);
1053
1054      //
1055      // Y := inv(L2)*L12*inv(L1)
1056      //
1057      ablas.rmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, false, true, 0, ref a, offs + n1, offs);
1058      ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, false, true, 0, ref a, offs + n1, offs);
1059
1060      //
1061      // W := inv(L1*U1)+X*Y
1062      //
1063      rmatrixluinverserec(ref a, offs, n1, ref work, ref info, ref rep);
1064      if (info <= 0) {
1065        return;
1066      }
1067      ablas.rmatrixgemm(n1, n1, n2, 1.0, ref a, offs, offs + n1, 0, ref a, offs + n1, offs, 0, 1.0, ref a, offs, offs);
1068
1069      //
1070      // X := -X*inv(L2)
1071      // Y := -inv(U2)*Y
1072      //
1073      ablas.rmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, false, true, 0, ref a, offs, offs + n1);
1074      for (i = 0; i <= n1 - 1; i++) {
1075        for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
1076          a[offs + i, i_] = -1 * a[offs + i, i_];
1077        }
1078      }
1079      ablas.rmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, true, false, 0, ref a, offs + n1, offs);
1080      for (i = 0; i <= n2 - 1; i++) {
1081        for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
1082          a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
1083        }
1084      }
1085
1086      //
1087      // Z := inv(L2*U2)
1088      //
1089      rmatrixluinverserec(ref a, offs + n1, n2, ref work, ref info, ref rep);
1090    }
1091
1092
1093    private static void cmatrixluinverserec(ref AP.Complex[,] a,
1094        int offs,
1095        int n,
1096        ref AP.Complex[] work,
1097        ref int info,
1098        ref matinvreport rep) {
1099      int i = 0;
1100      int iws = 0;
1101      int j = 0;
1102      int jb = 0;
1103      int jj = 0;
1104      int jp = 0;
1105      int k = 0;
1106      AP.Complex v = 0;
1107      int n1 = 0;
1108      int n2 = 0;
1109      int i_ = 0;
1110      int i1_ = 0;
1111
1112      if (n < 1) {
1113        info = -1;
1114        return;
1115      }
1116
1117      //
1118      // Base case
1119      //
1120      if (n <= ablas.ablascomplexblocksize(ref a)) {
1121
1122        //
1123        // Form inv(U)
1124        //
1125        cmatrixtrinverserec(ref a, offs, n, true, false, ref work, ref info, ref rep);
1126        if (info <= 0) {
1127          return;
1128        }
1129
1130        //
1131        // Solve the equation inv(A)*L = inv(U) for inv(A).
1132        //
1133        for (j = n - 1; j >= 0; j--) {
1134
1135          //
1136          // Copy current column of L to WORK and replace with zeros.
1137          //
1138          for (i = j + 1; i <= n - 1; i++) {
1139            work[i] = a[offs + i, offs + j];
1140            a[offs + i, offs + j] = 0;
1141          }
1142
1143          //
1144          // Compute current column of inv(A).
1145          //
1146          if (j < n - 1) {
1147            for (i = 0; i <= n - 1; i++) {
1148              i1_ = (j + 1) - (offs + j + 1);
1149              v = 0.0;
1150              for (i_ = offs + j + 1; i_ <= offs + n - 1; i_++) {
1151                v += a[offs + i, i_] * work[i_ + i1_];
1152              }
1153              a[offs + i, offs + j] = a[offs + i, offs + j] - v;
1154            }
1155          }
1156        }
1157        return;
1158      }
1159
1160      //
1161      // Recursive code:
1162      //
1163      //         ( L1      )   ( U1  U12 )
1164      // A    =  (         ) * (         )
1165      //         ( L12  L2 )   (     U2  )
1166      //
1167      //         ( W   X )
1168      // A^-1 =  (       )
1169      //         ( Y   Z )
1170      //
1171      ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1172      System.Diagnostics.Debug.Assert(n2 > 0, "LUInverseRec: internal error!");
1173
1174      //
1175      // X := inv(U1)*U12*inv(U2)
1176      //
1177      ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, true, false, 0, ref a, offs, offs + n1);
1178      ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, true, false, 0, ref a, offs, offs + n1);
1179
1180      //
1181      // Y := inv(L2)*L12*inv(L1)
1182      //
1183      ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, false, true, 0, ref a, offs + n1, offs);
1184      ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, false, true, 0, ref a, offs + n1, offs);
1185
1186      //
1187      // W := inv(L1*U1)+X*Y
1188      //
1189      cmatrixluinverserec(ref a, offs, n1, ref work, ref info, ref rep);
1190      if (info <= 0) {
1191        return;
1192      }
1193      ablas.cmatrixgemm(n1, n1, n2, 1.0, ref a, offs, offs + n1, 0, ref a, offs + n1, offs, 0, 1.0, ref a, offs, offs);
1194
1195      //
1196      // X := -X*inv(L2)
1197      // Y := -inv(U2)*Y
1198      //
1199      ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, false, true, 0, ref a, offs, offs + n1);
1200      for (i = 0; i <= n1 - 1; i++) {
1201        for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
1202          a[offs + i, i_] = -1 * a[offs + i, i_];
1203        }
1204      }
1205      ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, true, false, 0, ref a, offs + n1, offs);
1206      for (i = 0; i <= n2 - 1; i++) {
1207        for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
1208          a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
1209        }
1210      }
1211
1212      //
1213      // Z := inv(L2*U2)
1214      //
1215      cmatrixluinverserec(ref a, offs + n1, n2, ref work, ref info, ref rep);
1216    }
1217
1218
1219    /*************************************************************************
1220    Recursive subroutine for SPD inversion.
1221
1222      -- ALGLIB routine --
1223         10.02.2010
1224         Bochkanov Sergey
1225    *************************************************************************/
1226    private static void spdmatrixcholeskyinverserec(ref double[,] a,
1227        int offs,
1228        int n,
1229        bool isupper,
1230        ref double[] tmp) {
1231      int i = 0;
1232      int j = 0;
1233      double v = 0;
1234      int n1 = 0;
1235      int n2 = 0;
1236      int info2 = 0;
1237      matinvreport rep2 = new matinvreport();
1238      int i_ = 0;
1239      int i1_ = 0;
1240
1241      if (n < 1) {
1242        return;
1243      }
1244
1245      //
1246      // Base case
1247      //
1248      if (n <= ablas.ablasblocksize(ref a)) {
1249        rmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
1250        if (isupper) {
1251
1252          //
1253          // Compute the product U * U'.
1254          // NOTE: we never assume that diagonal of U is real
1255          //
1256          for (i = 0; i <= n - 1; i++) {
1257            if (i == 0) {
1258
1259              //
1260              // 1x1 matrix
1261              //
1262              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i]);
1263            } else {
1264
1265              //
1266              // (I+1)x(I+1) matrix,
1267              //
1268              // ( A11  A12 )   ( A11^H        )   ( A11*A11^H+A12*A12^H  A12*A22^H )
1269              // (          ) * (              ) = (                                )
1270              // (      A22 )   ( A12^H  A22^H )   ( A22*A12^H            A22*A22^H )
1271              //
1272              // A11 is IxI, A22 is 1x1.
1273              //
1274              i1_ = (offs) - (0);
1275              for (i_ = 0; i_ <= i - 1; i_++) {
1276                tmp[i_] = a[i_ + i1_, offs + i];
1277              }
1278              for (j = 0; j <= i - 1; j++) {
1279                v = a[offs + j, offs + i];
1280                i1_ = (j) - (offs + j);
1281                for (i_ = offs + j; i_ <= offs + i - 1; i_++) {
1282                  a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
1283                }
1284              }
1285              v = a[offs + i, offs + i];
1286              for (i_ = offs; i_ <= offs + i - 1; i_++) {
1287                a[i_, offs + i] = v * a[i_, offs + i];
1288              }
1289              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i]);
1290            }
1291          }
1292        } else {
1293
1294          //
1295          // Compute the product L' * L
1296          // NOTE: we never assume that diagonal of L is real
1297          //
1298          for (i = 0; i <= n - 1; i++) {
1299            if (i == 0) {
1300
1301              //
1302              // 1x1 matrix
1303              //
1304              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i]);
1305            } else {
1306
1307              //
1308              // (I+1)x(I+1) matrix,
1309              //
1310              // ( A11^H  A21^H )   ( A11      )   ( A11^H*A11+A21^H*A21  A21^H*A22 )
1311              // (              ) * (          ) = (                                )
1312              // (        A22^H )   ( A21  A22 )   ( A22^H*A21            A22^H*A22 )
1313              //
1314              // A11 is IxI, A22 is 1x1.
1315              //
1316              i1_ = (offs) - (0);
1317              for (i_ = 0; i_ <= i - 1; i_++) {
1318                tmp[i_] = a[offs + i, i_ + i1_];
1319              }
1320              for (j = 0; j <= i - 1; j++) {
1321                v = a[offs + i, offs + j];
1322                i1_ = (0) - (offs);
1323                for (i_ = offs; i_ <= offs + j; i_++) {
1324                  a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
1325                }
1326              }
1327              v = a[offs + i, offs + i];
1328              for (i_ = offs; i_ <= offs + i - 1; i_++) {
1329                a[offs + i, i_] = v * a[offs + i, i_];
1330              }
1331              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i]);
1332            }
1333          }
1334        }
1335        return;
1336      }
1337
1338      //
1339      // Recursive code: triangular factor inversion merged with
1340      // UU' or L'L multiplication
1341      //
1342      ablas.ablassplitlength(ref a, n, ref n1, ref n2);
1343
1344      //
1345      // form off-diagonal block of trangular inverse
1346      //
1347      if (isupper) {
1348        for (i = 0; i <= n1 - 1; i++) {
1349          for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
1350            a[offs + i, i_] = -1 * a[offs + i, i_];
1351          }
1352        }
1353        ablas.rmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 0, ref a, offs, offs + n1);
1354        ablas.rmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, false, 0, ref a, offs, offs + n1);
1355      } else {
1356        for (i = 0; i <= n2 - 1; i++) {
1357          for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
1358            a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
1359          }
1360        }
1361        ablas.rmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 0, ref a, offs + n1, offs);
1362        ablas.rmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, false, 0, ref a, offs + n1, offs);
1363      }
1364
1365      //
1366      // invert first diagonal block
1367      //
1368      spdmatrixcholeskyinverserec(ref a, offs, n1, isupper, ref tmp);
1369
1370      //
1371      // update first diagonal block with off-diagonal block,
1372      // update off-diagonal block
1373      //
1374      if (isupper) {
1375        ablas.rmatrixsyrk(n1, n2, 1.0, ref a, offs, offs + n1, 0, 1.0, ref a, offs, offs, isupper);
1376        ablas.rmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, false, 1, ref a, offs, offs + n1);
1377      } else {
1378        ablas.rmatrixsyrk(n1, n2, 1.0, ref a, offs + n1, offs, 1, 1.0, ref a, offs, offs, isupper);
1379        ablas.rmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, false, 1, ref a, offs + n1, offs);
1380      }
1381
1382      //
1383      // invert second diagonal block
1384      //
1385      spdmatrixcholeskyinverserec(ref a, offs + n1, n2, isupper, ref tmp);
1386    }
1387
1388
1389    /*************************************************************************
1390    Recursive subroutine for HPD inversion.
1391
1392      -- ALGLIB routine --
1393         10.02.2010
1394         Bochkanov Sergey
1395    *************************************************************************/
1396    private static void hpdmatrixcholeskyinverserec(ref AP.Complex[,] a,
1397        int offs,
1398        int n,
1399        bool isupper,
1400        ref AP.Complex[] tmp) {
1401      int i = 0;
1402      int j = 0;
1403      AP.Complex v = 0;
1404      int n1 = 0;
1405      int n2 = 0;
1406      int info2 = 0;
1407      matinvreport rep2 = new matinvreport();
1408      int i_ = 0;
1409      int i1_ = 0;
1410
1411      if (n < 1) {
1412        return;
1413      }
1414
1415      //
1416      // Base case
1417      //
1418      if (n <= ablas.ablascomplexblocksize(ref a)) {
1419        cmatrixtrinverserec(ref a, offs, n, isupper, false, ref tmp, ref info2, ref rep2);
1420        if (isupper) {
1421
1422          //
1423          // Compute the product U * U'.
1424          // NOTE: we never assume that diagonal of U is real
1425          //
1426          for (i = 0; i <= n - 1; i++) {
1427            if (i == 0) {
1428
1429              //
1430              // 1x1 matrix
1431              //
1432              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
1433            } else {
1434
1435              //
1436              // (I+1)x(I+1) matrix,
1437              //
1438              // ( A11  A12 )   ( A11^H        )   ( A11*A11^H+A12*A12^H  A12*A22^H )
1439              // (          ) * (              ) = (                                )
1440              // (      A22 )   ( A12^H  A22^H )   ( A22*A12^H            A22*A22^H )
1441              //
1442              // A11 is IxI, A22 is 1x1.
1443              //
1444              i1_ = (offs) - (0);
1445              for (i_ = 0; i_ <= i - 1; i_++) {
1446                tmp[i_] = AP.Math.Conj(a[i_ + i1_, offs + i]);
1447              }
1448              for (j = 0; j <= i - 1; j++) {
1449                v = a[offs + j, offs + i];
1450                i1_ = (j) - (offs + j);
1451                for (i_ = offs + j; i_ <= offs + i - 1; i_++) {
1452                  a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
1453                }
1454              }
1455              v = AP.Math.Conj(a[offs + i, offs + i]);
1456              for (i_ = offs; i_ <= offs + i - 1; i_++) {
1457                a[i_, offs + i] = v * a[i_, offs + i];
1458              }
1459              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
1460            }
1461          }
1462        } else {
1463
1464          //
1465          // Compute the product L' * L
1466          // NOTE: we never assume that diagonal of L is real
1467          //
1468          for (i = 0; i <= n - 1; i++) {
1469            if (i == 0) {
1470
1471              //
1472              // 1x1 matrix
1473              //
1474              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
1475            } else {
1476
1477              //
1478              // (I+1)x(I+1) matrix,
1479              //
1480              // ( A11^H  A21^H )   ( A11      )   ( A11^H*A11+A21^H*A21  A21^H*A22 )
1481              // (              ) * (          ) = (                                )
1482              // (        A22^H )   ( A21  A22 )   ( A22^H*A21            A22^H*A22 )
1483              //
1484              // A11 is IxI, A22 is 1x1.
1485              //
1486              i1_ = (offs) - (0);
1487              for (i_ = 0; i_ <= i - 1; i_++) {
1488                tmp[i_] = a[offs + i, i_ + i1_];
1489              }
1490              for (j = 0; j <= i - 1; j++) {
1491                v = AP.Math.Conj(a[offs + i, offs + j]);
1492                i1_ = (0) - (offs);
1493                for (i_ = offs; i_ <= offs + j; i_++) {
1494                  a[offs + j, i_] = a[offs + j, i_] + v * tmp[i_ + i1_];
1495                }
1496              }
1497              v = AP.Math.Conj(a[offs + i, offs + i]);
1498              for (i_ = offs; i_ <= offs + i - 1; i_++) {
1499                a[offs + i, i_] = v * a[offs + i, i_];
1500              }
1501              a[offs + i, offs + i] = AP.Math.Sqr(a[offs + i, offs + i].x) + AP.Math.Sqr(a[offs + i, offs + i].y);
1502            }
1503          }
1504        }
1505        return;
1506      }
1507
1508      //
1509      // Recursive code: triangular factor inversion merged with
1510      // UU' or L'L multiplication
1511      //
1512      ablas.ablascomplexsplitlength(ref a, n, ref n1, ref n2);
1513
1514      //
1515      // form off-diagonal block of trangular inverse
1516      //
1517      if (isupper) {
1518        for (i = 0; i <= n1 - 1; i++) {
1519          for (i_ = offs + n1; i_ <= offs + n - 1; i_++) {
1520            a[offs + i, i_] = -1 * a[offs + i, i_];
1521          }
1522        }
1523        ablas.cmatrixlefttrsm(n1, n2, ref a, offs, offs, isupper, false, 0, ref a, offs, offs + n1);
1524        ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, false, 0, ref a, offs, offs + n1);
1525      } else {
1526        for (i = 0; i <= n2 - 1; i++) {
1527          for (i_ = offs; i_ <= offs + n1 - 1; i_++) {
1528            a[offs + n1 + i, i_] = -1 * a[offs + n1 + i, i_];
1529          }
1530        }
1531        ablas.cmatrixrighttrsm(n2, n1, ref a, offs, offs, isupper, false, 0, ref a, offs + n1, offs);
1532        ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, false, 0, ref a, offs + n1, offs);
1533      }
1534
1535      //
1536      // invert first diagonal block
1537      //
1538      hpdmatrixcholeskyinverserec(ref a, offs, n1, isupper, ref tmp);
1539
1540      //
1541      // update first diagonal block with off-diagonal block,
1542      // update off-diagonal block
1543      //
1544      if (isupper) {
1545        ablas.cmatrixsyrk(n1, n2, 1.0, ref a, offs, offs + n1, 0, 1.0, ref a, offs, offs, isupper);
1546        ablas.cmatrixrighttrsm(n1, n2, ref a, offs + n1, offs + n1, isupper, false, 2, ref a, offs, offs + n1);
1547      } else {
1548        ablas.cmatrixsyrk(n1, n2, 1.0, ref a, offs + n1, offs, 2, 1.0, ref a, offs, offs, isupper);
1549        ablas.cmatrixlefttrsm(n2, n1, ref a, offs + n1, offs + n1, isupper, false, 2, ref a, offs + n1, offs);
1550      }
1551
1552      //
1553      // invert second diagonal block
1554      //
1555      hpdmatrixcholeskyinverserec(ref a, offs + n1, n2, isupper, ref tmp);
1556    }
1557  }
1558}
Note: See TracBrowser for help on using the repository browser.