Free cookie consent management tool by TermsFeed Policy Generator

Ignore:
Timestamp:
07/22/10 00:44:01 (14 years ago)
Author:
swagner
Message:

Sorted usings and removed unused usings in entire solution (#1094)

Location:
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/Properties/AssemblyInfo.cs

    r3839 r4068  
    11using System.Reflection;
    2 using System.Runtime.CompilerServices;
    32using System.Runtime.InteropServices;
    43
  • trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/ablas.cs

    r3839 r4068  
    1919*************************************************************************/
    2020
    21 using System;
    22 
    23 namespace alglib
    24 {
    25     public class ablas
    26     {
    27         /*************************************************************************
    28         Splits matrix length in two parts, left part should match ABLAS block size
    29 
    30         INPUT PARAMETERS
    31             A   -   real matrix, is passed to ensure that we didn't split
    32                     complex matrix using real splitting subroutine.
    33                     matrix itself is not changed.
    34             N   -   length, N>0
    35 
    36         OUTPUT PARAMETERS
    37             N1  -   length
    38             N2  -   length
    39 
    40         N1+N2=N, N1>=N2, N2 may be zero
    41 
    42           -- ALGLIB routine --
    43              15.12.2009
    44              Bochkanov Sergey
    45         *************************************************************************/
    46         public static void ablassplitlength(ref double[,] a,
    47             int n,
    48             ref int n1,
    49             ref int n2)
    50         {
    51             if( n>ablasblocksize(ref a) )
    52             {
    53                 ablasinternalsplitlength(n, ablasblocksize(ref a), ref n1, ref n2);
    54             }
    55             else
    56             {
    57                 ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
    58             }
    59         }
    60 
    61 
    62         /*************************************************************************
    63         Complex ABLASSplitLength
    64 
    65           -- ALGLIB routine --
    66              15.12.2009
    67              Bochkanov Sergey
    68         *************************************************************************/
    69         public static void ablascomplexsplitlength(ref AP.Complex[,] a,
    70             int n,
    71             ref int n1,
    72             ref int n2)
    73         {
    74             if( n>ablascomplexblocksize(ref a) )
    75             {
    76                 ablasinternalsplitlength(n, ablascomplexblocksize(ref a), ref n1, ref n2);
    77             }
    78             else
    79             {
    80                 ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
    81             }
    82         }
    83 
    84 
    85         /*************************************************************************
    86         Returns block size - subdivision size where  cache-oblivious  soubroutines
    87         switch to the optimized kernel.
    88 
    89         INPUT PARAMETERS
    90             A   -   real matrix, is passed to ensure that we didn't split
    91                     complex matrix using real splitting subroutine.
    92                     matrix itself is not changed.
    93 
    94           -- ALGLIB routine --
    95              15.12.2009
    96              Bochkanov Sergey
    97         *************************************************************************/
    98         public static int ablasblocksize(ref double[,] a)
    99         {
    100             int result = 0;
    101 
    102             result = 32;
    103             return result;
    104         }
    105 
    106 
    107         /*************************************************************************
    108         Block size for complex subroutines.
    109 
    110           -- ALGLIB routine --
    111              15.12.2009
    112              Bochkanov Sergey
    113         *************************************************************************/
    114         public static int ablascomplexblocksize(ref AP.Complex[,] a)
    115         {
    116             int result = 0;
    117 
    118             result = 24;
    119             return result;
    120         }
    121 
    122 
    123         /*************************************************************************
    124         Microblock size
    125 
    126           -- ALGLIB routine --
    127              15.12.2009
    128              Bochkanov Sergey
    129         *************************************************************************/
    130         public static int ablasmicroblocksize()
    131         {
    132             int result = 0;
    133 
    134             result = 8;
    135             return result;
    136         }
    137 
    138 
    139         /*************************************************************************
    140         Cache-oblivous complex "copy-and-transpose"
    141 
    142         Input parameters:
    143             M   -   number of rows
    144             N   -   number of columns
    145             A   -   source matrix, MxN submatrix is copied and transposed
    146             IA  -   submatrix offset (row index)
    147             JA  -   submatrix offset (column index)
    148             A   -   destination matrix
    149             IB  -   submatrix offset (row index)
    150             JB  -   submatrix offset (column index)
    151         *************************************************************************/
    152         public static void cmatrixtranspose(int m,
    153             int n,
    154             ref AP.Complex[,] a,
    155             int ia,
    156             int ja,
    157             ref AP.Complex[,] b,
    158             int ib,
    159             int jb)
    160         {
    161             int i = 0;
    162             int s1 = 0;
    163             int s2 = 0;
    164             int i_ = 0;
    165             int i1_ = 0;
    166 
    167             if( m<=2*ablascomplexblocksize(ref a) & n<=2*ablascomplexblocksize(ref a) )
    168             {
    169                
    170                 //
    171                 // base case
    172                 //
    173                 for(i=0; i<=m-1; i++)
    174                 {
    175                     i1_ = (ja) - (ib);
    176                     for(i_=ib; i_<=ib+n-1;i_++)
    177                     {
    178                         b[i_,jb+i] = a[ia+i,i_+i1_];
    179                     }
     21
     22namespace alglib {
     23  public class ablas {
     24    /*************************************************************************
     25    Splits matrix length in two parts, left part should match ABLAS block size
     26
     27    INPUT PARAMETERS
     28        A   -   real matrix, is passed to ensure that we didn't split
     29                complex matrix using real splitting subroutine.
     30                matrix itself is not changed.
     31        N   -   length, N>0
     32
     33    OUTPUT PARAMETERS
     34        N1  -   length
     35        N2  -   length
     36
     37    N1+N2=N, N1>=N2, N2 may be zero
     38
     39      -- ALGLIB routine --
     40         15.12.2009
     41         Bochkanov Sergey
     42    *************************************************************************/
     43    public static void ablassplitlength(ref double[,] a,
     44        int n,
     45        ref int n1,
     46        ref int n2) {
     47      if (n > ablasblocksize(ref a)) {
     48        ablasinternalsplitlength(n, ablasblocksize(ref a), ref n1, ref n2);
     49      } else {
     50        ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
     51      }
     52    }
     53
     54
     55    /*************************************************************************
     56    Complex ABLASSplitLength
     57
     58      -- ALGLIB routine --
     59         15.12.2009
     60         Bochkanov Sergey
     61    *************************************************************************/
     62    public static void ablascomplexsplitlength(ref AP.Complex[,] a,
     63        int n,
     64        ref int n1,
     65        ref int n2) {
     66      if (n > ablascomplexblocksize(ref a)) {
     67        ablasinternalsplitlength(n, ablascomplexblocksize(ref a), ref n1, ref n2);
     68      } else {
     69        ablasinternalsplitlength(n, ablasmicroblocksize(), ref n1, ref n2);
     70      }
     71    }
     72
     73
     74    /*************************************************************************
     75    Returns block size - subdivision size where  cache-oblivious  soubroutines
     76    switch to the optimized kernel.
     77
     78    INPUT PARAMETERS
     79        A   -   real matrix, is passed to ensure that we didn't split
     80                complex matrix using real splitting subroutine.
     81                matrix itself is not changed.
     82
     83      -- ALGLIB routine --
     84         15.12.2009
     85         Bochkanov Sergey
     86    *************************************************************************/
     87    public static int ablasblocksize(ref double[,] a) {
     88      int result = 0;
     89
     90      result = 32;
     91      return result;
     92    }
     93
     94
     95    /*************************************************************************
     96    Block size for complex subroutines.
     97
     98      -- ALGLIB routine --
     99         15.12.2009
     100         Bochkanov Sergey
     101    *************************************************************************/
     102    public static int ablascomplexblocksize(ref AP.Complex[,] a) {
     103      int result = 0;
     104
     105      result = 24;
     106      return result;
     107    }
     108
     109
     110    /*************************************************************************
     111    Microblock size
     112
     113      -- ALGLIB routine --
     114         15.12.2009
     115         Bochkanov Sergey
     116    *************************************************************************/
     117    public static int ablasmicroblocksize() {
     118      int result = 0;
     119
     120      result = 8;
     121      return result;
     122    }
     123
     124
     125    /*************************************************************************
     126    Cache-oblivous complex "copy-and-transpose"
     127
     128    Input parameters:
     129        M   -   number of rows
     130        N   -   number of columns
     131        A   -   source matrix, MxN submatrix is copied and transposed
     132        IA  -   submatrix offset (row index)
     133        JA  -   submatrix offset (column index)
     134        A   -   destination matrix
     135        IB  -   submatrix offset (row index)
     136        JB  -   submatrix offset (column index)
     137    *************************************************************************/
     138    public static void cmatrixtranspose(int m,
     139        int n,
     140        ref AP.Complex[,] a,
     141        int ia,
     142        int ja,
     143        ref AP.Complex[,] b,
     144        int ib,
     145        int jb) {
     146      int i = 0;
     147      int s1 = 0;
     148      int s2 = 0;
     149      int i_ = 0;
     150      int i1_ = 0;
     151
     152      if (m <= 2 * ablascomplexblocksize(ref a) & n <= 2 * ablascomplexblocksize(ref a)) {
     153
     154        //
     155        // base case
     156        //
     157        for (i = 0; i <= m - 1; i++) {
     158          i1_ = (ja) - (ib);
     159          for (i_ = ib; i_ <= ib + n - 1; i_++) {
     160            b[i_, jb + i] = a[ia + i, i_ + i1_];
     161          }
     162        }
     163      } else {
     164
     165        //
     166        // Cache-oblivious recursion
     167        //
     168        if (m > n) {
     169          ablascomplexsplitlength(ref a, m, ref s1, ref s2);
     170          cmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
     171          cmatrixtranspose(s2, n, ref a, ia + s1, ja, ref b, ib, jb + s1);
     172        } else {
     173          ablascomplexsplitlength(ref a, n, ref s1, ref s2);
     174          cmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
     175          cmatrixtranspose(m, s2, ref a, ia, ja + s1, ref b, ib + s1, jb);
     176        }
     177      }
     178    }
     179
     180
     181    /*************************************************************************
     182    Cache-oblivous real "copy-and-transpose"
     183
     184    Input parameters:
     185        M   -   number of rows
     186        N   -   number of columns
     187        A   -   source matrix, MxN submatrix is copied and transposed
     188        IA  -   submatrix offset (row index)
     189        JA  -   submatrix offset (column index)
     190        A   -   destination matrix
     191        IB  -   submatrix offset (row index)
     192        JB  -   submatrix offset (column index)
     193    *************************************************************************/
     194    public static void rmatrixtranspose(int m,
     195        int n,
     196        ref double[,] a,
     197        int ia,
     198        int ja,
     199        ref double[,] b,
     200        int ib,
     201        int jb) {
     202      int i = 0;
     203      int s1 = 0;
     204      int s2 = 0;
     205      int i_ = 0;
     206      int i1_ = 0;
     207
     208      if (m <= 2 * ablasblocksize(ref a) & n <= 2 * ablasblocksize(ref a)) {
     209
     210        //
     211        // base case
     212        //
     213        for (i = 0; i <= m - 1; i++) {
     214          i1_ = (ja) - (ib);
     215          for (i_ = ib; i_ <= ib + n - 1; i_++) {
     216            b[i_, jb + i] = a[ia + i, i_ + i1_];
     217          }
     218        }
     219      } else {
     220
     221        //
     222        // Cache-oblivious recursion
     223        //
     224        if (m > n) {
     225          ablassplitlength(ref a, m, ref s1, ref s2);
     226          rmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
     227          rmatrixtranspose(s2, n, ref a, ia + s1, ja, ref b, ib, jb + s1);
     228        } else {
     229          ablassplitlength(ref a, n, ref s1, ref s2);
     230          rmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
     231          rmatrixtranspose(m, s2, ref a, ia, ja + s1, ref b, ib + s1, jb);
     232        }
     233      }
     234    }
     235
     236
     237    /*************************************************************************
     238    Copy
     239
     240    Input parameters:
     241        M   -   number of rows
     242        N   -   number of columns
     243        A   -   source matrix, MxN submatrix is copied and transposed
     244        IA  -   submatrix offset (row index)
     245        JA  -   submatrix offset (column index)
     246        B   -   destination matrix
     247        IB  -   submatrix offset (row index)
     248        JB  -   submatrix offset (column index)
     249    *************************************************************************/
     250    public static void cmatrixcopy(int m,
     251        int n,
     252        ref AP.Complex[,] a,
     253        int ia,
     254        int ja,
     255        ref AP.Complex[,] b,
     256        int ib,
     257        int jb) {
     258      int i = 0;
     259      int i_ = 0;
     260      int i1_ = 0;
     261
     262      for (i = 0; i <= m - 1; i++) {
     263        i1_ = (ja) - (jb);
     264        for (i_ = jb; i_ <= jb + n - 1; i_++) {
     265          b[ib + i, i_] = a[ia + i, i_ + i1_];
     266        }
     267      }
     268    }
     269
     270
     271    /*************************************************************************
     272    Copy
     273
     274    Input parameters:
     275        M   -   number of rows
     276        N   -   number of columns
     277        A   -   source matrix, MxN submatrix is copied and transposed
     278        IA  -   submatrix offset (row index)
     279        JA  -   submatrix offset (column index)
     280        B   -   destination matrix
     281        IB  -   submatrix offset (row index)
     282        JB  -   submatrix offset (column index)
     283    *************************************************************************/
     284    public static void rmatrixcopy(int m,
     285        int n,
     286        ref double[,] a,
     287        int ia,
     288        int ja,
     289        ref double[,] b,
     290        int ib,
     291        int jb) {
     292      int i = 0;
     293      int i_ = 0;
     294      int i1_ = 0;
     295
     296      for (i = 0; i <= m - 1; i++) {
     297        i1_ = (ja) - (jb);
     298        for (i_ = jb; i_ <= jb + n - 1; i_++) {
     299          b[ib + i, i_] = a[ia + i, i_ + i1_];
     300        }
     301      }
     302    }
     303
     304
     305    /*************************************************************************
     306    Rank-1 correction: A := A + u*v'
     307
     308    INPUT PARAMETERS:
     309        M   -   number of rows
     310        N   -   number of columns
     311        A   -   target matrix, MxN submatrix is updated
     312        IA  -   submatrix offset (row index)
     313        JA  -   submatrix offset (column index)
     314        U   -   vector #1
     315        IU  -   subvector offset
     316        V   -   vector #2
     317        IV  -   subvector offset
     318    *************************************************************************/
     319    public static void cmatrixrank1(int m,
     320        int n,
     321        ref AP.Complex[,] a,
     322        int ia,
     323        int ja,
     324        ref AP.Complex[] u,
     325        int iu,
     326        ref AP.Complex[] v,
     327        int iv) {
     328      int i = 0;
     329      AP.Complex s = 0;
     330      int i_ = 0;
     331      int i1_ = 0;
     332
     333      if (m == 0 | n == 0) {
     334        return;
     335      }
     336      if (ablasf.cmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv)) {
     337        return;
     338      }
     339      for (i = 0; i <= m - 1; i++) {
     340        s = u[iu + i];
     341        i1_ = (iv) - (ja);
     342        for (i_ = ja; i_ <= ja + n - 1; i_++) {
     343          a[ia + i, i_] = a[ia + i, i_] + s * v[i_ + i1_];
     344        }
     345      }
     346    }
     347
     348
     349    /*************************************************************************
     350    Rank-1 correction: A := A + u*v'
     351
     352    INPUT PARAMETERS:
     353        M   -   number of rows
     354        N   -   number of columns
     355        A   -   target matrix, MxN submatrix is updated
     356        IA  -   submatrix offset (row index)
     357        JA  -   submatrix offset (column index)
     358        U   -   vector #1
     359        IU  -   subvector offset
     360        V   -   vector #2
     361        IV  -   subvector offset
     362    *************************************************************************/
     363    public static void rmatrixrank1(int m,
     364        int n,
     365        ref double[,] a,
     366        int ia,
     367        int ja,
     368        ref double[] u,
     369        int iu,
     370        ref double[] v,
     371        int iv) {
     372      int i = 0;
     373      double s = 0;
     374      int i_ = 0;
     375      int i1_ = 0;
     376
     377      if (m == 0 | n == 0) {
     378        return;
     379      }
     380      if (ablasf.rmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv)) {
     381        return;
     382      }
     383      for (i = 0; i <= m - 1; i++) {
     384        s = u[iu + i];
     385        i1_ = (iv) - (ja);
     386        for (i_ = ja; i_ <= ja + n - 1; i_++) {
     387          a[ia + i, i_] = a[ia + i, i_] + s * v[i_ + i1_];
     388        }
     389      }
     390    }
     391
     392
     393    /*************************************************************************
     394    Matrix-vector product: y := op(A)*x
     395
     396    INPUT PARAMETERS:
     397        M   -   number of rows of op(A)
     398                M>=0
     399        N   -   number of columns of op(A)
     400                N>=0
     401        A   -   target matrix
     402        IA  -   submatrix offset (row index)
     403        JA  -   submatrix offset (column index)
     404        OpA -   operation type:
     405                * OpA=0     =>  op(A) = A
     406                * OpA=1     =>  op(A) = A^T
     407                * OpA=2     =>  op(A) = A^H
     408        X   -   input vector
     409        IX  -   subvector offset
     410        IY  -   subvector offset
     411
     412    OUTPUT PARAMETERS:
     413        Y   -   vector which stores result
     414
     415    if M=0, then subroutine does nothing.
     416    if N=0, Y is filled by zeros.
     417
     418
     419      -- ALGLIB routine --
     420
     421         28.01.2010
     422         Bochkanov Sergey
     423    *************************************************************************/
     424    public static void cmatrixmv(int m,
     425        int n,
     426        ref AP.Complex[,] a,
     427        int ia,
     428        int ja,
     429        int opa,
     430        ref AP.Complex[] x,
     431        int ix,
     432        ref AP.Complex[] y,
     433        int iy) {
     434      int i = 0;
     435      AP.Complex v = 0;
     436      int i_ = 0;
     437      int i1_ = 0;
     438
     439      if (m == 0) {
     440        return;
     441      }
     442      if (n == 0) {
     443        for (i = 0; i <= m - 1; i++) {
     444          y[iy + i] = 0;
     445        }
     446        return;
     447      }
     448      if (ablasf.cmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy)) {
     449        return;
     450      }
     451      if (opa == 0) {
     452
     453        //
     454        // y = A*x
     455        //
     456        for (i = 0; i <= m - 1; i++) {
     457          i1_ = (ix) - (ja);
     458          v = 0.0;
     459          for (i_ = ja; i_ <= ja + n - 1; i_++) {
     460            v += a[ia + i, i_] * x[i_ + i1_];
     461          }
     462          y[iy + i] = v;
     463        }
     464        return;
     465      }
     466      if (opa == 1) {
     467
     468        //
     469        // y = A^T*x
     470        //
     471        for (i = 0; i <= m - 1; i++) {
     472          y[iy + i] = 0;
     473        }
     474        for (i = 0; i <= n - 1; i++) {
     475          v = x[ix + i];
     476          i1_ = (ja) - (iy);
     477          for (i_ = iy; i_ <= iy + m - 1; i_++) {
     478            y[i_] = y[i_] + v * a[ia + i, i_ + i1_];
     479          }
     480        }
     481        return;
     482      }
     483      if (opa == 2) {
     484
     485        //
     486        // y = A^H*x
     487        //
     488        for (i = 0; i <= m - 1; i++) {
     489          y[iy + i] = 0;
     490        }
     491        for (i = 0; i <= n - 1; i++) {
     492          v = x[ix + i];
     493          i1_ = (ja) - (iy);
     494          for (i_ = iy; i_ <= iy + m - 1; i_++) {
     495            y[i_] = y[i_] + v * AP.Math.Conj(a[ia + i, i_ + i1_]);
     496          }
     497        }
     498        return;
     499      }
     500    }
     501
     502
     503    /*************************************************************************
     504    Matrix-vector product: y := op(A)*x
     505
     506    INPUT PARAMETERS:
     507        M   -   number of rows of op(A)
     508        N   -   number of columns of op(A)
     509        A   -   target matrix
     510        IA  -   submatrix offset (row index)
     511        JA  -   submatrix offset (column index)
     512        OpA -   operation type:
     513                * OpA=0     =>  op(A) = A
     514                * OpA=1     =>  op(A) = A^T
     515        X   -   input vector
     516        IX  -   subvector offset
     517        IY  -   subvector offset
     518
     519    OUTPUT PARAMETERS:
     520        Y   -   vector which stores result
     521
     522    if M=0, then subroutine does nothing.
     523    if N=0, Y is filled by zeros.
     524
     525
     526      -- ALGLIB routine --
     527
     528         28.01.2010
     529         Bochkanov Sergey
     530    *************************************************************************/
     531    public static void rmatrixmv(int m,
     532        int n,
     533        ref double[,] a,
     534        int ia,
     535        int ja,
     536        int opa,
     537        ref double[] x,
     538        int ix,
     539        ref double[] y,
     540        int iy) {
     541      int i = 0;
     542      double v = 0;
     543      int i_ = 0;
     544      int i1_ = 0;
     545
     546      if (m == 0) {
     547        return;
     548      }
     549      if (n == 0) {
     550        for (i = 0; i <= m - 1; i++) {
     551          y[iy + i] = 0;
     552        }
     553        return;
     554      }
     555      if (ablasf.rmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy)) {
     556        return;
     557      }
     558      if (opa == 0) {
     559
     560        //
     561        // y = A*x
     562        //
     563        for (i = 0; i <= m - 1; i++) {
     564          i1_ = (ix) - (ja);
     565          v = 0.0;
     566          for (i_ = ja; i_ <= ja + n - 1; i_++) {
     567            v += a[ia + i, i_] * x[i_ + i1_];
     568          }
     569          y[iy + i] = v;
     570        }
     571        return;
     572      }
     573      if (opa == 1) {
     574
     575        //
     576        // y = A^T*x
     577        //
     578        for (i = 0; i <= m - 1; i++) {
     579          y[iy + i] = 0;
     580        }
     581        for (i = 0; i <= n - 1; i++) {
     582          v = x[ix + i];
     583          i1_ = (ja) - (iy);
     584          for (i_ = iy; i_ <= iy + m - 1; i_++) {
     585            y[i_] = y[i_] + v * a[ia + i, i_ + i1_];
     586          }
     587        }
     588        return;
     589      }
     590    }
     591
     592
     593    /*************************************************************************
     594    This subroutine calculates X*op(A^-1) where:
     595    * X is MxN general matrix
     596    * A is NxN upper/lower triangular/unitriangular matrix
     597    * "op" may be identity transformation, transposition, conjugate transposition
     598
     599    Multiplication result replaces X.
     600    Cache-oblivious algorithm is used.
     601
     602    INPUT PARAMETERS
     603        N   -   matrix size, N>=0
     604        M   -   matrix size, N>=0
     605        A       -   matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1]
     606        I1      -   submatrix offset
     607        J1      -   submatrix offset
     608        IsUpper -   whether matrix is upper triangular
     609        IsUnit  -   whether matrix is unitriangular
     610        OpType  -   transformation type:
     611                    * 0 - no transformation
     612                    * 1 - transposition
     613                    * 2 - conjugate transposition
     614        C   -   matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
     615        I2  -   submatrix offset
     616        J2  -   submatrix offset
     617
     618      -- ALGLIB routine --
     619         15.12.2009
     620         Bochkanov Sergey
     621    *************************************************************************/
     622    public static void cmatrixrighttrsm(int m,
     623        int n,
     624        ref AP.Complex[,] a,
     625        int i1,
     626        int j1,
     627        bool isupper,
     628        bool isunit,
     629        int optype,
     630        ref AP.Complex[,] x,
     631        int i2,
     632        int j2) {
     633      int s1 = 0;
     634      int s2 = 0;
     635      int bs = 0;
     636
     637      bs = ablascomplexblocksize(ref a);
     638      if (m <= bs & n <= bs) {
     639        cmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     640        return;
     641      }
     642      if (m >= n) {
     643
     644        //
     645        // Split X: X*A = (X1 X2)^T*A
     646        //
     647        ablascomplexsplitlength(ref a, m, ref s1, ref s2);
     648        cmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     649        cmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2 + s1, j2);
     650      } else {
     651
     652        //
     653        // Split A:
     654        //               (A1  A12)
     655        // X*op(A) = X*op(       )
     656        //               (     A2)
     657        //
     658        // Different variants depending on
     659        // IsUpper/OpType combinations
     660        //
     661        ablascomplexsplitlength(ref a, n, ref s1, ref s2);
     662        if (isupper & optype == 0) {
     663
     664          //
     665          //                  (A1  A12)-1
     666          // X*A^-1 = (X1 X2)*(       )
     667          //                  (     A2)
     668          //
     669          cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     670          cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1 + s1, 0, 1.0, ref x, i2, j2 + s1);
     671          cmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     672          return;
     673        }
     674        if (isupper & optype != 0) {
     675
     676          //
     677          //                  (A1'     )-1
     678          // X*A^-1 = (X1 X2)*(        )
     679          //                  (A12' A2')
     680          //
     681          cmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     682          cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2 + s1, 0, ref a, i1, j1 + s1, optype, 1.0, ref x, i2, j2);
     683          cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     684          return;
     685        }
     686        if (!isupper & optype == 0) {
     687
     688          //
     689          //                  (A1     )-1
     690          // X*A^-1 = (X1 X2)*(       )
     691          //                  (A21  A2)
     692          //
     693          cmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     694          cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2 + s1, 0, ref a, i1 + s1, j1, 0, 1.0, ref x, i2, j2);
     695          cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     696          return;
     697        }
     698        if (!isupper & optype != 0) {
     699
     700          //
     701          //                  (A1' A21')-1
     702          // X*A^-1 = (X1 X2)*(        )
     703          //                  (     A2')
     704          //
     705          cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     706          cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1 + s1, j1, optype, 1.0, ref x, i2, j2 + s1);
     707          cmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     708          return;
     709        }
     710      }
     711    }
     712
     713
     714    /*************************************************************************
     715    This subroutine calculates op(A^-1)*X where:
     716    * X is MxN general matrix
     717    * A is MxM upper/lower triangular/unitriangular matrix
     718    * "op" may be identity transformation, transposition, conjugate transposition
     719
     720    Multiplication result replaces X.
     721    Cache-oblivious algorithm is used.
     722
     723    INPUT PARAMETERS
     724        N   -   matrix size, N>=0
     725        M   -   matrix size, N>=0
     726        A       -   matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1]
     727        I1      -   submatrix offset
     728        J1      -   submatrix offset
     729        IsUpper -   whether matrix is upper triangular
     730        IsUnit  -   whether matrix is unitriangular
     731        OpType  -   transformation type:
     732                    * 0 - no transformation
     733                    * 1 - transposition
     734                    * 2 - conjugate transposition
     735        C   -   matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
     736        I2  -   submatrix offset
     737        J2  -   submatrix offset
     738
     739      -- ALGLIB routine --
     740         15.12.2009
     741         Bochkanov Sergey
     742    *************************************************************************/
     743    public static void cmatrixlefttrsm(int m,
     744        int n,
     745        ref AP.Complex[,] a,
     746        int i1,
     747        int j1,
     748        bool isupper,
     749        bool isunit,
     750        int optype,
     751        ref AP.Complex[,] x,
     752        int i2,
     753        int j2) {
     754      int s1 = 0;
     755      int s2 = 0;
     756      int bs = 0;
     757
     758      bs = ablascomplexblocksize(ref a);
     759      if (m <= bs & n <= bs) {
     760        cmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     761        return;
     762      }
     763      if (n >= m) {
     764
     765        //
     766        // Split X: op(A)^-1*X = op(A)^-1*(X1 X2)
     767        //
     768        ablascomplexsplitlength(ref x, n, ref s1, ref s2);
     769        cmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     770        cmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2 + s1);
     771      } else {
     772
     773        //
     774        // Split A
     775        //
     776        ablascomplexsplitlength(ref a, m, ref s1, ref s2);
     777        if (isupper & optype == 0) {
     778
     779          //
     780          //           (A1  A12)-1  ( X1 )
     781          // A^-1*X* = (       )   *(    )
     782          //           (     A2)    ( X2 )
     783          //
     784          cmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     785          cmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1 + s1, 0, ref x, i2 + s1, j2, 0, 1.0, ref x, i2, j2);
     786          cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     787          return;
     788        }
     789        if (isupper & optype != 0) {
     790
     791          //
     792          //          (A1'     )-1 ( X1 )
     793          // A^-1*X = (        )  *(    )
     794          //          (A12' A2')   ( X2 )
     795          //
     796          cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     797          cmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1 + s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2 + s1, j2);
     798          cmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     799          return;
     800        }
     801        if (!isupper & optype == 0) {
     802
     803          //
     804          //          (A1     )-1 ( X1 )
     805          // A^-1*X = (       )  *(    )
     806          //          (A21  A2)   ( X2 )
     807          //
     808          cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     809          cmatrixgemm(s2, n, s1, -1.0, ref a, i1 + s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2 + s1, j2);
     810          cmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     811          return;
     812        }
     813        if (!isupper & optype != 0) {
     814
     815          //
     816          //          (A1' A21')-1 ( X1 )
     817          // A^-1*X = (        )  *(    )
     818          //          (     A2')   ( X2 )
     819          //
     820          cmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     821          cmatrixgemm(s1, n, s2, -1.0, ref a, i1 + s1, j1, optype, ref x, i2 + s1, j2, 0, 1.0, ref x, i2, j2);
     822          cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     823          return;
     824        }
     825      }
     826    }
     827
     828
     829    /*************************************************************************
     830    Same as CMatrixRightTRSM, but for real matrices
     831
     832    OpType may be only 0 or 1.
     833
     834      -- ALGLIB routine --
     835         15.12.2009
     836         Bochkanov Sergey
     837    *************************************************************************/
     838    public static void rmatrixrighttrsm(int m,
     839        int n,
     840        ref double[,] a,
     841        int i1,
     842        int j1,
     843        bool isupper,
     844        bool isunit,
     845        int optype,
     846        ref double[,] x,
     847        int i2,
     848        int j2) {
     849      int s1 = 0;
     850      int s2 = 0;
     851      int bs = 0;
     852
     853      bs = ablasblocksize(ref a);
     854      if (m <= bs & n <= bs) {
     855        rmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     856        return;
     857      }
     858      if (m >= n) {
     859
     860        //
     861        // Split X: X*A = (X1 X2)^T*A
     862        //
     863        ablassplitlength(ref a, m, ref s1, ref s2);
     864        rmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     865        rmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2 + s1, j2);
     866      } else {
     867
     868        //
     869        // Split A:
     870        //               (A1  A12)
     871        // X*op(A) = X*op(       )
     872        //               (     A2)
     873        //
     874        // Different variants depending on
     875        // IsUpper/OpType combinations
     876        //
     877        ablassplitlength(ref a, n, ref s1, ref s2);
     878        if (isupper & optype == 0) {
     879
     880          //
     881          //                  (A1  A12)-1
     882          // X*A^-1 = (X1 X2)*(       )
     883          //                  (     A2)
     884          //
     885          rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     886          rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1 + s1, 0, 1.0, ref x, i2, j2 + s1);
     887          rmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     888          return;
     889        }
     890        if (isupper & optype != 0) {
     891
     892          //
     893          //                  (A1'     )-1
     894          // X*A^-1 = (X1 X2)*(        )
     895          //                  (A12' A2')
     896          //
     897          rmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     898          rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2 + s1, 0, ref a, i1, j1 + s1, optype, 1.0, ref x, i2, j2);
     899          rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     900          return;
     901        }
     902        if (!isupper & optype == 0) {
     903
     904          //
     905          //                  (A1     )-1
     906          // X*A^-1 = (X1 X2)*(       )
     907          //                  (A21  A2)
     908          //
     909          rmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     910          rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2 + s1, 0, ref a, i1 + s1, j1, 0, 1.0, ref x, i2, j2);
     911          rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     912          return;
     913        }
     914        if (!isupper & optype != 0) {
     915
     916          //
     917          //                  (A1' A21')-1
     918          // X*A^-1 = (X1 X2)*(        )
     919          //                  (     A2')
     920          //
     921          rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     922          rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1 + s1, j1, optype, 1.0, ref x, i2, j2 + s1);
     923          rmatrixrighttrsm(m, s2, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2, j2 + s1);
     924          return;
     925        }
     926      }
     927    }
     928
     929
     930    /*************************************************************************
     931    Same as CMatrixLeftTRSM, but for real matrices
     932
     933    OpType may be only 0 or 1.
     934
     935      -- ALGLIB routine --
     936         15.12.2009
     937         Bochkanov Sergey
     938    *************************************************************************/
     939    public static void rmatrixlefttrsm(int m,
     940        int n,
     941        ref double[,] a,
     942        int i1,
     943        int j1,
     944        bool isupper,
     945        bool isunit,
     946        int optype,
     947        ref double[,] x,
     948        int i2,
     949        int j2) {
     950      int s1 = 0;
     951      int s2 = 0;
     952      int bs = 0;
     953
     954      bs = ablasblocksize(ref a);
     955      if (m <= bs & n <= bs) {
     956        rmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     957        return;
     958      }
     959      if (n >= m) {
     960
     961        //
     962        // Split X: op(A)^-1*X = op(A)^-1*(X1 X2)
     963        //
     964        ablassplitlength(ref x, n, ref s1, ref s2);
     965        rmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     966        rmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2 + s1);
     967      } else {
     968
     969        //
     970        // Split A
     971        //
     972        ablassplitlength(ref a, m, ref s1, ref s2);
     973        if (isupper & optype == 0) {
     974
     975          //
     976          //           (A1  A12)-1  ( X1 )
     977          // A^-1*X* = (       )   *(    )
     978          //           (     A2)    ( X2 )
     979          //
     980          rmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     981          rmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1 + s1, 0, ref x, i2 + s1, j2, 0, 1.0, ref x, i2, j2);
     982          rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     983          return;
     984        }
     985        if (isupper & optype != 0) {
     986
     987          //
     988          //          (A1'     )-1 ( X1 )
     989          // A^-1*X = (        )  *(    )
     990          //          (A12' A2')   ( X2 )
     991          //
     992          rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     993          rmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1 + s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2 + s1, j2);
     994          rmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     995          return;
     996        }
     997        if (!isupper & optype == 0) {
     998
     999          //
     1000          //          (A1     )-1 ( X1 )
     1001          // A^-1*X = (       )  *(    )
     1002          //          (A21  A2)   ( X2 )
     1003          //
     1004          rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     1005          rmatrixgemm(s2, n, s1, -1.0, ref a, i1 + s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2 + s1, j2);
     1006          rmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     1007          return;
     1008        }
     1009        if (!isupper & optype != 0) {
     1010
     1011          //
     1012          //          (A1' A21')-1 ( X1 )
     1013          // A^-1*X = (        )  *(    )
     1014          //          (     A2')   ( X2 )
     1015          //
     1016          rmatrixlefttrsm(s2, n, ref a, i1 + s1, j1 + s1, isupper, isunit, optype, ref x, i2 + s1, j2);
     1017          rmatrixgemm(s1, n, s2, -1.0, ref a, i1 + s1, j1, optype, ref x, i2 + s1, j2, 0, 1.0, ref x, i2, j2);
     1018          rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
     1019          return;
     1020        }
     1021      }
     1022    }
     1023
     1024
     1025    /*************************************************************************
     1026    This subroutine calculates  C=alpha*A*A^H+beta*C  or  C=alpha*A^H*A+beta*C
     1027    where:
     1028    * C is NxN Hermitian matrix given by its upper/lower triangle
     1029    * A is NxK matrix when A*A^H is calculated, KxN matrix otherwise
     1030
     1031    Additional info:
     1032    * cache-oblivious algorithm is used.
     1033    * multiplication result replaces C. If Beta=0, C elements are not used in
     1034      calculations (not multiplied by zero - just not referenced)
     1035    * if Alpha=0, A is not used (not multiplied by zero - just not referenced)
     1036    * if both Beta and Alpha are zero, C is filled by zeros.
     1037
     1038    INPUT PARAMETERS
     1039        N       -   matrix size, N>=0
     1040        K       -   matrix size, K>=0
     1041        Alpha   -   coefficient
     1042        A       -   matrix
     1043        IA      -   submatrix offset
     1044        JA      -   submatrix offset
     1045        OpTypeA -   multiplication type:
     1046                    * 0 - A*A^H is calculated
     1047                    * 2 - A^H*A is calculated
     1048        Beta    -   coefficient
     1049        C       -   matrix
     1050        IC      -   submatrix offset
     1051        JC      -   submatrix offset
     1052        IsUpper -   whether C is upper triangular or lower triangular
     1053
     1054      -- ALGLIB routine --
     1055         16.12.2009
     1056         Bochkanov Sergey
     1057    *************************************************************************/
     1058    public static void cmatrixsyrk(int n,
     1059        int k,
     1060        double alpha,
     1061        ref AP.Complex[,] a,
     1062        int ia,
     1063        int ja,
     1064        int optypea,
     1065        double beta,
     1066        ref AP.Complex[,] c,
     1067        int ic,
     1068        int jc,
     1069        bool isupper) {
     1070      int s1 = 0;
     1071      int s2 = 0;
     1072      int bs = 0;
     1073
     1074      bs = ablascomplexblocksize(ref a);
     1075      if (n <= bs & k <= bs) {
     1076        cmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1077        return;
     1078      }
     1079      if (k >= n) {
     1080
     1081        //
     1082        // Split K
     1083        //
     1084        ablascomplexsplitlength(ref a, k, ref s1, ref s2);
     1085        if (optypea == 0) {
     1086          cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1087          cmatrixsyrk(n, s2, alpha, ref a, ia, ja + s1, optypea, 1.0, ref c, ic, jc, isupper);
     1088        } else {
     1089          cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1090          cmatrixsyrk(n, s2, alpha, ref a, ia + s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
     1091        }
     1092      } else {
     1093
     1094        //
     1095        // Split N
     1096        //
     1097        ablascomplexsplitlength(ref a, n, ref s1, ref s2);
     1098        if (optypea == 0 & isupper) {
     1099          cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1100          cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia + s1, ja, 2, beta, ref c, ic, jc + s1);
     1101          cmatrixsyrk(s2, k, alpha, ref a, ia + s1, ja, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1102          return;
     1103        }
     1104        if (optypea == 0 & !isupper) {
     1105          cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1106          cmatrixgemm(s2, s1, k, alpha, ref a, ia + s1, ja, 0, ref a, ia, ja, 2, beta, ref c, ic + s1, jc);
     1107          cmatrixsyrk(s2, k, alpha, ref a, ia + s1, ja, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1108          return;
     1109        }
     1110        if (optypea != 0 & isupper) {
     1111          cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1112          cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 2, ref a, ia, ja + s1, 0, beta, ref c, ic, jc + s1);
     1113          cmatrixsyrk(s2, k, alpha, ref a, ia, ja + s1, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1114          return;
     1115        }
     1116        if (optypea != 0 & !isupper) {
     1117          cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1118          cmatrixgemm(s2, s1, k, alpha, ref a, ia, ja + s1, 2, ref a, ia, ja, 0, beta, ref c, ic + s1, jc);
     1119          cmatrixsyrk(s2, k, alpha, ref a, ia, ja + s1, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1120          return;
     1121        }
     1122      }
     1123    }
     1124
     1125
     1126    /*************************************************************************
     1127    Same as CMatrixSYRK, but for real matrices
     1128
     1129    OpType may be only 0 or 1.
     1130
     1131      -- ALGLIB routine --
     1132         16.12.2009
     1133         Bochkanov Sergey
     1134    *************************************************************************/
     1135    public static void rmatrixsyrk(int n,
     1136        int k,
     1137        double alpha,
     1138        ref double[,] a,
     1139        int ia,
     1140        int ja,
     1141        int optypea,
     1142        double beta,
     1143        ref double[,] c,
     1144        int ic,
     1145        int jc,
     1146        bool isupper) {
     1147      int s1 = 0;
     1148      int s2 = 0;
     1149      int bs = 0;
     1150
     1151      bs = ablasblocksize(ref a);
     1152      if (n <= bs & k <= bs) {
     1153        rmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1154        return;
     1155      }
     1156      if (k >= n) {
     1157
     1158        //
     1159        // Split K
     1160        //
     1161        ablassplitlength(ref a, k, ref s1, ref s2);
     1162        if (optypea == 0) {
     1163          rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1164          rmatrixsyrk(n, s2, alpha, ref a, ia, ja + s1, optypea, 1.0, ref c, ic, jc, isupper);
     1165        } else {
     1166          rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1167          rmatrixsyrk(n, s2, alpha, ref a, ia + s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
     1168        }
     1169      } else {
     1170
     1171        //
     1172        // Split N
     1173        //
     1174        ablassplitlength(ref a, n, ref s1, ref s2);
     1175        if (optypea == 0 & isupper) {
     1176          rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1177          rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia + s1, ja, 1, beta, ref c, ic, jc + s1);
     1178          rmatrixsyrk(s2, k, alpha, ref a, ia + s1, ja, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1179          return;
     1180        }
     1181        if (optypea == 0 & !isupper) {
     1182          rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1183          rmatrixgemm(s2, s1, k, alpha, ref a, ia + s1, ja, 0, ref a, ia, ja, 1, beta, ref c, ic + s1, jc);
     1184          rmatrixsyrk(s2, k, alpha, ref a, ia + s1, ja, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1185          return;
     1186        }
     1187        if (optypea != 0 & isupper) {
     1188          rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1189          rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 1, ref a, ia, ja + s1, 0, beta, ref c, ic, jc + s1);
     1190          rmatrixsyrk(s2, k, alpha, ref a, ia, ja + s1, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1191          return;
     1192        }
     1193        if (optypea != 0 & !isupper) {
     1194          rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
     1195          rmatrixgemm(s2, s1, k, alpha, ref a, ia, ja + s1, 1, ref a, ia, ja, 0, beta, ref c, ic + s1, jc);
     1196          rmatrixsyrk(s2, k, alpha, ref a, ia, ja + s1, optypea, beta, ref c, ic + s1, jc + s1, isupper);
     1197          return;
     1198        }
     1199      }
     1200    }
     1201
     1202
     1203    /*************************************************************************
     1204    This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where:
     1205    * C is MxN general matrix
     1206    * op1(A) is MxK matrix
     1207    * op2(B) is KxN matrix
     1208    * "op" may be identity transformation, transposition, conjugate transposition
     1209
     1210    Additional info:
     1211    * cache-oblivious algorithm is used.
     1212    * multiplication result replaces C. If Beta=0, C elements are not used in
     1213      calculations (not multiplied by zero - just not referenced)
     1214    * if Alpha=0, A is not used (not multiplied by zero - just not referenced)
     1215    * if both Beta and Alpha are zero, C is filled by zeros.
     1216
     1217    INPUT PARAMETERS
     1218        N       -   matrix size, N>0
     1219        M       -   matrix size, N>0
     1220        K       -   matrix size, K>0
     1221        Alpha   -   coefficient
     1222        A       -   matrix
     1223        IA      -   submatrix offset
     1224        JA      -   submatrix offset
     1225        OpTypeA -   transformation type:
     1226                    * 0 - no transformation
     1227                    * 1 - transposition
     1228                    * 2 - conjugate transposition
     1229        B       -   matrix
     1230        IB      -   submatrix offset
     1231        JB      -   submatrix offset
     1232        OpTypeB -   transformation type:
     1233                    * 0 - no transformation
     1234                    * 1 - transposition
     1235                    * 2 - conjugate transposition
     1236        Beta    -   coefficient
     1237        C       -   matrix
     1238        IC      -   submatrix offset
     1239        JC      -   submatrix offset
     1240
     1241      -- ALGLIB routine --
     1242         16.12.2009
     1243         Bochkanov Sergey
     1244    *************************************************************************/
     1245    public static void cmatrixgemm(int m,
     1246        int n,
     1247        int k,
     1248        AP.Complex alpha,
     1249        ref AP.Complex[,] a,
     1250        int ia,
     1251        int ja,
     1252        int optypea,
     1253        ref AP.Complex[,] b,
     1254        int ib,
     1255        int jb,
     1256        int optypeb,
     1257        AP.Complex beta,
     1258        ref AP.Complex[,] c,
     1259        int ic,
     1260        int jc) {
     1261      int s1 = 0;
     1262      int s2 = 0;
     1263      int bs = 0;
     1264
     1265      bs = ablascomplexblocksize(ref a);
     1266      if (m <= bs & n <= bs & k <= bs) {
     1267        cmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1268        return;
     1269      }
     1270      if (m >= n & m >= k) {
     1271
     1272        //
     1273        // A*B = (A1 A2)^T*B
     1274        //
     1275        ablascomplexsplitlength(ref a, m, ref s1, ref s2);
     1276        cmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1277        if (optypea == 0) {
     1278          cmatrixgemm(s2, n, k, alpha, ref a, ia + s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic + s1, jc);
     1279        } else {
     1280          cmatrixgemm(s2, n, k, alpha, ref a, ia, ja + s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic + s1, jc);
     1281        }
     1282        return;
     1283      }
     1284      if (n >= m & n >= k) {
     1285
     1286        //
     1287        // A*B = A*(B1 B2)
     1288        //
     1289        ablascomplexsplitlength(ref a, n, ref s1, ref s2);
     1290        if (optypeb == 0) {
     1291          cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1292          cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb + s1, optypeb, beta, ref c, ic, jc + s1);
     1293        } else {
     1294          cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1295          cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib + s1, jb, optypeb, beta, ref c, ic, jc + s1);
     1296        }
     1297        return;
     1298      }
     1299      if (k >= m & k >= n) {
     1300
     1301        //
     1302        // A*B = (A1 A2)*(B1 B2)^T
     1303        //
     1304        ablascomplexsplitlength(ref a, k, ref s1, ref s2);
     1305        if (optypea == 0 & optypeb == 0) {
     1306          cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1307          cmatrixgemm(m, n, s2, alpha, ref a, ia, ja + s1, optypea, ref b, ib + s1, jb, optypeb, 1.0, ref c, ic, jc);
     1308        }
     1309        if (optypea == 0 & optypeb != 0) {
     1310          cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1311          cmatrixgemm(m, n, s2, alpha, ref a, ia, ja + s1, optypea, ref b, ib, jb + s1, optypeb, 1.0, ref c, ic, jc);
     1312        }
     1313        if (optypea != 0 & optypeb == 0) {
     1314          cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1315          cmatrixgemm(m, n, s2, alpha, ref a, ia + s1, ja, optypea, ref b, ib + s1, jb, optypeb, 1.0, ref c, ic, jc);
     1316        }
     1317        if (optypea != 0 & optypeb != 0) {
     1318          cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1319          cmatrixgemm(m, n, s2, alpha, ref a, ia + s1, ja, optypea, ref b, ib, jb + s1, optypeb, 1.0, ref c, ic, jc);
     1320        }
     1321        return;
     1322      }
     1323    }
     1324
     1325
     1326    /*************************************************************************
     1327    Same as CMatrixGEMM, but for real numbers.
     1328    OpType may be only 0 or 1.
     1329
     1330      -- ALGLIB routine --
     1331         16.12.2009
     1332         Bochkanov Sergey
     1333    *************************************************************************/
     1334    public static void rmatrixgemm(int m,
     1335        int n,
     1336        int k,
     1337        double alpha,
     1338        ref double[,] a,
     1339        int ia,
     1340        int ja,
     1341        int optypea,
     1342        ref double[,] b,
     1343        int ib,
     1344        int jb,
     1345        int optypeb,
     1346        double beta,
     1347        ref double[,] c,
     1348        int ic,
     1349        int jc) {
     1350      int s1 = 0;
     1351      int s2 = 0;
     1352      int bs = 0;
     1353
     1354      bs = ablasblocksize(ref a);
     1355      if (m <= bs & n <= bs & k <= bs) {
     1356        rmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1357        return;
     1358      }
     1359      if (m >= n & m >= k) {
     1360
     1361        //
     1362        // A*B = (A1 A2)^T*B
     1363        //
     1364        ablassplitlength(ref a, m, ref s1, ref s2);
     1365        if (optypea == 0) {
     1366          rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1367          rmatrixgemm(s2, n, k, alpha, ref a, ia + s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic + s1, jc);
     1368        } else {
     1369          rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1370          rmatrixgemm(s2, n, k, alpha, ref a, ia, ja + s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic + s1, jc);
     1371        }
     1372        return;
     1373      }
     1374      if (n >= m & n >= k) {
     1375
     1376        //
     1377        // A*B = A*(B1 B2)
     1378        //
     1379        ablassplitlength(ref a, n, ref s1, ref s2);
     1380        if (optypeb == 0) {
     1381          rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1382          rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb + s1, optypeb, beta, ref c, ic, jc + s1);
     1383        } else {
     1384          rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1385          rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib + s1, jb, optypeb, beta, ref c, ic, jc + s1);
     1386        }
     1387        return;
     1388      }
     1389      if (k >= m & k >= n) {
     1390
     1391        //
     1392        // A*B = (A1 A2)*(B1 B2)^T
     1393        //
     1394        ablassplitlength(ref a, k, ref s1, ref s2);
     1395        if (optypea == 0 & optypeb == 0) {
     1396          rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1397          rmatrixgemm(m, n, s2, alpha, ref a, ia, ja + s1, optypea, ref b, ib + s1, jb, optypeb, 1.0, ref c, ic, jc);
     1398        }
     1399        if (optypea == 0 & optypeb != 0) {
     1400          rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1401          rmatrixgemm(m, n, s2, alpha, ref a, ia, ja + s1, optypea, ref b, ib, jb + s1, optypeb, 1.0, ref c, ic, jc);
     1402        }
     1403        if (optypea != 0 & optypeb == 0) {
     1404          rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1405          rmatrixgemm(m, n, s2, alpha, ref a, ia + s1, ja, optypea, ref b, ib + s1, jb, optypeb, 1.0, ref c, ic, jc);
     1406        }
     1407        if (optypea != 0 & optypeb != 0) {
     1408          rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
     1409          rmatrixgemm(m, n, s2, alpha, ref a, ia + s1, ja, optypea, ref b, ib, jb + s1, optypeb, 1.0, ref c, ic, jc);
     1410        }
     1411        return;
     1412      }
     1413    }
     1414
     1415
     1416    /*************************************************************************
     1417    Complex ABLASSplitLength
     1418
     1419      -- ALGLIB routine --
     1420         15.12.2009
     1421         Bochkanov Sergey
     1422    *************************************************************************/
     1423    private static void ablasinternalsplitlength(int n,
     1424        int nb,
     1425        ref int n1,
     1426        ref int n2) {
     1427      int r = 0;
     1428
     1429      if (n <= nb) {
     1430
     1431        //
     1432        // Block size, no further splitting
     1433        //
     1434        n1 = n;
     1435        n2 = 0;
     1436      } else {
     1437
     1438        //
     1439        // Greater than block size
     1440        //
     1441        if (n % nb != 0) {
     1442
     1443          //
     1444          // Split remainder
     1445          //
     1446          n2 = n % nb;
     1447          n1 = n - n2;
     1448        } else {
     1449
     1450          //
     1451          // Split on block boundaries
     1452          //
     1453          n2 = n / 2;
     1454          n1 = n - n2;
     1455          if (n1 % nb == 0) {
     1456            return;
     1457          }
     1458          r = nb - n1 % nb;
     1459          n1 = n1 + r;
     1460          n2 = n2 - r;
     1461        }
     1462      }
     1463    }
     1464
     1465
     1466    /*************************************************************************
     1467    Level 2 variant of CMatrixRightTRSM
     1468    *************************************************************************/
     1469    private static void cmatrixrighttrsm2(int m,
     1470        int n,
     1471        ref AP.Complex[,] a,
     1472        int i1,
     1473        int j1,
     1474        bool isupper,
     1475        bool isunit,
     1476        int optype,
     1477        ref AP.Complex[,] x,
     1478        int i2,
     1479        int j2) {
     1480      int i = 0;
     1481      int j = 0;
     1482      AP.Complex vc = 0;
     1483      AP.Complex vd = 0;
     1484      int i_ = 0;
     1485      int i1_ = 0;
     1486
     1487
     1488      //
     1489      // Special case
     1490      //
     1491      if (n * m == 0) {
     1492        return;
     1493      }
     1494
     1495      //
     1496      // Try to call fast TRSM
     1497      //
     1498      if (ablasf.cmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2)) {
     1499        return;
     1500      }
     1501
     1502      //
     1503      // General case
     1504      //
     1505      if (isupper) {
     1506
     1507        //
     1508        // Upper triangular matrix
     1509        //
     1510        if (optype == 0) {
     1511
     1512          //
     1513          // X*A^(-1)
     1514          //
     1515          for (i = 0; i <= m - 1; i++) {
     1516            for (j = 0; j <= n - 1; j++) {
     1517              if (isunit) {
     1518                vd = 1;
     1519              } else {
     1520                vd = a[i1 + j, j1 + j];
     1521              }
     1522              x[i2 + i, j2 + j] = x[i2 + i, j2 + j] / vd;
     1523              if (j < n - 1) {
     1524                vc = x[i2 + i, j2 + j];
     1525                i1_ = (j1 + j + 1) - (j2 + j + 1);
     1526                for (i_ = j2 + j + 1; i_ <= j2 + n - 1; i_++) {
     1527                  x[i2 + i, i_] = x[i2 + i, i_] - vc * a[i1 + j, i_ + i1_];
    1801528                }
    181             }
    182             else
    183             {
    184                
    185                 //
    186                 // Cache-oblivious recursion
    187                 //
    188                 if( m>n )
    189                 {
    190                     ablascomplexsplitlength(ref a, m, ref s1, ref s2);
    191                     cmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
    192                     cmatrixtranspose(s2, n, ref a, ia+s1, ja, ref b, ib, jb+s1);
     1529              }
     1530            }
     1531          }
     1532          return;
     1533        }
     1534        if (optype == 1) {
     1535
     1536          //
     1537          // X*A^(-T)
     1538          //
     1539          for (i = 0; i <= m - 1; i++) {
     1540            for (j = n - 1; j >= 0; j--) {
     1541              vc = 0;
     1542              vd = 1;
     1543              if (j < n - 1) {
     1544                i1_ = (j1 + j + 1) - (j2 + j + 1);
     1545                vc = 0.0;
     1546                for (i_ = j2 + j + 1; i_ <= j2 + n - 1; i_++) {
     1547                  vc += x[i2 + i, i_] * a[i1 + j, i_ + i1_];
    1931548                }
    194                 else
    195                 {
    196                     ablascomplexsplitlength(ref a, n, ref s1, ref s2);
    197                     cmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
    198                     cmatrixtranspose(m, s2, ref a, ia, ja+s1, ref b, ib+s1, jb);
     1549              }
     1550              if (!isunit) {
     1551                vd = a[i1 + j, j1 + j];
     1552              }
     1553              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vc) / vd;
     1554            }
     1555          }
     1556          return;
     1557        }
     1558        if (optype == 2) {
     1559
     1560          //
     1561          // X*A^(-H)
     1562          //
     1563          for (i = 0; i <= m - 1; i++) {
     1564            for (j = n - 1; j >= 0; j--) {
     1565              vc = 0;
     1566              vd = 1;
     1567              if (j < n - 1) {
     1568                i1_ = (j1 + j + 1) - (j2 + j + 1);
     1569                vc = 0.0;
     1570                for (i_ = j2 + j + 1; i_ <= j2 + n - 1; i_++) {
     1571                  vc += x[i2 + i, i_] * AP.Math.Conj(a[i1 + j, i_ + i1_]);
    1991572                }
    200             }
    201         }
    202 
    203 
    204         /*************************************************************************
    205         Cache-oblivous real "copy-and-transpose"
    206 
    207         Input parameters:
    208             M   -   number of rows
    209             N   -   number of columns
    210             A   -   source matrix, MxN submatrix is copied and transposed
    211             IA  -   submatrix offset (row index)
    212             JA  -   submatrix offset (column index)
    213             A   -   destination matrix
    214             IB  -   submatrix offset (row index)
    215             JB  -   submatrix offset (column index)
    216         *************************************************************************/
    217         public static void rmatrixtranspose(int m,
    218             int n,
    219             ref double[,] a,
    220             int ia,
    221             int ja,
    222             ref double[,] b,
    223             int ib,
    224             int jb)
    225         {
    226             int i = 0;
    227             int s1 = 0;
    228             int s2 = 0;
    229             int i_ = 0;
    230             int i1_ = 0;
    231 
    232             if( m<=2*ablasblocksize(ref a) & n<=2*ablasblocksize(ref a) )
    233             {
    234                
    235                 //
    236                 // base case
    237                 //
    238                 for(i=0; i<=m-1; i++)
    239                 {
    240                     i1_ = (ja) - (ib);
    241                     for(i_=ib; i_<=ib+n-1;i_++)
    242                     {
    243                         b[i_,jb+i] = a[ia+i,i_+i1_];
    244                     }
     1573              }
     1574              if (!isunit) {
     1575                vd = AP.Math.Conj(a[i1 + j, j1 + j]);
     1576              }
     1577              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vc) / vd;
     1578            }
     1579          }
     1580          return;
     1581        }
     1582      } else {
     1583
     1584        //
     1585        // Lower triangular matrix
     1586        //
     1587        if (optype == 0) {
     1588
     1589          //
     1590          // X*A^(-1)
     1591          //
     1592          for (i = 0; i <= m - 1; i++) {
     1593            for (j = n - 1; j >= 0; j--) {
     1594              if (isunit) {
     1595                vd = 1;
     1596              } else {
     1597                vd = a[i1 + j, j1 + j];
     1598              }
     1599              x[i2 + i, j2 + j] = x[i2 + i, j2 + j] / vd;
     1600              if (j > 0) {
     1601                vc = x[i2 + i, j2 + j];
     1602                i1_ = (j1) - (j2);
     1603                for (i_ = j2; i_ <= j2 + j - 1; i_++) {
     1604                  x[i2 + i, i_] = x[i2 + i, i_] - vc * a[i1 + j, i_ + i1_];
    2451605                }
    246             }
    247             else
    248             {
    249                
    250                 //
    251                 // Cache-oblivious recursion
    252                 //
    253                 if( m>n )
    254                 {
    255                     ablassplitlength(ref a, m, ref s1, ref s2);
    256                     rmatrixtranspose(s1, n, ref a, ia, ja, ref b, ib, jb);
    257                     rmatrixtranspose(s2, n, ref a, ia+s1, ja, ref b, ib, jb+s1);
     1606              }
     1607            }
     1608          }
     1609          return;
     1610        }
     1611        if (optype == 1) {
     1612
     1613          //
     1614          // X*A^(-T)
     1615          //
     1616          for (i = 0; i <= m - 1; i++) {
     1617            for (j = 0; j <= n - 1; j++) {
     1618              vc = 0;
     1619              vd = 1;
     1620              if (j > 0) {
     1621                i1_ = (j1) - (j2);
     1622                vc = 0.0;
     1623                for (i_ = j2; i_ <= j2 + j - 1; i_++) {
     1624                  vc += x[i2 + i, i_] * a[i1 + j, i_ + i1_];
    2581625                }
    259                 else
    260                 {
    261                     ablassplitlength(ref a, n, ref s1, ref s2);
    262                     rmatrixtranspose(m, s1, ref a, ia, ja, ref b, ib, jb);
    263                     rmatrixtranspose(m, s2, ref a, ia, ja+s1, ref b, ib+s1, jb);
     1626              }
     1627              if (!isunit) {
     1628                vd = a[i1 + j, j1 + j];
     1629              }
     1630              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vc) / vd;
     1631            }
     1632          }
     1633          return;
     1634        }
     1635        if (optype == 2) {
     1636
     1637          //
     1638          // X*A^(-H)
     1639          //
     1640          for (i = 0; i <= m - 1; i++) {
     1641            for (j = 0; j <= n - 1; j++) {
     1642              vc = 0;
     1643              vd = 1;
     1644              if (j > 0) {
     1645                i1_ = (j1) - (j2);
     1646                vc = 0.0;
     1647                for (i_ = j2; i_ <= j2 + j - 1; i_++) {
     1648                  vc += x[i2 + i, i_] * AP.Math.Conj(a[i1 + j, i_ + i1_]);
    2641649                }
    265             }
    266         }
    267 
    268 
    269         /*************************************************************************
    270         Copy
    271 
    272         Input parameters:
    273             M   -   number of rows
    274             N   -   number of columns
    275             A   -   source matrix, MxN submatrix is copied and transposed
    276             IA  -   submatrix offset (row index)
    277             JA  -   submatrix offset (column index)
    278             B   -   destination matrix
    279             IB  -   submatrix offset (row index)
    280             JB  -   submatrix offset (column index)
    281         *************************************************************************/
    282         public static void cmatrixcopy(int m,
    283             int n,
    284             ref AP.Complex[,] a,
    285             int ia,
    286             int ja,
    287             ref AP.Complex[,] b,
    288             int ib,
    289             int jb)
    290         {
    291             int i = 0;
    292             int i_ = 0;
    293             int i1_ = 0;
    294 
    295             for(i=0; i<=m-1; i++)
    296             {
    297                 i1_ = (ja) - (jb);
    298                 for(i_=jb; i_<=jb+n-1;i_++)
    299                 {
    300                     b[ib+i,i_] = a[ia+i,i_+i1_];
     1650              }
     1651              if (!isunit) {
     1652                vd = AP.Math.Conj(a[i1 + j, j1 + j]);
     1653              }
     1654              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vc) / vd;
     1655            }
     1656          }
     1657          return;
     1658        }
     1659      }
     1660    }
     1661
     1662
     1663    /*************************************************************************
     1664    Level-2 subroutine
     1665    *************************************************************************/
     1666    private static void cmatrixlefttrsm2(int m,
     1667        int n,
     1668        ref AP.Complex[,] a,
     1669        int i1,
     1670        int j1,
     1671        bool isupper,
     1672        bool isunit,
     1673        int optype,
     1674        ref AP.Complex[,] x,
     1675        int i2,
     1676        int j2) {
     1677      int i = 0;
     1678      int j = 0;
     1679      AP.Complex vc = 0;
     1680      AP.Complex vd = 0;
     1681      int i_ = 0;
     1682
     1683
     1684      //
     1685      // Special case
     1686      //
     1687      if (n * m == 0) {
     1688        return;
     1689      }
     1690
     1691      //
     1692      // Try to call fast TRSM
     1693      //
     1694      if (ablasf.cmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2)) {
     1695        return;
     1696      }
     1697
     1698      //
     1699      // General case
     1700      //
     1701      if (isupper) {
     1702
     1703        //
     1704        // Upper triangular matrix
     1705        //
     1706        if (optype == 0) {
     1707
     1708          //
     1709          // A^(-1)*X
     1710          //
     1711          for (i = m - 1; i >= 0; i--) {
     1712            for (j = i + 1; j <= m - 1; j++) {
     1713              vc = a[i1 + i, j1 + j];
     1714              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1715                x[i2 + i, i_] = x[i2 + i, i_] - vc * x[i2 + j, i_];
     1716              }
     1717            }
     1718            if (!isunit) {
     1719              vd = 1 / a[i1 + i, j1 + i];
     1720              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1721                x[i2 + i, i_] = vd * x[i2 + i, i_];
     1722              }
     1723            }
     1724          }
     1725          return;
     1726        }
     1727        if (optype == 1) {
     1728
     1729          //
     1730          // A^(-T)*X
     1731          //
     1732          for (i = 0; i <= m - 1; i++) {
     1733            if (isunit) {
     1734              vd = 1;
     1735            } else {
     1736              vd = 1 / a[i1 + i, j1 + i];
     1737            }
     1738            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1739              x[i2 + i, i_] = vd * x[i2 + i, i_];
     1740            }
     1741            for (j = i + 1; j <= m - 1; j++) {
     1742              vc = a[i1 + i, j1 + j];
     1743              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1744                x[i2 + j, i_] = x[i2 + j, i_] - vc * x[i2 + i, i_];
     1745              }
     1746            }
     1747          }
     1748          return;
     1749        }
     1750        if (optype == 2) {
     1751
     1752          //
     1753          // A^(-H)*X
     1754          //
     1755          for (i = 0; i <= m - 1; i++) {
     1756            if (isunit) {
     1757              vd = 1;
     1758            } else {
     1759              vd = 1 / AP.Math.Conj(a[i1 + i, j1 + i]);
     1760            }
     1761            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1762              x[i2 + i, i_] = vd * x[i2 + i, i_];
     1763            }
     1764            for (j = i + 1; j <= m - 1; j++) {
     1765              vc = AP.Math.Conj(a[i1 + i, j1 + j]);
     1766              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1767                x[i2 + j, i_] = x[i2 + j, i_] - vc * x[i2 + i, i_];
     1768              }
     1769            }
     1770          }
     1771          return;
     1772        }
     1773      } else {
     1774
     1775        //
     1776        // Lower triangular matrix
     1777        //
     1778        if (optype == 0) {
     1779
     1780          //
     1781          // A^(-1)*X
     1782          //
     1783          for (i = 0; i <= m - 1; i++) {
     1784            for (j = 0; j <= i - 1; j++) {
     1785              vc = a[i1 + i, j1 + j];
     1786              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1787                x[i2 + i, i_] = x[i2 + i, i_] - vc * x[i2 + j, i_];
     1788              }
     1789            }
     1790            if (isunit) {
     1791              vd = 1;
     1792            } else {
     1793              vd = 1 / a[i1 + j, j1 + j];
     1794            }
     1795            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1796              x[i2 + i, i_] = vd * x[i2 + i, i_];
     1797            }
     1798          }
     1799          return;
     1800        }
     1801        if (optype == 1) {
     1802
     1803          //
     1804          // A^(-T)*X
     1805          //
     1806          for (i = m - 1; i >= 0; i--) {
     1807            if (isunit) {
     1808              vd = 1;
     1809            } else {
     1810              vd = 1 / a[i1 + i, j1 + i];
     1811            }
     1812            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1813              x[i2 + i, i_] = vd * x[i2 + i, i_];
     1814            }
     1815            for (j = i - 1; j >= 0; j--) {
     1816              vc = a[i1 + i, j1 + j];
     1817              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1818                x[i2 + j, i_] = x[i2 + j, i_] - vc * x[i2 + i, i_];
     1819              }
     1820            }
     1821          }
     1822          return;
     1823        }
     1824        if (optype == 2) {
     1825
     1826          //
     1827          // A^(-H)*X
     1828          //
     1829          for (i = m - 1; i >= 0; i--) {
     1830            if (isunit) {
     1831              vd = 1;
     1832            } else {
     1833              vd = 1 / AP.Math.Conj(a[i1 + i, j1 + i]);
     1834            }
     1835            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1836              x[i2 + i, i_] = vd * x[i2 + i, i_];
     1837            }
     1838            for (j = i - 1; j >= 0; j--) {
     1839              vc = AP.Math.Conj(a[i1 + i, j1 + j]);
     1840              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     1841                x[i2 + j, i_] = x[i2 + j, i_] - vc * x[i2 + i, i_];
     1842              }
     1843            }
     1844          }
     1845          return;
     1846        }
     1847      }
     1848    }
     1849
     1850
     1851    /*************************************************************************
     1852    Level 2 subroutine
     1853
     1854      -- ALGLIB routine --
     1855         15.12.2009
     1856         Bochkanov Sergey
     1857    *************************************************************************/
     1858    private static void rmatrixrighttrsm2(int m,
     1859        int n,
     1860        ref double[,] a,
     1861        int i1,
     1862        int j1,
     1863        bool isupper,
     1864        bool isunit,
     1865        int optype,
     1866        ref double[,] x,
     1867        int i2,
     1868        int j2) {
     1869      int i = 0;
     1870      int j = 0;
     1871      double vr = 0;
     1872      double vd = 0;
     1873      int i_ = 0;
     1874      int i1_ = 0;
     1875
     1876
     1877      //
     1878      // Special case
     1879      //
     1880      if (n * m == 0) {
     1881        return;
     1882      }
     1883
     1884      //
     1885      // Try to use "fast" code
     1886      //
     1887      if (ablasf.rmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2)) {
     1888        return;
     1889      }
     1890
     1891      //
     1892      // General case
     1893      //
     1894      if (isupper) {
     1895
     1896        //
     1897        // Upper triangular matrix
     1898        //
     1899        if (optype == 0) {
     1900
     1901          //
     1902          // X*A^(-1)
     1903          //
     1904          for (i = 0; i <= m - 1; i++) {
     1905            for (j = 0; j <= n - 1; j++) {
     1906              if (isunit) {
     1907                vd = 1;
     1908              } else {
     1909                vd = a[i1 + j, j1 + j];
     1910              }
     1911              x[i2 + i, j2 + j] = x[i2 + i, j2 + j] / vd;
     1912              if (j < n - 1) {
     1913                vr = x[i2 + i, j2 + j];
     1914                i1_ = (j1 + j + 1) - (j2 + j + 1);
     1915                for (i_ = j2 + j + 1; i_ <= j2 + n - 1; i_++) {
     1916                  x[i2 + i, i_] = x[i2 + i, i_] - vr * a[i1 + j, i_ + i1_];
    3011917                }
    302             }
    303         }
    304 
    305 
    306         /*************************************************************************
    307         Copy
    308 
    309         Input parameters:
    310             M   -   number of rows
    311             N   -   number of columns
    312             A   -   source matrix, MxN submatrix is copied and transposed
    313             IA  -   submatrix offset (row index)
    314             JA  -   submatrix offset (column index)
    315             B   -   destination matrix
    316             IB  -   submatrix offset (row index)
    317             JB  -   submatrix offset (column index)
    318         *************************************************************************/
    319         public static void rmatrixcopy(int m,
    320             int n,
    321             ref double[,] a,
    322             int ia,
    323             int ja,
    324             ref double[,] b,
    325             int ib,
    326             int jb)
    327         {
    328             int i = 0;
    329             int i_ = 0;
    330             int i1_ = 0;
    331 
    332             for(i=0; i<=m-1; i++)
    333             {
    334                 i1_ = (ja) - (jb);
    335                 for(i_=jb; i_<=jb+n-1;i_++)
    336                 {
    337                     b[ib+i,i_] = a[ia+i,i_+i1_];
     1918              }
     1919            }
     1920          }
     1921          return;
     1922        }
     1923        if (optype == 1) {
     1924
     1925          //
     1926          // X*A^(-T)
     1927          //
     1928          for (i = 0; i <= m - 1; i++) {
     1929            for (j = n - 1; j >= 0; j--) {
     1930              vr = 0;
     1931              vd = 1;
     1932              if (j < n - 1) {
     1933                i1_ = (j1 + j + 1) - (j2 + j + 1);
     1934                vr = 0.0;
     1935                for (i_ = j2 + j + 1; i_ <= j2 + n - 1; i_++) {
     1936                  vr += x[i2 + i, i_] * a[i1 + j, i_ + i1_];
    3381937                }
    339             }
    340         }
    341 
    342 
    343         /*************************************************************************
    344         Rank-1 correction: A := A + u*v'
    345 
    346         INPUT PARAMETERS:
    347             M   -   number of rows
    348             N   -   number of columns
    349             A   -   target matrix, MxN submatrix is updated
    350             IA  -   submatrix offset (row index)
    351             JA  -   submatrix offset (column index)
    352             U   -   vector #1
    353             IU  -   subvector offset
    354             V   -   vector #2
    355             IV  -   subvector offset
    356         *************************************************************************/
    357         public static void cmatrixrank1(int m,
    358             int n,
    359             ref AP.Complex[,] a,
    360             int ia,
    361             int ja,
    362             ref AP.Complex[] u,
    363             int iu,
    364             ref AP.Complex[] v,
    365             int iv)
    366         {
    367             int i = 0;
    368             AP.Complex s = 0;
    369             int i_ = 0;
    370             int i1_ = 0;
    371 
    372             if( m==0 | n==0 )
    373             {
    374                 return;
    375             }
    376             if( ablasf.cmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv) )
    377             {
    378                 return;
    379             }
    380             for(i=0; i<=m-1; i++)
    381             {
    382                 s = u[iu+i];
    383                 i1_ = (iv) - (ja);
    384                 for(i_=ja; i_<=ja+n-1;i_++)
    385                 {
    386                     a[ia+i,i_] = a[ia+i,i_] + s*v[i_+i1_];
     1938              }
     1939              if (!isunit) {
     1940                vd = a[i1 + j, j1 + j];
     1941              }
     1942              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vr) / vd;
     1943            }
     1944          }
     1945          return;
     1946        }
     1947      } else {
     1948
     1949        //
     1950        // Lower triangular matrix
     1951        //
     1952        if (optype == 0) {
     1953
     1954          //
     1955          // X*A^(-1)
     1956          //
     1957          for (i = 0; i <= m - 1; i++) {
     1958            for (j = n - 1; j >= 0; j--) {
     1959              if (isunit) {
     1960                vd = 1;
     1961              } else {
     1962                vd = a[i1 + j, j1 + j];
     1963              }
     1964              x[i2 + i, j2 + j] = x[i2 + i, j2 + j] / vd;
     1965              if (j > 0) {
     1966                vr = x[i2 + i, j2 + j];
     1967                i1_ = (j1) - (j2);
     1968                for (i_ = j2; i_ <= j2 + j - 1; i_++) {
     1969                  x[i2 + i, i_] = x[i2 + i, i_] - vr * a[i1 + j, i_ + i1_];
    3871970                }
    388             }
    389         }
    390 
    391 
    392         /*************************************************************************
    393         Rank-1 correction: A := A + u*v'
    394 
    395         INPUT PARAMETERS:
    396             M   -   number of rows
    397             N   -   number of columns
    398             A   -   target matrix, MxN submatrix is updated
    399             IA  -   submatrix offset (row index)
    400             JA  -   submatrix offset (column index)
    401             U   -   vector #1
    402             IU  -   subvector offset
    403             V   -   vector #2
    404             IV  -   subvector offset
    405         *************************************************************************/
    406         public static void rmatrixrank1(int m,
    407             int n,
    408             ref double[,] a,
    409             int ia,
    410             int ja,
    411             ref double[] u,
    412             int iu,
    413             ref double[] v,
    414             int iv)
    415         {
    416             int i = 0;
    417             double s = 0;
    418             int i_ = 0;
    419             int i1_ = 0;
    420 
    421             if( m==0 | n==0 )
    422             {
    423                 return;
    424             }
    425             if( ablasf.rmatrixrank1f(m, n, ref a, ia, ja, ref u, iu, ref v, iv) )
    426             {
    427                 return;
    428             }
    429             for(i=0; i<=m-1; i++)
    430             {
    431                 s = u[iu+i];
    432                 i1_ = (iv) - (ja);
    433                 for(i_=ja; i_<=ja+n-1;i_++)
    434                 {
    435                     a[ia+i,i_] = a[ia+i,i_] + s*v[i_+i1_];
     1971              }
     1972            }
     1973          }
     1974          return;
     1975        }
     1976        if (optype == 1) {
     1977
     1978          //
     1979          // X*A^(-T)
     1980          //
     1981          for (i = 0; i <= m - 1; i++) {
     1982            for (j = 0; j <= n - 1; j++) {
     1983              vr = 0;
     1984              vd = 1;
     1985              if (j > 0) {
     1986                i1_ = (j1) - (j2);
     1987                vr = 0.0;
     1988                for (i_ = j2; i_ <= j2 + j - 1; i_++) {
     1989                  vr += x[i2 + i, i_] * a[i1 + j, i_ + i1_];
    4361990                }
    437             }
    438         }
    439 
    440 
    441         /*************************************************************************
    442         Matrix-vector product: y := op(A)*x
    443 
    444         INPUT PARAMETERS:
    445             M   -   number of rows of op(A)
    446                     M>=0
    447             N   -   number of columns of op(A)
    448                     N>=0
    449             A   -   target matrix
    450             IA  -   submatrix offset (row index)
    451             JA  -   submatrix offset (column index)
    452             OpA -   operation type:
    453                     * OpA=0     =>  op(A) = A
    454                     * OpA=1     =>  op(A) = A^T
    455                     * OpA=2     =>  op(A) = A^H
    456             X   -   input vector
    457             IX  -   subvector offset
    458             IY  -   subvector offset
    459 
    460         OUTPUT PARAMETERS:
    461             Y   -   vector which stores result
    462 
    463         if M=0, then subroutine does nothing.
    464         if N=0, Y is filled by zeros.
    465 
    466 
    467           -- ALGLIB routine --
    468 
    469              28.01.2010
    470              Bochkanov Sergey
    471         *************************************************************************/
    472         public static void cmatrixmv(int m,
    473             int n,
    474             ref AP.Complex[,] a,
    475             int ia,
    476             int ja,
    477             int opa,
    478             ref AP.Complex[] x,
    479             int ix,
    480             ref AP.Complex[] y,
    481             int iy)
    482         {
    483             int i = 0;
    484             AP.Complex v = 0;
    485             int i_ = 0;
    486             int i1_ = 0;
    487 
    488             if( m==0 )
    489             {
    490                 return;
    491             }
    492             if( n==0 )
    493             {
    494                 for(i=0; i<=m-1; i++)
    495                 {
    496                     y[iy+i] = 0;
     1991              }
     1992              if (!isunit) {
     1993                vd = a[i1 + j, j1 + j];
     1994              }
     1995              x[i2 + i, j2 + j] = (x[i2 + i, j2 + j] - vr) / vd;
     1996            }
     1997          }
     1998          return;
     1999        }
     2000      }
     2001    }
     2002
     2003
     2004    /*************************************************************************
     2005    Level 2 subroutine
     2006    *************************************************************************/
     2007    private static void rmatrixlefttrsm2(int m,
     2008        int n,
     2009        ref double[,] a,
     2010        int i1,
     2011        int j1,
     2012        bool isupper,
     2013        bool isunit,
     2014        int optype,
     2015        ref double[,] x,
     2016        int i2,
     2017        int j2) {
     2018      int i = 0;
     2019      int j = 0;
     2020      double vr = 0;
     2021      double vd = 0;
     2022      int i_ = 0;
     2023
     2024
     2025      //
     2026      // Special case
     2027      //
     2028      if (n * m == 0) {
     2029        return;
     2030      }
     2031
     2032      //
     2033      // Try fast code
     2034      //
     2035      if (ablasf.rmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2)) {
     2036        return;
     2037      }
     2038
     2039      //
     2040      // General case
     2041      //
     2042      if (isupper) {
     2043
     2044        //
     2045        // Upper triangular matrix
     2046        //
     2047        if (optype == 0) {
     2048
     2049          //
     2050          // A^(-1)*X
     2051          //
     2052          for (i = m - 1; i >= 0; i--) {
     2053            for (j = i + 1; j <= m - 1; j++) {
     2054              vr = a[i1 + i, j1 + j];
     2055              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2056                x[i2 + i, i_] = x[i2 + i, i_] - vr * x[i2 + j, i_];
     2057              }
     2058            }
     2059            if (!isunit) {
     2060              vd = 1 / a[i1 + i, j1 + i];
     2061              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2062                x[i2 + i, i_] = vd * x[i2 + i, i_];
     2063              }
     2064            }
     2065          }
     2066          return;
     2067        }
     2068        if (optype == 1) {
     2069
     2070          //
     2071          // A^(-T)*X
     2072          //
     2073          for (i = 0; i <= m - 1; i++) {
     2074            if (isunit) {
     2075              vd = 1;
     2076            } else {
     2077              vd = 1 / a[i1 + i, j1 + i];
     2078            }
     2079            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2080              x[i2 + i, i_] = vd * x[i2 + i, i_];
     2081            }
     2082            for (j = i + 1; j <= m - 1; j++) {
     2083              vr = a[i1 + i, j1 + j];
     2084              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2085                x[i2 + j, i_] = x[i2 + j, i_] - vr * x[i2 + i, i_];
     2086              }
     2087            }
     2088          }
     2089          return;
     2090        }
     2091      } else {
     2092
     2093        //
     2094        // Lower triangular matrix
     2095        //
     2096        if (optype == 0) {
     2097
     2098          //
     2099          // A^(-1)*X
     2100          //
     2101          for (i = 0; i <= m - 1; i++) {
     2102            for (j = 0; j <= i - 1; j++) {
     2103              vr = a[i1 + i, j1 + j];
     2104              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2105                x[i2 + i, i_] = x[i2 + i, i_] - vr * x[i2 + j, i_];
     2106              }
     2107            }
     2108            if (isunit) {
     2109              vd = 1;
     2110            } else {
     2111              vd = 1 / a[i1 + j, j1 + j];
     2112            }
     2113            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2114              x[i2 + i, i_] = vd * x[i2 + i, i_];
     2115            }
     2116          }
     2117          return;
     2118        }
     2119        if (optype == 1) {
     2120
     2121          //
     2122          // A^(-T)*X
     2123          //
     2124          for (i = m - 1; i >= 0; i--) {
     2125            if (isunit) {
     2126              vd = 1;
     2127            } else {
     2128              vd = 1 / a[i1 + i, j1 + i];
     2129            }
     2130            for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2131              x[i2 + i, i_] = vd * x[i2 + i, i_];
     2132            }
     2133            for (j = i - 1; j >= 0; j--) {
     2134              vr = a[i1 + i, j1 + j];
     2135              for (i_ = j2; i_ <= j2 + n - 1; i_++) {
     2136                x[i2 + j, i_] = x[i2 + j, i_] - vr * x[i2 + i, i_];
     2137              }
     2138            }
     2139          }
     2140          return;
     2141        }
     2142      }
     2143    }
     2144
     2145
     2146    /*************************************************************************
     2147    Level 2 subroutine
     2148    *************************************************************************/
     2149    private static void cmatrixsyrk2(int n,
     2150        int k,
     2151        double alpha,
     2152        ref AP.Complex[,] a,
     2153        int ia,
     2154        int ja,
     2155        int optypea,
     2156        double beta,
     2157        ref AP.Complex[,] c,
     2158        int ic,
     2159        int jc,
     2160        bool isupper) {
     2161      int i = 0;
     2162      int j = 0;
     2163      int j1 = 0;
     2164      int j2 = 0;
     2165      AP.Complex v = 0;
     2166      int i_ = 0;
     2167      int i1_ = 0;
     2168
     2169
     2170      //
     2171      // Fast exit (nothing to be done)
     2172      //
     2173      if (((double)(alpha) == (double)(0) | k == 0) & (double)(beta) == (double)(1)) {
     2174        return;
     2175      }
     2176
     2177      //
     2178      // Try to call fast SYRK
     2179      //
     2180      if (ablasf.cmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper)) {
     2181        return;
     2182      }
     2183
     2184      //
     2185      // SYRK
     2186      //
     2187      if (optypea == 0) {
     2188
     2189        //
     2190        // C=alpha*A*A^H+beta*C
     2191        //
     2192        for (i = 0; i <= n - 1; i++) {
     2193          if (isupper) {
     2194            j1 = i;
     2195            j2 = n - 1;
     2196          } else {
     2197            j1 = 0;
     2198            j2 = i;
     2199          }
     2200          for (j = j1; j <= j2; j++) {
     2201            if ((double)(alpha) != (double)(0) & k > 0) {
     2202              v = 0.0;
     2203              for (i_ = ja; i_ <= ja + k - 1; i_++) {
     2204                v += a[ia + i, i_] * AP.Math.Conj(a[ia + j, i_]);
     2205              }
     2206            } else {
     2207              v = 0;
     2208            }
     2209            if ((double)(beta) == (double)(0)) {
     2210              c[ic + i, jc + j] = alpha * v;
     2211            } else {
     2212              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
     2213            }
     2214          }
     2215        }
     2216        return;
     2217      } else {
     2218
     2219        //
     2220        // C=alpha*A^H*A+beta*C
     2221        //
     2222        for (i = 0; i <= n - 1; i++) {
     2223          if (isupper) {
     2224            j1 = i;
     2225            j2 = n - 1;
     2226          } else {
     2227            j1 = 0;
     2228            j2 = i;
     2229          }
     2230          if ((double)(beta) == (double)(0)) {
     2231            for (j = j1; j <= j2; j++) {
     2232              c[ic + i, jc + j] = 0;
     2233            }
     2234          } else {
     2235            for (i_ = jc + j1; i_ <= jc + j2; i_++) {
     2236              c[ic + i, i_] = beta * c[ic + i, i_];
     2237            }
     2238          }
     2239        }
     2240        for (i = 0; i <= k - 1; i++) {
     2241          for (j = 0; j <= n - 1; j++) {
     2242            if (isupper) {
     2243              j1 = j;
     2244              j2 = n - 1;
     2245            } else {
     2246              j1 = 0;
     2247              j2 = j;
     2248            }
     2249            v = alpha * AP.Math.Conj(a[ia + i, ja + j]);
     2250            i1_ = (ja + j1) - (jc + j1);
     2251            for (i_ = jc + j1; i_ <= jc + j2; i_++) {
     2252              c[ic + j, i_] = c[ic + j, i_] + v * a[ia + i, i_ + i1_];
     2253            }
     2254          }
     2255        }
     2256        return;
     2257      }
     2258    }
     2259
     2260
     2261    /*************************************************************************
     2262    Level 2 subrotuine
     2263    *************************************************************************/
     2264    private static void rmatrixsyrk2(int n,
     2265        int k,
     2266        double alpha,
     2267        ref double[,] a,
     2268        int ia,
     2269        int ja,
     2270        int optypea,
     2271        double beta,
     2272        ref double[,] c,
     2273        int ic,
     2274        int jc,
     2275        bool isupper) {
     2276      int i = 0;
     2277      int j = 0;
     2278      int j1 = 0;
     2279      int j2 = 0;
     2280      double v = 0;
     2281      int i_ = 0;
     2282      int i1_ = 0;
     2283
     2284
     2285      //
     2286      // Fast exit (nothing to be done)
     2287      //
     2288      if (((double)(alpha) == (double)(0) | k == 0) & (double)(beta) == (double)(1)) {
     2289        return;
     2290      }
     2291
     2292      //
     2293      // Try to call fast SYRK
     2294      //
     2295      if (ablasf.rmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper)) {
     2296        return;
     2297      }
     2298
     2299      //
     2300      // SYRK
     2301      //
     2302      if (optypea == 0) {
     2303
     2304        //
     2305        // C=alpha*A*A^H+beta*C
     2306        //
     2307        for (i = 0; i <= n - 1; i++) {
     2308          if (isupper) {
     2309            j1 = i;
     2310            j2 = n - 1;
     2311          } else {
     2312            j1 = 0;
     2313            j2 = i;
     2314          }
     2315          for (j = j1; j <= j2; j++) {
     2316            if ((double)(alpha) != (double)(0) & k > 0) {
     2317              v = 0.0;
     2318              for (i_ = ja; i_ <= ja + k - 1; i_++) {
     2319                v += a[ia + i, i_] * a[ia + j, i_];
     2320              }
     2321            } else {
     2322              v = 0;
     2323            }
     2324            if ((double)(beta) == (double)(0)) {
     2325              c[ic + i, jc + j] = alpha * v;
     2326            } else {
     2327              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
     2328            }
     2329          }
     2330        }
     2331        return;
     2332      } else {
     2333
     2334        //
     2335        // C=alpha*A^H*A+beta*C
     2336        //
     2337        for (i = 0; i <= n - 1; i++) {
     2338          if (isupper) {
     2339            j1 = i;
     2340            j2 = n - 1;
     2341          } else {
     2342            j1 = 0;
     2343            j2 = i;
     2344          }
     2345          if ((double)(beta) == (double)(0)) {
     2346            for (j = j1; j <= j2; j++) {
     2347              c[ic + i, jc + j] = 0;
     2348            }
     2349          } else {
     2350            for (i_ = jc + j1; i_ <= jc + j2; i_++) {
     2351              c[ic + i, i_] = beta * c[ic + i, i_];
     2352            }
     2353          }
     2354        }
     2355        for (i = 0; i <= k - 1; i++) {
     2356          for (j = 0; j <= n - 1; j++) {
     2357            if (isupper) {
     2358              j1 = j;
     2359              j2 = n - 1;
     2360            } else {
     2361              j1 = 0;
     2362              j2 = j;
     2363            }
     2364            v = alpha * a[ia + i, ja + j];
     2365            i1_ = (ja + j1) - (jc + j1);
     2366            for (i_ = jc + j1; i_ <= jc + j2; i_++) {
     2367              c[ic + j, i_] = c[ic + j, i_] + v * a[ia + i, i_ + i1_];
     2368            }
     2369          }
     2370        }
     2371        return;
     2372      }
     2373    }
     2374
     2375
     2376    /*************************************************************************
     2377    GEMM kernel
     2378
     2379      -- ALGLIB routine --
     2380         16.12.2009
     2381         Bochkanov Sergey
     2382    *************************************************************************/
     2383    private static void cmatrixgemmk(int m,
     2384        int n,
     2385        int k,
     2386        AP.Complex alpha,
     2387        ref AP.Complex[,] a,
     2388        int ia,
     2389        int ja,
     2390        int optypea,
     2391        ref AP.Complex[,] b,
     2392        int ib,
     2393        int jb,
     2394        int optypeb,
     2395        AP.Complex beta,
     2396        ref AP.Complex[,] c,
     2397        int ic,
     2398        int jc) {
     2399      int i = 0;
     2400      int j = 0;
     2401      AP.Complex v = 0;
     2402      int i_ = 0;
     2403      int i1_ = 0;
     2404
     2405
     2406      //
     2407      // Special case
     2408      //
     2409      if (m * n == 0) {
     2410        return;
     2411      }
     2412
     2413      //
     2414      // Try optimized code
     2415      //
     2416      if (ablasf.cmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc)) {
     2417        return;
     2418      }
     2419
     2420      //
     2421      // Another special case
     2422      //
     2423      if (k == 0) {
     2424        if (beta != 0) {
     2425          for (i = 0; i <= m - 1; i++) {
     2426            for (j = 0; j <= n - 1; j++) {
     2427              c[ic + i, jc + j] = beta * c[ic + i, jc + j];
     2428            }
     2429          }
     2430        } else {
     2431          for (i = 0; i <= m - 1; i++) {
     2432            for (j = 0; j <= n - 1; j++) {
     2433              c[ic + i, jc + j] = 0;
     2434            }
     2435          }
     2436        }
     2437        return;
     2438      }
     2439
     2440      //
     2441      // General case
     2442      //
     2443      if (optypea == 0 & optypeb != 0) {
     2444
     2445        //
     2446        // A*B'
     2447        //
     2448        for (i = 0; i <= m - 1; i++) {
     2449          for (j = 0; j <= n - 1; j++) {
     2450            if (k == 0 | alpha == 0) {
     2451              v = 0;
     2452            } else {
     2453              if (optypeb == 1) {
     2454                i1_ = (jb) - (ja);
     2455                v = 0.0;
     2456                for (i_ = ja; i_ <= ja + k - 1; i_++) {
     2457                  v += a[ia + i, i_] * b[ib + j, i_ + i1_];
    4972458                }
    498                 return;
    499             }
    500             if( ablasf.cmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy) )
    501             {
    502                 return;
    503             }
    504             if( opa==0 )
    505             {
    506                
    507                 //
    508                 // y = A*x
    509                 //
    510                 for(i=0; i<=m-1; i++)
    511                 {
    512                     i1_ = (ix)-(ja);
    513                     v = 0.0;
    514                     for(i_=ja; i_<=ja+n-1;i_++)
    515                     {
    516                         v += a[ia+i,i_]*x[i_+i1_];
    517                     }
    518                     y[iy+i] = v;
     2459              } else {
     2460                i1_ = (jb) - (ja);
     2461                v = 0.0;
     2462                for (i_ = ja; i_ <= ja + k - 1; i_++) {
     2463                  v += a[ia + i, i_] * AP.Math.Conj(b[ib + j, i_ + i1_]);
    5192464                }
    520                 return;
    521             }
    522             if( opa==1 )
    523             {
    524                
    525                 //
    526                 // y = A^T*x
    527                 //
    528                 for(i=0; i<=m-1; i++)
    529                 {
    530                     y[iy+i] = 0;
     2465              }
     2466            }
     2467            if (beta == 0) {
     2468              c[ic + i, jc + j] = alpha * v;
     2469            } else {
     2470              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
     2471            }
     2472          }
     2473        }
     2474        return;
     2475      }
     2476      if (optypea == 0 & optypeb == 0) {
     2477
     2478        //
     2479        // A*B
     2480        //
     2481        for (i = 0; i <= m - 1; i++) {
     2482          if (beta != 0) {
     2483            for (i_ = jc; i_ <= jc + n - 1; i_++) {
     2484              c[ic + i, i_] = beta * c[ic + i, i_];
     2485            }
     2486          } else {
     2487            for (j = 0; j <= n - 1; j++) {
     2488              c[ic + i, jc + j] = 0;
     2489            }
     2490          }
     2491          if (alpha != 0) {
     2492            for (j = 0; j <= k - 1; j++) {
     2493              v = alpha * a[ia + i, ja + j];
     2494              i1_ = (jb) - (jc);
     2495              for (i_ = jc; i_ <= jc + n - 1; i_++) {
     2496                c[ic + i, i_] = c[ic + i, i_] + v * b[ib + j, i_ + i1_];
     2497              }
     2498            }
     2499          }
     2500        }
     2501        return;
     2502      }
     2503      if (optypea != 0 & optypeb != 0) {
     2504
     2505        //
     2506        // A'*B'
     2507        //
     2508        for (i = 0; i <= m - 1; i++) {
     2509          for (j = 0; j <= n - 1; j++) {
     2510            if (alpha == 0) {
     2511              v = 0;
     2512            } else {
     2513              if (optypea == 1) {
     2514                if (optypeb == 1) {
     2515                  i1_ = (jb) - (ia);
     2516                  v = 0.0;
     2517                  for (i_ = ia; i_ <= ia + k - 1; i_++) {
     2518                    v += a[i_, ja + i] * b[ib + j, i_ + i1_];
     2519                  }
     2520                } else {
     2521                  i1_ = (jb) - (ia);
     2522                  v = 0.0;
     2523                  for (i_ = ia; i_ <= ia + k - 1; i_++) {
     2524                    v += a[i_, ja + i] * AP.Math.Conj(b[ib + j, i_ + i1_]);
     2525                  }
    5312526                }
    532                 for(i=0; i<=n-1; i++)
    533                 {
    534                     v = x[ix+i];
    535                     i1_ = (ja) - (iy);
    536                     for(i_=iy; i_<=iy+m-1;i_++)
    537                     {
    538                         y[i_] = y[i_] + v*a[ia+i,i_+i1_];
    539                     }
     2527              } else {
     2528                if (optypeb == 1) {
     2529                  i1_ = (jb) - (ia);
     2530                  v = 0.0;
     2531                  for (i_ = ia; i_ <= ia + k - 1; i_++) {
     2532                    v += AP.Math.Conj(a[i_, ja + i]) * b[ib + j, i_ + i1_];
     2533                  }
     2534                } else {
     2535                  i1_ = (jb) - (ia);
     2536                  v = 0.0;
     2537                  for (i_ = ia; i_ <= ia + k - 1; i_++) {
     2538                    v += AP.Math.Conj(a[i_, ja + i]) * AP.Math.Conj(b[ib + j, i_ + i1_]);
     2539                  }
    5402540                }
    541                 return;
    542             }
    543             if( opa==2 )
    544             {
    545                
    546                 //
    547                 // y = A^H*x
    548                 //
    549                 for(i=0; i<=m-1; i++)
    550                 {
    551                     y[iy+i] = 0;
    552                 }
    553                 for(i=0; i<=n-1; i++)
    554                 {
    555                     v = x[ix+i];
    556                     i1_ = (ja) - (iy);
    557                     for(i_=iy; i_<=iy+m-1;i_++)
    558                     {
    559                         y[i_] = y[i_] + v*AP.Math.Conj(a[ia+i,i_+i1_]);
    560                     }
    561                 }
    562                 return;
    563             }
    564         }
    565 
    566 
    567         /*************************************************************************
    568         Matrix-vector product: y := op(A)*x
    569 
    570         INPUT PARAMETERS:
    571             M   -   number of rows of op(A)
    572             N   -   number of columns of op(A)
    573             A   -   target matrix
    574             IA  -   submatrix offset (row index)
    575             JA  -   submatrix offset (column index)
    576             OpA -   operation type:
    577                     * OpA=0     =>  op(A) = A
    578                     * OpA=1     =>  op(A) = A^T
    579             X   -   input vector
    580             IX  -   subvector offset
    581             IY  -   subvector offset
    582 
    583         OUTPUT PARAMETERS:
    584             Y   -   vector which stores result
    585 
    586         if M=0, then subroutine does nothing.
    587         if N=0, Y is filled by zeros.
    588 
    589 
    590           -- ALGLIB routine --
    591 
    592              28.01.2010
    593              Bochkanov Sergey
    594         *************************************************************************/
    595         public static void rmatrixmv(int m,
    596             int n,
    597             ref double[,] a,
    598             int ia,
    599             int ja,
    600             int opa,
    601             ref double[] x,
    602             int ix,
    603             ref double[] y,
    604             int iy)
    605         {
    606             int i = 0;
    607             double v = 0;
    608             int i_ = 0;
    609             int i1_ = 0;
    610 
    611             if( m==0 )
    612             {
    613                 return;
    614             }
    615             if( n==0 )
    616             {
    617                 for(i=0; i<=m-1; i++)
    618                 {
    619                     y[iy+i] = 0;
    620                 }
    621                 return;
    622             }
    623             if( ablasf.rmatrixmvf(m, n, ref a, ia, ja, opa, ref x, ix, ref y, iy) )
    624             {
    625                 return;
    626             }
    627             if( opa==0 )
    628             {
    629                
    630                 //
    631                 // y = A*x
    632                 //
    633                 for(i=0; i<=m-1; i++)
    634                 {
    635                     i1_ = (ix)-(ja);
    636                     v = 0.0;
    637                     for(i_=ja; i_<=ja+n-1;i_++)
    638                     {
    639                         v += a[ia+i,i_]*x[i_+i1_];
    640                     }
    641                     y[iy+i] = v;
    642                 }
    643                 return;
    644             }
    645             if( opa==1 )
    646             {
    647                
    648                 //
    649                 // y = A^T*x
    650                 //
    651                 for(i=0; i<=m-1; i++)
    652                 {
    653                     y[iy+i] = 0;
    654                 }
    655                 for(i=0; i<=n-1; i++)
    656                 {
    657                     v = x[ix+i];
    658                     i1_ = (ja) - (iy);
    659                     for(i_=iy; i_<=iy+m-1;i_++)
    660                     {
    661                         y[i_] = y[i_] + v*a[ia+i,i_+i1_];
    662                     }
    663                 }
    664                 return;
    665             }
    666         }
    667 
    668 
    669         /*************************************************************************
    670         This subroutine calculates X*op(A^-1) where:
    671         * X is MxN general matrix
    672         * A is NxN upper/lower triangular/unitriangular matrix
    673         * "op" may be identity transformation, transposition, conjugate transposition
    674 
    675         Multiplication result replaces X.
    676         Cache-oblivious algorithm is used.
    677 
    678         INPUT PARAMETERS
    679             N   -   matrix size, N>=0
    680             M   -   matrix size, N>=0
    681             A       -   matrix, actial matrix is stored in A[I1:I1+N-1,J1:J1+N-1]
    682             I1      -   submatrix offset
    683             J1      -   submatrix offset
    684             IsUpper -   whether matrix is upper triangular
    685             IsUnit  -   whether matrix is unitriangular
    686             OpType  -   transformation type:
    687                         * 0 - no transformation
    688                         * 1 - transposition
    689                         * 2 - conjugate transposition
    690             C   -   matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
    691             I2  -   submatrix offset
    692             J2  -   submatrix offset
    693 
    694           -- ALGLIB routine --
    695              15.12.2009
    696              Bochkanov Sergey
    697         *************************************************************************/
    698         public static void cmatrixrighttrsm(int m,
    699             int n,
    700             ref AP.Complex[,] a,
    701             int i1,
    702             int j1,
    703             bool isupper,
    704             bool isunit,
    705             int optype,
    706             ref AP.Complex[,] x,
    707             int i2,
    708             int j2)
    709         {
    710             int s1 = 0;
    711             int s2 = 0;
    712             int bs = 0;
    713 
    714             bs = ablascomplexblocksize(ref a);
    715             if( m<=bs & n<=bs )
    716             {
    717                 cmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    718                 return;
    719             }
    720             if( m>=n )
    721             {
    722                
    723                 //
    724                 // Split X: X*A = (X1 X2)^T*A
    725                 //
    726                 ablascomplexsplitlength(ref a, m, ref s1, ref s2);
    727                 cmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    728                 cmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2+s1, j2);
    729             }
    730             else
    731             {
    732                
    733                 //
    734                 // Split A:
    735                 //               (A1  A12)
    736                 // X*op(A) = X*op(       )
    737                 //               (     A2)
    738                 //
    739                 // Different variants depending on
    740                 // IsUpper/OpType combinations
    741                 //
    742                 ablascomplexsplitlength(ref a, n, ref s1, ref s2);
    743                 if( isupper & optype==0 )
    744                 {
    745                    
    746                     //
    747                     //                  (A1  A12)-1
    748                     // X*A^-1 = (X1 X2)*(       )
    749                     //                  (     A2)
    750                     //
    751                     cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    752                     cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1+s1, 0, 1.0, ref x, i2, j2+s1);
    753                     cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    754                     return;
    755                 }
    756                 if( isupper & optype!=0 )
    757                 {
    758                    
    759                     //
    760                     //                  (A1'     )-1
    761                     // X*A^-1 = (X1 X2)*(        )
    762                     //                  (A12' A2')
    763                     //
    764                     cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    765                     cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1, j1+s1, optype, 1.0, ref x, i2, j2);
    766                     cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    767                     return;
    768                 }
    769                 if( !isupper & optype==0 )
    770                 {
    771                    
    772                     //
    773                     //                  (A1     )-1
    774                     // X*A^-1 = (X1 X2)*(       )
    775                     //                  (A21  A2)
    776                     //
    777                     cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    778                     cmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1+s1, j1, 0, 1.0, ref x, i2, j2);
    779                     cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    780                     return;
    781                 }
    782                 if( !isupper & optype!=0 )
    783                 {
    784                    
    785                     //
    786                     //                  (A1' A21')-1
    787                     // X*A^-1 = (X1 X2)*(        )
    788                     //                  (     A2')
    789                     //
    790                     cmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    791                     cmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1+s1, j1, optype, 1.0, ref x, i2, j2+s1);
    792                     cmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    793                     return;
    794                 }
    795             }
    796         }
    797 
    798 
    799         /*************************************************************************
    800         This subroutine calculates op(A^-1)*X where:
    801         * X is MxN general matrix
    802         * A is MxM upper/lower triangular/unitriangular matrix
    803         * "op" may be identity transformation, transposition, conjugate transposition
    804 
    805         Multiplication result replaces X.
    806         Cache-oblivious algorithm is used.
    807 
    808         INPUT PARAMETERS
    809             N   -   matrix size, N>=0
    810             M   -   matrix size, N>=0
    811             A       -   matrix, actial matrix is stored in A[I1:I1+M-1,J1:J1+M-1]
    812             I1      -   submatrix offset
    813             J1      -   submatrix offset
    814             IsUpper -   whether matrix is upper triangular
    815             IsUnit  -   whether matrix is unitriangular
    816             OpType  -   transformation type:
    817                         * 0 - no transformation
    818                         * 1 - transposition
    819                         * 2 - conjugate transposition
    820             C   -   matrix, actial matrix is stored in C[I2:I2+M-1,J2:J2+N-1]
    821             I2  -   submatrix offset
    822             J2  -   submatrix offset
    823 
    824           -- ALGLIB routine --
    825              15.12.2009
    826              Bochkanov Sergey
    827         *************************************************************************/
    828         public static void cmatrixlefttrsm(int m,
    829             int n,
    830             ref AP.Complex[,] a,
    831             int i1,
    832             int j1,
    833             bool isupper,
    834             bool isunit,
    835             int optype,
    836             ref AP.Complex[,] x,
    837             int i2,
    838             int j2)
    839         {
    840             int s1 = 0;
    841             int s2 = 0;
    842             int bs = 0;
    843 
    844             bs = ablascomplexblocksize(ref a);
    845             if( m<=bs & n<=bs )
    846             {
    847                 cmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    848                 return;
    849             }
    850             if( n>=m )
    851             {
    852                
    853                 //
    854                 // Split X: op(A)^-1*X = op(A)^-1*(X1 X2)
    855                 //
    856                 ablascomplexsplitlength(ref x, n, ref s1, ref s2);
    857                 cmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    858                 cmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2+s1);
    859             }
    860             else
    861             {
    862                
    863                 //
    864                 // Split A
    865                 //
    866                 ablascomplexsplitlength(ref a, m, ref s1, ref s2);
    867                 if( isupper & optype==0 )
    868                 {
    869                    
    870                     //
    871                     //           (A1  A12)-1  ( X1 )
    872                     // A^-1*X* = (       )   *(    )
    873                     //           (     A2)    ( X2 )
    874                     //
    875                     cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    876                     cmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1+s1, 0, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
    877                     cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    878                     return;
    879                 }
    880                 if( isupper & optype!=0 )
    881                 {
    882                    
    883                     //
    884                     //          (A1'     )-1 ( X1 )
    885                     // A^-1*X = (        )  *(    )
    886                     //          (A12' A2')   ( X2 )
    887                     //
    888                     cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    889                     cmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1+s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
    890                     cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    891                     return;
    892                 }
    893                 if( !isupper & optype==0 )
    894                 {
    895                    
    896                     //
    897                     //          (A1     )-1 ( X1 )
    898                     // A^-1*X = (       )  *(    )
    899                     //          (A21  A2)   ( X2 )
    900                     //
    901                     cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    902                     cmatrixgemm(s2, n, s1, -1.0, ref a, i1+s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
    903                     cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    904                     return;
    905                 }
    906                 if( !isupper & optype!=0 )
    907                 {
    908                    
    909                     //
    910                     //          (A1' A21')-1 ( X1 )
    911                     // A^-1*X = (        )  *(    )
    912                     //          (     A2')   ( X2 )
    913                     //
    914                     cmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    915                     cmatrixgemm(s1, n, s2, -1.0, ref a, i1+s1, j1, optype, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
    916                     cmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    917                     return;
    918                 }
    919             }
    920         }
    921 
    922 
    923         /*************************************************************************
    924         Same as CMatrixRightTRSM, but for real matrices
    925 
    926         OpType may be only 0 or 1.
    927 
    928           -- ALGLIB routine --
    929              15.12.2009
    930              Bochkanov Sergey
    931         *************************************************************************/
    932         public static void rmatrixrighttrsm(int m,
    933             int n,
    934             ref double[,] a,
    935             int i1,
    936             int j1,
    937             bool isupper,
    938             bool isunit,
    939             int optype,
    940             ref double[,] x,
    941             int i2,
    942             int j2)
    943         {
    944             int s1 = 0;
    945             int s2 = 0;
    946             int bs = 0;
    947 
    948             bs = ablasblocksize(ref a);
    949             if( m<=bs & n<=bs )
    950             {
    951                 rmatrixrighttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    952                 return;
    953             }
    954             if( m>=n )
    955             {
    956                
    957                 //
    958                 // Split X: X*A = (X1 X2)^T*A
    959                 //
    960                 ablassplitlength(ref a, m, ref s1, ref s2);
    961                 rmatrixrighttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    962                 rmatrixrighttrsm(s2, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2+s1, j2);
    963             }
    964             else
    965             {
    966                
    967                 //
    968                 // Split A:
    969                 //               (A1  A12)
    970                 // X*op(A) = X*op(       )
    971                 //               (     A2)
    972                 //
    973                 // Different variants depending on
    974                 // IsUpper/OpType combinations
    975                 //
    976                 ablassplitlength(ref a, n, ref s1, ref s2);
    977                 if( isupper & optype==0 )
    978                 {
    979                    
    980                     //
    981                     //                  (A1  A12)-1
    982                     // X*A^-1 = (X1 X2)*(       )
    983                     //                  (     A2)
    984                     //
    985                     rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    986                     rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1, j1+s1, 0, 1.0, ref x, i2, j2+s1);
    987                     rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    988                     return;
    989                 }
    990                 if( isupper & optype!=0 )
    991                 {
    992                    
    993                     //
    994                     //                  (A1'     )-1
    995                     // X*A^-1 = (X1 X2)*(        )
    996                     //                  (A12' A2')
    997                     //
    998                     rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    999                     rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1, j1+s1, optype, 1.0, ref x, i2, j2);
    1000                     rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1001                     return;
    1002                 }
    1003                 if( !isupper & optype==0 )
    1004                 {
    1005                    
    1006                     //
    1007                     //                  (A1     )-1
    1008                     // X*A^-1 = (X1 X2)*(       )
    1009                     //                  (A21  A2)
    1010                     //
    1011                     rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    1012                     rmatrixgemm(m, s1, s2, -1.0, ref x, i2, j2+s1, 0, ref a, i1+s1, j1, 0, 1.0, ref x, i2, j2);
    1013                     rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1014                     return;
    1015                 }
    1016                 if( !isupper & optype!=0 )
    1017                 {
    1018                    
    1019                     //
    1020                     //                  (A1' A21')-1
    1021                     // X*A^-1 = (X1 X2)*(        )
    1022                     //                  (     A2')
    1023                     //
    1024                     rmatrixrighttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1025                     rmatrixgemm(m, s2, s1, -1.0, ref x, i2, j2, 0, ref a, i1+s1, j1, optype, 1.0, ref x, i2, j2+s1);
    1026                     rmatrixrighttrsm(m, s2, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2, j2+s1);
    1027                     return;
    1028                 }
    1029             }
    1030         }
    1031 
    1032 
    1033         /*************************************************************************
    1034         Same as CMatrixLeftTRSM, but for real matrices
    1035 
    1036         OpType may be only 0 or 1.
    1037 
    1038           -- ALGLIB routine --
    1039              15.12.2009
    1040              Bochkanov Sergey
    1041         *************************************************************************/
    1042         public static void rmatrixlefttrsm(int m,
    1043             int n,
    1044             ref double[,] a,
    1045             int i1,
    1046             int j1,
    1047             bool isupper,
    1048             bool isunit,
    1049             int optype,
    1050             ref double[,] x,
    1051             int i2,
    1052             int j2)
    1053         {
    1054             int s1 = 0;
    1055             int s2 = 0;
    1056             int bs = 0;
    1057 
    1058             bs = ablasblocksize(ref a);
    1059             if( m<=bs & n<=bs )
    1060             {
    1061                 rmatrixlefttrsm2(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1062                 return;
    1063             }
    1064             if( n>=m )
    1065             {
    1066                
    1067                 //
    1068                 // Split X: op(A)^-1*X = op(A)^-1*(X1 X2)
    1069                 //
    1070                 ablassplitlength(ref x, n, ref s1, ref s2);
    1071                 rmatrixlefttrsm(m, s1, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1072                 rmatrixlefttrsm(m, s2, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2+s1);
    1073             }
    1074             else
    1075             {
    1076                
    1077                 //
    1078                 // Split A
    1079                 //
    1080                 ablassplitlength(ref a, m, ref s1, ref s2);
    1081                 if( isupper & optype==0 )
    1082                 {
    1083                    
    1084                     //
    1085                     //           (A1  A12)-1  ( X1 )
    1086                     // A^-1*X* = (       )   *(    )
    1087                     //           (     A2)    ( X2 )
    1088                     //
    1089                     rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    1090                     rmatrixgemm(s1, n, s2, -1.0, ref a, i1, j1+s1, 0, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
    1091                     rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1092                     return;
    1093                 }
    1094                 if( isupper & optype!=0 )
    1095                 {
    1096                    
    1097                     //
    1098                     //          (A1'     )-1 ( X1 )
    1099                     // A^-1*X = (        )  *(    )
    1100                     //          (A12' A2')   ( X2 )
    1101                     //
    1102                     rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1103                     rmatrixgemm(s2, n, s1, -1.0, ref a, i1, j1+s1, optype, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
    1104                     rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    1105                     return;
    1106                 }
    1107                 if( !isupper & optype==0 )
    1108                 {
    1109                    
    1110                     //
    1111                     //          (A1     )-1 ( X1 )
    1112                     // A^-1*X = (       )  *(    )
    1113                     //          (A21  A2)   ( X2 )
    1114                     //
    1115                     rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1116                     rmatrixgemm(s2, n, s1, -1.0, ref a, i1+s1, j1, 0, ref x, i2, j2, 0, 1.0, ref x, i2+s1, j2);
    1117                     rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    1118                     return;
    1119                 }
    1120                 if( !isupper & optype!=0 )
    1121                 {
    1122                    
    1123                     //
    1124                     //          (A1' A21')-1 ( X1 )
    1125                     // A^-1*X = (        )  *(    )
    1126                     //          (     A2')   ( X2 )
    1127                     //
    1128                     rmatrixlefttrsm(s2, n, ref a, i1+s1, j1+s1, isupper, isunit, optype, ref x, i2+s1, j2);
    1129                     rmatrixgemm(s1, n, s2, -1.0, ref a, i1+s1, j1, optype, ref x, i2+s1, j2, 0, 1.0, ref x, i2, j2);
    1130                     rmatrixlefttrsm(s1, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2);
    1131                     return;
    1132                 }
    1133             }
    1134         }
    1135 
    1136 
    1137         /*************************************************************************
    1138         This subroutine calculates  C=alpha*A*A^H+beta*C  or  C=alpha*A^H*A+beta*C
    1139         where:
    1140         * C is NxN Hermitian matrix given by its upper/lower triangle
    1141         * A is NxK matrix when A*A^H is calculated, KxN matrix otherwise
    1142 
    1143         Additional info:
    1144         * cache-oblivious algorithm is used.
    1145         * multiplication result replaces C. If Beta=0, C elements are not used in
    1146           calculations (not multiplied by zero - just not referenced)
    1147         * if Alpha=0, A is not used (not multiplied by zero - just not referenced)
    1148         * if both Beta and Alpha are zero, C is filled by zeros.
    1149 
    1150         INPUT PARAMETERS
    1151             N       -   matrix size, N>=0
    1152             K       -   matrix size, K>=0
    1153             Alpha   -   coefficient
    1154             A       -   matrix
    1155             IA      -   submatrix offset
    1156             JA      -   submatrix offset
    1157             OpTypeA -   multiplication type:
    1158                         * 0 - A*A^H is calculated
    1159                         * 2 - A^H*A is calculated
    1160             Beta    -   coefficient
    1161             C       -   matrix
    1162             IC      -   submatrix offset
    1163             JC      -   submatrix offset
    1164             IsUpper -   whether C is upper triangular or lower triangular
    1165 
    1166           -- ALGLIB routine --
    1167              16.12.2009
    1168              Bochkanov Sergey
    1169         *************************************************************************/
    1170         public static void cmatrixsyrk(int n,
    1171             int k,
    1172             double alpha,
    1173             ref AP.Complex[,] a,
    1174             int ia,
    1175             int ja,
    1176             int optypea,
    1177             double beta,
    1178             ref AP.Complex[,] c,
    1179             int ic,
    1180             int jc,
    1181             bool isupper)
    1182         {
    1183             int s1 = 0;
    1184             int s2 = 0;
    1185             int bs = 0;
    1186 
    1187             bs = ablascomplexblocksize(ref a);
    1188             if( n<=bs & k<=bs )
    1189             {
    1190                 cmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1191                 return;
    1192             }
    1193             if( k>=n )
    1194             {
    1195                
    1196                 //
    1197                 // Split K
    1198                 //
    1199                 ablascomplexsplitlength(ref a, k, ref s1, ref s2);
    1200                 if( optypea==0 )
    1201                 {
    1202                     cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1203                     cmatrixsyrk(n, s2, alpha, ref a, ia, ja+s1, optypea, 1.0, ref c, ic, jc, isupper);
    1204                 }
    1205                 else
    1206                 {
    1207                     cmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1208                     cmatrixsyrk(n, s2, alpha, ref a, ia+s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
    1209                 }
    1210             }
    1211             else
    1212             {
    1213                
    1214                 //
    1215                 // Split N
    1216                 //
    1217                 ablascomplexsplitlength(ref a, n, ref s1, ref s2);
    1218                 if( optypea==0 & isupper )
    1219                 {
    1220                     cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1221                     cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia+s1, ja, 2, beta, ref c, ic, jc+s1);
    1222                     cmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1223                     return;
    1224                 }
    1225                 if( optypea==0 & !isupper )
    1226                 {
    1227                     cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1228                     cmatrixgemm(s2, s1, k, alpha, ref a, ia+s1, ja, 0, ref a, ia, ja, 2, beta, ref c, ic+s1, jc);
    1229                     cmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1230                     return;
    1231                 }
    1232                 if( optypea!=0 & isupper )
    1233                 {
    1234                     cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1235                     cmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 2, ref a, ia, ja+s1, 0, beta, ref c, ic, jc+s1);
    1236                     cmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1237                     return;
    1238                 }
    1239                 if( optypea!=0 & !isupper )
    1240                 {
    1241                     cmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1242                     cmatrixgemm(s2, s1, k, alpha, ref a, ia, ja+s1, 2, ref a, ia, ja, 0, beta, ref c, ic+s1, jc);
    1243                     cmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1244                     return;
    1245                 }
    1246             }
    1247         }
    1248 
    1249 
    1250         /*************************************************************************
    1251         Same as CMatrixSYRK, but for real matrices
    1252 
    1253         OpType may be only 0 or 1.
    1254 
    1255           -- ALGLIB routine --
    1256              16.12.2009
    1257              Bochkanov Sergey
    1258         *************************************************************************/
    1259         public static void rmatrixsyrk(int n,
    1260             int k,
    1261             double alpha,
    1262             ref double[,] a,
    1263             int ia,
    1264             int ja,
    1265             int optypea,
    1266             double beta,
    1267             ref double[,] c,
    1268             int ic,
    1269             int jc,
    1270             bool isupper)
    1271         {
    1272             int s1 = 0;
    1273             int s2 = 0;
    1274             int bs = 0;
    1275 
    1276             bs = ablasblocksize(ref a);
    1277             if( n<=bs & k<=bs )
    1278             {
    1279                 rmatrixsyrk2(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1280                 return;
    1281             }
    1282             if( k>=n )
    1283             {
    1284                
    1285                 //
    1286                 // Split K
    1287                 //
    1288                 ablassplitlength(ref a, k, ref s1, ref s2);
    1289                 if( optypea==0 )
    1290                 {
    1291                     rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1292                     rmatrixsyrk(n, s2, alpha, ref a, ia, ja+s1, optypea, 1.0, ref c, ic, jc, isupper);
    1293                 }
    1294                 else
    1295                 {
    1296                     rmatrixsyrk(n, s1, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1297                     rmatrixsyrk(n, s2, alpha, ref a, ia+s1, ja, optypea, 1.0, ref c, ic, jc, isupper);
    1298                 }
    1299             }
    1300             else
    1301             {
    1302                
    1303                 //
    1304                 // Split N
    1305                 //
    1306                 ablassplitlength(ref a, n, ref s1, ref s2);
    1307                 if( optypea==0 & isupper )
    1308                 {
    1309                     rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1310                     rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 0, ref a, ia+s1, ja, 1, beta, ref c, ic, jc+s1);
    1311                     rmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1312                     return;
    1313                 }
    1314                 if( optypea==0 & !isupper )
    1315                 {
    1316                     rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1317                     rmatrixgemm(s2, s1, k, alpha, ref a, ia+s1, ja, 0, ref a, ia, ja, 1, beta, ref c, ic+s1, jc);
    1318                     rmatrixsyrk(s2, k, alpha, ref a, ia+s1, ja, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1319                     return;
    1320                 }
    1321                 if( optypea!=0 & isupper )
    1322                 {
    1323                     rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1324                     rmatrixgemm(s1, s2, k, alpha, ref a, ia, ja, 1, ref a, ia, ja+s1, 0, beta, ref c, ic, jc+s1);
    1325                     rmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1326                     return;
    1327                 }
    1328                 if( optypea!=0 & !isupper )
    1329                 {
    1330                     rmatrixsyrk(s1, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper);
    1331                     rmatrixgemm(s2, s1, k, alpha, ref a, ia, ja+s1, 1, ref a, ia, ja, 0, beta, ref c, ic+s1, jc);
    1332                     rmatrixsyrk(s2, k, alpha, ref a, ia, ja+s1, optypea, beta, ref c, ic+s1, jc+s1, isupper);
    1333                     return;
    1334                 }
    1335             }
    1336         }
    1337 
    1338 
    1339         /*************************************************************************
    1340         This subroutine calculates C = alpha*op1(A)*op2(B) +beta*C where:
    1341         * C is MxN general matrix
    1342         * op1(A) is MxK matrix
    1343         * op2(B) is KxN matrix
    1344         * "op" may be identity transformation, transposition, conjugate transposition
    1345 
    1346         Additional info:
    1347         * cache-oblivious algorithm is used.
    1348         * multiplication result replaces C. If Beta=0, C elements are not used in
    1349           calculations (not multiplied by zero - just not referenced)
    1350         * if Alpha=0, A is not used (not multiplied by zero - just not referenced)
    1351         * if both Beta and Alpha are zero, C is filled by zeros.
    1352 
    1353         INPUT PARAMETERS
    1354             N       -   matrix size, N>0
    1355             M       -   matrix size, N>0
    1356             K       -   matrix size, K>0
    1357             Alpha   -   coefficient
    1358             A       -   matrix
    1359             IA      -   submatrix offset
    1360             JA      -   submatrix offset
    1361             OpTypeA -   transformation type:
    1362                         * 0 - no transformation
    1363                         * 1 - transposition
    1364                         * 2 - conjugate transposition
    1365             B       -   matrix
    1366             IB      -   submatrix offset
    1367             JB      -   submatrix offset
    1368             OpTypeB -   transformation type:
    1369                         * 0 - no transformation
    1370                         * 1 - transposition
    1371                         * 2 - conjugate transposition
    1372             Beta    -   coefficient
    1373             C       -   matrix
    1374             IC      -   submatrix offset
    1375             JC      -   submatrix offset
    1376 
    1377           -- ALGLIB routine --
    1378              16.12.2009
    1379              Bochkanov Sergey
    1380         *************************************************************************/
    1381         public static void cmatrixgemm(int m,
    1382             int n,
    1383             int k,
    1384             AP.Complex alpha,
    1385             ref AP.Complex[,] a,
    1386             int ia,
    1387             int ja,
    1388             int optypea,
    1389             ref AP.Complex[,] b,
    1390             int ib,
    1391             int jb,
    1392             int optypeb,
    1393             AP.Complex beta,
    1394             ref AP.Complex[,] c,
    1395             int ic,
    1396             int jc)
    1397         {
    1398             int s1 = 0;
    1399             int s2 = 0;
    1400             int bs = 0;
    1401 
    1402             bs = ablascomplexblocksize(ref a);
    1403             if( m<=bs & n<=bs & k<=bs )
    1404             {
    1405                 cmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1406                 return;
    1407             }
    1408             if( m>=n & m>=k )
    1409             {
    1410                
    1411                 //
    1412                 // A*B = (A1 A2)^T*B
    1413                 //
    1414                 ablascomplexsplitlength(ref a, m, ref s1, ref s2);
    1415                 cmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1416                 if( optypea==0 )
    1417                 {
    1418                     cmatrixgemm(s2, n, k, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
    1419                 }
    1420                 else
    1421                 {
    1422                     cmatrixgemm(s2, n, k, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
    1423                 }
    1424                 return;
    1425             }
    1426             if( n>=m & n>=k )
    1427             {
    1428                
    1429                 //
    1430                 // A*B = A*(B1 B2)
    1431                 //
    1432                 ablascomplexsplitlength(ref a, n, ref s1, ref s2);
    1433                 if( optypeb==0 )
    1434                 {
    1435                     cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1436                     cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb+s1, optypeb, beta, ref c, ic, jc+s1);
    1437                 }
    1438                 else
    1439                 {
    1440                     cmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1441                     cmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib+s1, jb, optypeb, beta, ref c, ic, jc+s1);
    1442                 }
    1443                 return;
    1444             }
    1445             if( k>=m & k>=n )
    1446             {
    1447                
    1448                 //
    1449                 // A*B = (A1 A2)*(B1 B2)^T
    1450                 //
    1451                 ablascomplexsplitlength(ref a, k, ref s1, ref s2);
    1452                 if( optypea==0 & optypeb==0 )
    1453                 {
    1454                     cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1455                     cmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
    1456                 }
    1457                 if( optypea==0 & optypeb!=0 )
    1458                 {
    1459                     cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1460                     cmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
    1461                 }
    1462                 if( optypea!=0 & optypeb==0 )
    1463                 {
    1464                     cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1465                     cmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
    1466                 }
    1467                 if( optypea!=0 & optypeb!=0 )
    1468                 {
    1469                     cmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1470                     cmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
    1471                 }
    1472                 return;
    1473             }
    1474         }
    1475 
    1476 
    1477         /*************************************************************************
    1478         Same as CMatrixGEMM, but for real numbers.
    1479         OpType may be only 0 or 1.
    1480 
    1481           -- ALGLIB routine --
    1482              16.12.2009
    1483              Bochkanov Sergey
    1484         *************************************************************************/
    1485         public static void rmatrixgemm(int m,
    1486             int n,
    1487             int k,
    1488             double alpha,
    1489             ref double[,] a,
    1490             int ia,
    1491             int ja,
    1492             int optypea,
    1493             ref double[,] b,
    1494             int ib,
    1495             int jb,
    1496             int optypeb,
    1497             double beta,
    1498             ref double[,] c,
    1499             int ic,
    1500             int jc)
    1501         {
    1502             int s1 = 0;
    1503             int s2 = 0;
    1504             int bs = 0;
    1505 
    1506             bs = ablasblocksize(ref a);
    1507             if( m<=bs & n<=bs & k<=bs )
    1508             {
    1509                 rmatrixgemmk(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1510                 return;
    1511             }
    1512             if( m>=n & m>=k )
    1513             {
    1514                
    1515                 //
    1516                 // A*B = (A1 A2)^T*B
    1517                 //
    1518                 ablassplitlength(ref a, m, ref s1, ref s2);
    1519                 if( optypea==0 )
    1520                 {
    1521                     rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1522                     rmatrixgemm(s2, n, k, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
    1523                 }
    1524                 else
    1525                 {
    1526                     rmatrixgemm(s1, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1527                     rmatrixgemm(s2, n, k, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb, optypeb, beta, ref c, ic+s1, jc);
    1528                 }
    1529                 return;
    1530             }
    1531             if( n>=m & n>=k )
    1532             {
    1533                
    1534                 //
    1535                 // A*B = A*(B1 B2)
    1536                 //
    1537                 ablassplitlength(ref a, n, ref s1, ref s2);
    1538                 if( optypeb==0 )
    1539                 {
    1540                     rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1541                     rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb+s1, optypeb, beta, ref c, ic, jc+s1);
    1542                 }
    1543                 else
    1544                 {
    1545                     rmatrixgemm(m, s1, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1546                     rmatrixgemm(m, s2, k, alpha, ref a, ia, ja, optypea, ref b, ib+s1, jb, optypeb, beta, ref c, ic, jc+s1);
    1547                 }
    1548                 return;
    1549             }
    1550             if( k>=m & k>=n )
    1551             {
    1552                
    1553                 //
    1554                 // A*B = (A1 A2)*(B1 B2)^T
    1555                 //
    1556                 ablassplitlength(ref a, k, ref s1, ref s2);
    1557                 if( optypea==0 & optypeb==0 )
    1558                 {
    1559                     rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1560                     rmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
    1561                 }
    1562                 if( optypea==0 & optypeb!=0 )
    1563                 {
    1564                     rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1565                     rmatrixgemm(m, n, s2, alpha, ref a, ia, ja+s1, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
    1566                 }
    1567                 if( optypea!=0 & optypeb==0 )
    1568                 {
    1569                     rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1570                     rmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib+s1, jb, optypeb, 1.0, ref c, ic, jc);
    1571                 }
    1572                 if( optypea!=0 & optypeb!=0 )
    1573                 {
    1574                     rmatrixgemm(m, n, s1, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc);
    1575                     rmatrixgemm(m, n, s2, alpha, ref a, ia+s1, ja, optypea, ref b, ib, jb+s1, optypeb, 1.0, ref c, ic, jc);
    1576                 }
    1577                 return;
    1578             }
    1579         }
    1580 
    1581 
    1582         /*************************************************************************
    1583         Complex ABLASSplitLength
    1584 
    1585           -- ALGLIB routine --
    1586              15.12.2009
    1587              Bochkanov Sergey
    1588         *************************************************************************/
    1589         private static void ablasinternalsplitlength(int n,
    1590             int nb,
    1591             ref int n1,
    1592             ref int n2)
    1593         {
    1594             int r = 0;
    1595 
    1596             if( n<=nb )
    1597             {
    1598                
    1599                 //
    1600                 // Block size, no further splitting
    1601                 //
    1602                 n1 = n;
    1603                 n2 = 0;
    1604             }
    1605             else
    1606             {
    1607                
    1608                 //
    1609                 // Greater than block size
    1610                 //
    1611                 if( n%nb!=0 )
    1612                 {
    1613                    
    1614                     //
    1615                     // Split remainder
    1616                     //
    1617                     n2 = n%nb;
    1618                     n1 = n-n2;
    1619                 }
    1620                 else
    1621                 {
    1622                    
    1623                     //
    1624                     // Split on block boundaries
    1625                     //
    1626                     n2 = n/2;
    1627                     n1 = n-n2;
    1628                     if( n1%nb==0 )
    1629                     {
    1630                         return;
    1631                     }
    1632                     r = nb-n1%nb;
    1633                     n1 = n1+r;
    1634                     n2 = n2-r;
    1635                 }
    1636             }
    1637         }
    1638 
    1639 
    1640         /*************************************************************************
    1641         Level 2 variant of CMatrixRightTRSM
    1642         *************************************************************************/
    1643         private static void cmatrixrighttrsm2(int m,
    1644             int n,
    1645             ref AP.Complex[,] a,
    1646             int i1,
    1647             int j1,
    1648             bool isupper,
    1649             bool isunit,
    1650             int optype,
    1651             ref AP.Complex[,] x,
    1652             int i2,
    1653             int j2)
    1654         {
    1655             int i = 0;
    1656             int j = 0;
    1657             AP.Complex vc = 0;
    1658             AP.Complex vd = 0;
    1659             int i_ = 0;
    1660             int i1_ = 0;
    1661 
    1662            
    1663             //
    1664             // Special case
    1665             //
    1666             if( n*m==0 )
    1667             {
    1668                 return;
    1669             }
    1670            
    1671             //
    1672             // Try to call fast TRSM
    1673             //
    1674             if( ablasf.cmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
    1675             {
    1676                 return;
    1677             }
    1678            
    1679             //
    1680             // General case
    1681             //
    1682             if( isupper )
    1683             {
    1684                
    1685                 //
    1686                 // Upper triangular matrix
    1687                 //
    1688                 if( optype==0 )
    1689                 {
    1690                    
    1691                     //
    1692                     // X*A^(-1)
    1693                     //
    1694                     for(i=0; i<=m-1; i++)
    1695                     {
    1696                         for(j=0; j<=n-1; j++)
    1697                         {
    1698                             if( isunit )
    1699                             {
    1700                                 vd = 1;
    1701                             }
    1702                             else
    1703                             {
    1704                                 vd = a[i1+j,j1+j];
    1705                             }
    1706                             x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
    1707                             if( j<n-1 )
    1708                             {
    1709                                 vc = x[i2+i,j2+j];
    1710                                 i1_ = (j1+j+1) - (j2+j+1);
    1711                                 for(i_=j2+j+1; i_<=j2+n-1;i_++)
    1712                                 {
    1713                                     x[i2+i,i_] = x[i2+i,i_] - vc*a[i1+j,i_+i1_];
    1714                                 }
    1715                             }
    1716                         }
    1717                     }
    1718                     return;
    1719                 }
    1720                 if( optype==1 )
    1721                 {
    1722                    
    1723                     //
    1724                     // X*A^(-T)
    1725                     //
    1726                     for(i=0; i<=m-1; i++)
    1727                     {
    1728                         for(j=n-1; j>=0; j--)
    1729                         {
    1730                             vc = 0;
    1731                             vd = 1;
    1732                             if( j<n-1 )
    1733                             {
    1734                                 i1_ = (j1+j+1)-(j2+j+1);
    1735                                 vc = 0.0;
    1736                                 for(i_=j2+j+1; i_<=j2+n-1;i_++)
    1737                                 {
    1738                                     vc += x[i2+i,i_]*a[i1+j,i_+i1_];
    1739                                 }
    1740                             }
    1741                             if( !isunit )
    1742                             {
    1743                                 vd = a[i1+j,j1+j];
    1744                             }
    1745                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
    1746                         }
    1747                     }
    1748                     return;
    1749                 }
    1750                 if( optype==2 )
    1751                 {
    1752                    
    1753                     //
    1754                     // X*A^(-H)
    1755                     //
    1756                     for(i=0; i<=m-1; i++)
    1757                     {
    1758                         for(j=n-1; j>=0; j--)
    1759                         {
    1760                             vc = 0;
    1761                             vd = 1;
    1762                             if( j<n-1 )
    1763                             {
    1764                                 i1_ = (j1+j+1)-(j2+j+1);
    1765                                 vc = 0.0;
    1766                                 for(i_=j2+j+1; i_<=j2+n-1;i_++)
    1767                                 {
    1768                                     vc += x[i2+i,i_]*AP.Math.Conj(a[i1+j,i_+i1_]);
    1769                                 }
    1770                             }
    1771                             if( !isunit )
    1772                             {
    1773                                 vd = AP.Math.Conj(a[i1+j,j1+j]);
    1774                             }
    1775                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
    1776                         }
    1777                     }
    1778                     return;
    1779                 }
    1780             }
    1781             else
    1782             {
    1783                
    1784                 //
    1785                 // Lower triangular matrix
    1786                 //
    1787                 if( optype==0 )
    1788                 {
    1789                    
    1790                     //
    1791                     // X*A^(-1)
    1792                     //
    1793                     for(i=0; i<=m-1; i++)
    1794                     {
    1795                         for(j=n-1; j>=0; j--)
    1796                         {
    1797                             if( isunit )
    1798                             {
    1799                                 vd = 1;
    1800                             }
    1801                             else
    1802                             {
    1803                                 vd = a[i1+j,j1+j];
    1804                             }
    1805                             x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
    1806                             if( j>0 )
    1807                             {
    1808                                 vc = x[i2+i,j2+j];
    1809                                 i1_ = (j1) - (j2);
    1810                                 for(i_=j2; i_<=j2+j-1;i_++)
    1811                                 {
    1812                                     x[i2+i,i_] = x[i2+i,i_] - vc*a[i1+j,i_+i1_];
    1813                                 }
    1814                             }
    1815                         }
    1816                     }
    1817                     return;
    1818                 }
    1819                 if( optype==1 )
    1820                 {
    1821                    
    1822                     //
    1823                     // X*A^(-T)
    1824                     //
    1825                     for(i=0; i<=m-1; i++)
    1826                     {
    1827                         for(j=0; j<=n-1; j++)
    1828                         {
    1829                             vc = 0;
    1830                             vd = 1;
    1831                             if( j>0 )
    1832                             {
    1833                                 i1_ = (j1)-(j2);
    1834                                 vc = 0.0;
    1835                                 for(i_=j2; i_<=j2+j-1;i_++)
    1836                                 {
    1837                                     vc += x[i2+i,i_]*a[i1+j,i_+i1_];
    1838                                 }
    1839                             }
    1840                             if( !isunit )
    1841                             {
    1842                                 vd = a[i1+j,j1+j];
    1843                             }
    1844                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
    1845                         }
    1846                     }
    1847                     return;
    1848                 }
    1849                 if( optype==2 )
    1850                 {
    1851                    
    1852                     //
    1853                     // X*A^(-H)
    1854                     //
    1855                     for(i=0; i<=m-1; i++)
    1856                     {
    1857                         for(j=0; j<=n-1; j++)
    1858                         {
    1859                             vc = 0;
    1860                             vd = 1;
    1861                             if( j>0 )
    1862                             {
    1863                                 i1_ = (j1)-(j2);
    1864                                 vc = 0.0;
    1865                                 for(i_=j2; i_<=j2+j-1;i_++)
    1866                                 {
    1867                                     vc += x[i2+i,i_]*AP.Math.Conj(a[i1+j,i_+i1_]);
    1868                                 }
    1869                             }
    1870                             if( !isunit )
    1871                             {
    1872                                 vd = AP.Math.Conj(a[i1+j,j1+j]);
    1873                             }
    1874                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vc)/vd;
    1875                         }
    1876                     }
    1877                     return;
    1878                 }
    1879             }
    1880         }
    1881 
    1882 
    1883         /*************************************************************************
    1884         Level-2 subroutine
    1885         *************************************************************************/
    1886         private static void cmatrixlefttrsm2(int m,
    1887             int n,
    1888             ref AP.Complex[,] a,
    1889             int i1,
    1890             int j1,
    1891             bool isupper,
    1892             bool isunit,
    1893             int optype,
    1894             ref AP.Complex[,] x,
    1895             int i2,
    1896             int j2)
    1897         {
    1898             int i = 0;
    1899             int j = 0;
    1900             AP.Complex vc = 0;
    1901             AP.Complex vd = 0;
    1902             int i_ = 0;
    1903 
    1904            
    1905             //
    1906             // Special case
    1907             //
    1908             if( n*m==0 )
    1909             {
    1910                 return;
    1911             }
    1912            
    1913             //
    1914             // Try to call fast TRSM
    1915             //
    1916             if( ablasf.cmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
    1917             {
    1918                 return;
    1919             }
    1920            
    1921             //
    1922             // General case
    1923             //
    1924             if( isupper )
    1925             {
    1926                
    1927                 //
    1928                 // Upper triangular matrix
    1929                 //
    1930                 if( optype==0 )
    1931                 {
    1932                    
    1933                     //
    1934                     // A^(-1)*X
    1935                     //
    1936                     for(i=m-1; i>=0; i--)
    1937                     {
    1938                         for(j=i+1; j<=m-1; j++)
    1939                         {
    1940                             vc = a[i1+i,j1+j];
    1941                             for(i_=j2; i_<=j2+n-1;i_++)
    1942                             {
    1943                                 x[i2+i,i_] = x[i2+i,i_] - vc*x[i2+j,i_];
    1944                             }
    1945                         }
    1946                         if( !isunit )
    1947                         {
    1948                             vd = 1/a[i1+i,j1+i];
    1949                             for(i_=j2; i_<=j2+n-1;i_++)
    1950                             {
    1951                                 x[i2+i,i_] = vd*x[i2+i,i_];
    1952                             }
    1953                         }
    1954                     }
    1955                     return;
    1956                 }
    1957                 if( optype==1 )
    1958                 {
    1959                    
    1960                     //
    1961                     // A^(-T)*X
    1962                     //
    1963                     for(i=0; i<=m-1; i++)
    1964                     {
    1965                         if( isunit )
    1966                         {
    1967                             vd = 1;
    1968                         }
    1969                         else
    1970                         {
    1971                             vd = 1/a[i1+i,j1+i];
    1972                         }
    1973                         for(i_=j2; i_<=j2+n-1;i_++)
    1974                         {
    1975                             x[i2+i,i_] = vd*x[i2+i,i_];
    1976                         }
    1977                         for(j=i+1; j<=m-1; j++)
    1978                         {
    1979                             vc = a[i1+i,j1+j];
    1980                             for(i_=j2; i_<=j2+n-1;i_++)
    1981                             {
    1982                                 x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
    1983                             }
    1984                         }
    1985                     }
    1986                     return;
    1987                 }
    1988                 if( optype==2 )
    1989                 {
    1990                    
    1991                     //
    1992                     // A^(-H)*X
    1993                     //
    1994                     for(i=0; i<=m-1; i++)
    1995                     {
    1996                         if( isunit )
    1997                         {
    1998                             vd = 1;
    1999                         }
    2000                         else
    2001                         {
    2002                             vd = 1/AP.Math.Conj(a[i1+i,j1+i]);
    2003                         }
    2004                         for(i_=j2; i_<=j2+n-1;i_++)
    2005                         {
    2006                             x[i2+i,i_] = vd*x[i2+i,i_];
    2007                         }
    2008                         for(j=i+1; j<=m-1; j++)
    2009                         {
    2010                             vc = AP.Math.Conj(a[i1+i,j1+j]);
    2011                             for(i_=j2; i_<=j2+n-1;i_++)
    2012                             {
    2013                                 x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
    2014                             }
    2015                         }
    2016                     }
    2017                     return;
    2018                 }
    2019             }
    2020             else
    2021             {
    2022                
    2023                 //
    2024                 // Lower triangular matrix
    2025                 //
    2026                 if( optype==0 )
    2027                 {
    2028                    
    2029                     //
    2030                     // A^(-1)*X
    2031                     //
    2032                     for(i=0; i<=m-1; i++)
    2033                     {
    2034                         for(j=0; j<=i-1; j++)
    2035                         {
    2036                             vc = a[i1+i,j1+j];
    2037                             for(i_=j2; i_<=j2+n-1;i_++)
    2038                             {
    2039                                 x[i2+i,i_] = x[i2+i,i_] - vc*x[i2+j,i_];
    2040                             }
    2041                         }
    2042                         if( isunit )
    2043                         {
    2044                             vd = 1;
    2045                         }
    2046                         else
    2047                         {
    2048                             vd = 1/a[i1+j,j1+j];
    2049                         }
    2050                         for(i_=j2; i_<=j2+n-1;i_++)
    2051                         {
    2052                             x[i2+i,i_] = vd*x[i2+i,i_];
    2053                         }
    2054                     }
    2055                     return;
    2056                 }
    2057                 if( optype==1 )
    2058                 {
    2059                    
    2060                     //
    2061                     // A^(-T)*X
    2062                     //
    2063                     for(i=m-1; i>=0; i--)
    2064                     {
    2065                         if( isunit )
    2066                         {
    2067                             vd = 1;
    2068                         }
    2069                         else
    2070                         {
    2071                             vd = 1/a[i1+i,j1+i];
    2072                         }
    2073                         for(i_=j2; i_<=j2+n-1;i_++)
    2074                         {
    2075                             x[i2+i,i_] = vd*x[i2+i,i_];
    2076                         }
    2077                         for(j=i-1; j>=0; j--)
    2078                         {
    2079                             vc = a[i1+i,j1+j];
    2080                             for(i_=j2; i_<=j2+n-1;i_++)
    2081                             {
    2082                                 x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
    2083                             }
    2084                         }
    2085                     }
    2086                     return;
    2087                 }
    2088                 if( optype==2 )
    2089                 {
    2090                    
    2091                     //
    2092                     // A^(-H)*X
    2093                     //
    2094                     for(i=m-1; i>=0; i--)
    2095                     {
    2096                         if( isunit )
    2097                         {
    2098                             vd = 1;
    2099                         }
    2100                         else
    2101                         {
    2102                             vd = 1/AP.Math.Conj(a[i1+i,j1+i]);
    2103                         }
    2104                         for(i_=j2; i_<=j2+n-1;i_++)
    2105                         {
    2106                             x[i2+i,i_] = vd*x[i2+i,i_];
    2107                         }
    2108                         for(j=i-1; j>=0; j--)
    2109                         {
    2110                             vc = AP.Math.Conj(a[i1+i,j1+j]);
    2111                             for(i_=j2; i_<=j2+n-1;i_++)
    2112                             {
    2113                                 x[i2+j,i_] = x[i2+j,i_] - vc*x[i2+i,i_];
    2114                             }
    2115                         }
    2116                     }
    2117                     return;
    2118                 }
    2119             }
    2120         }
    2121 
    2122 
    2123         /*************************************************************************
    2124         Level 2 subroutine
    2125 
    2126           -- ALGLIB routine --
    2127              15.12.2009
    2128              Bochkanov Sergey
    2129         *************************************************************************/
    2130         private static void rmatrixrighttrsm2(int m,
    2131             int n,
    2132             ref double[,] a,
    2133             int i1,
    2134             int j1,
    2135             bool isupper,
    2136             bool isunit,
    2137             int optype,
    2138             ref double[,] x,
    2139             int i2,
    2140             int j2)
    2141         {
    2142             int i = 0;
    2143             int j = 0;
    2144             double vr = 0;
    2145             double vd = 0;
    2146             int i_ = 0;
    2147             int i1_ = 0;
    2148 
    2149            
    2150             //
    2151             // Special case
    2152             //
    2153             if( n*m==0 )
    2154             {
    2155                 return;
    2156             }
    2157            
    2158             //
    2159             // Try to use "fast" code
    2160             //
    2161             if( ablasf.rmatrixrighttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
    2162             {
    2163                 return;
    2164             }
    2165            
    2166             //
    2167             // General case
    2168             //
    2169             if( isupper )
    2170             {
    2171                
    2172                 //
    2173                 // Upper triangular matrix
    2174                 //
    2175                 if( optype==0 )
    2176                 {
    2177                    
    2178                     //
    2179                     // X*A^(-1)
    2180                     //
    2181                     for(i=0; i<=m-1; i++)
    2182                     {
    2183                         for(j=0; j<=n-1; j++)
    2184                         {
    2185                             if( isunit )
    2186                             {
    2187                                 vd = 1;
    2188                             }
    2189                             else
    2190                             {
    2191                                 vd = a[i1+j,j1+j];
    2192                             }
    2193                             x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
    2194                             if( j<n-1 )
    2195                             {
    2196                                 vr = x[i2+i,j2+j];
    2197                                 i1_ = (j1+j+1) - (j2+j+1);
    2198                                 for(i_=j2+j+1; i_<=j2+n-1;i_++)
    2199                                 {
    2200                                     x[i2+i,i_] = x[i2+i,i_] - vr*a[i1+j,i_+i1_];
    2201                                 }
    2202                             }
    2203                         }
    2204                     }
    2205                     return;
    2206                 }
    2207                 if( optype==1 )
    2208                 {
    2209                    
    2210                     //
    2211                     // X*A^(-T)
    2212                     //
    2213                     for(i=0; i<=m-1; i++)
    2214                     {
    2215                         for(j=n-1; j>=0; j--)
    2216                         {
    2217                             vr = 0;
    2218                             vd = 1;
    2219                             if( j<n-1 )
    2220                             {
    2221                                 i1_ = (j1+j+1)-(j2+j+1);
    2222                                 vr = 0.0;
    2223                                 for(i_=j2+j+1; i_<=j2+n-1;i_++)
    2224                                 {
    2225                                     vr += x[i2+i,i_]*a[i1+j,i_+i1_];
    2226                                 }
    2227                             }
    2228                             if( !isunit )
    2229                             {
    2230                                 vd = a[i1+j,j1+j];
    2231                             }
    2232                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vr)/vd;
    2233                         }
    2234                     }
    2235                     return;
    2236                 }
    2237             }
    2238             else
    2239             {
    2240                
    2241                 //
    2242                 // Lower triangular matrix
    2243                 //
    2244                 if( optype==0 )
    2245                 {
    2246                    
    2247                     //
    2248                     // X*A^(-1)
    2249                     //
    2250                     for(i=0; i<=m-1; i++)
    2251                     {
    2252                         for(j=n-1; j>=0; j--)
    2253                         {
    2254                             if( isunit )
    2255                             {
    2256                                 vd = 1;
    2257                             }
    2258                             else
    2259                             {
    2260                                 vd = a[i1+j,j1+j];
    2261                             }
    2262                             x[i2+i,j2+j] = x[i2+i,j2+j]/vd;
    2263                             if( j>0 )
    2264                             {
    2265                                 vr = x[i2+i,j2+j];
    2266                                 i1_ = (j1) - (j2);
    2267                                 for(i_=j2; i_<=j2+j-1;i_++)
    2268                                 {
    2269                                     x[i2+i,i_] = x[i2+i,i_] - vr*a[i1+j,i_+i1_];
    2270                                 }
    2271                             }
    2272                         }
    2273                     }
    2274                     return;
    2275                 }
    2276                 if( optype==1 )
    2277                 {
    2278                    
    2279                     //
    2280                     // X*A^(-T)
    2281                     //
    2282                     for(i=0; i<=m-1; i++)
    2283                     {
    2284                         for(j=0; j<=n-1; j++)
    2285                         {
    2286                             vr = 0;
    2287                             vd = 1;
    2288                             if( j>0 )
    2289                             {
    2290                                 i1_ = (j1)-(j2);
    2291                                 vr = 0.0;
    2292                                 for(i_=j2; i_<=j2+j-1;i_++)
    2293                                 {
    2294                                     vr += x[i2+i,i_]*a[i1+j,i_+i1_];
    2295                                 }
    2296                             }
    2297                             if( !isunit )
    2298                             {
    2299                                 vd = a[i1+j,j1+j];
    2300                             }
    2301                             x[i2+i,j2+j] = (x[i2+i,j2+j]-vr)/vd;
    2302                         }
    2303                     }
    2304                     return;
    2305                 }
    2306             }
    2307         }
    2308 
    2309 
    2310         /*************************************************************************
    2311         Level 2 subroutine
    2312         *************************************************************************/
    2313         private static void rmatrixlefttrsm2(int m,
    2314             int n,
    2315             ref double[,] a,
    2316             int i1,
    2317             int j1,
    2318             bool isupper,
    2319             bool isunit,
    2320             int optype,
    2321             ref double[,] x,
    2322             int i2,
    2323             int j2)
    2324         {
    2325             int i = 0;
    2326             int j = 0;
    2327             double vr = 0;
    2328             double vd = 0;
    2329             int i_ = 0;
    2330 
    2331            
    2332             //
    2333             // Special case
    2334             //
    2335             if( n*m==0 )
    2336             {
    2337                 return;
    2338             }
    2339            
    2340             //
    2341             // Try fast code
    2342             //
    2343             if( ablasf.rmatrixlefttrsmf(m, n, ref a, i1, j1, isupper, isunit, optype, ref x, i2, j2) )
    2344             {
    2345                 return;
    2346             }
    2347            
    2348             //
    2349             // General case
    2350             //
    2351             if( isupper )
    2352             {
    2353                
    2354                 //
    2355                 // Upper triangular matrix
    2356                 //
    2357                 if( optype==0 )
    2358                 {
    2359                    
    2360                     //
    2361                     // A^(-1)*X
    2362                     //
    2363                     for(i=m-1; i>=0; i--)
    2364                     {
    2365                         for(j=i+1; j<=m-1; j++)
    2366                         {
    2367                             vr = a[i1+i,j1+j];
    2368                             for(i_=j2; i_<=j2+n-1;i_++)
    2369                             {
    2370                                 x[i2+i,i_] = x[i2+i,i_] - vr*x[i2+j,i_];
    2371                             }
    2372                         }
    2373                         if( !isunit )
    2374                         {
    2375                             vd = 1/a[i1+i,j1+i];
    2376                             for(i_=j2; i_<=j2+n-1;i_++)
    2377                             {
    2378                                 x[i2+i,i_] = vd*x[i2+i,i_];
    2379                             }
    2380                         }
    2381                     }
    2382                     return;
    2383                 }
    2384                 if( optype==1 )
    2385                 {
    2386                    
    2387                     //
    2388                     // A^(-T)*X
    2389                     //
    2390                     for(i=0; i<=m-1; i++)
    2391                     {
    2392                         if( isunit )
    2393                         {
    2394                             vd = 1;
    2395                         }
    2396                         else
    2397                         {
    2398                             vd = 1/a[i1+i,j1+i];
    2399                         }
    2400                         for(i_=j2; i_<=j2+n-1;i_++)
    2401                         {
    2402                             x[i2+i,i_] = vd*x[i2+i,i_];
    2403                         }
    2404                         for(j=i+1; j<=m-1; j++)
    2405                         {
    2406                             vr = a[i1+i,j1+j];
    2407                             for(i_=j2; i_<=j2+n-1;i_++)
    2408                             {
    2409                                 x[i2+j,i_] = x[i2+j,i_] - vr*x[i2+i,i_];
    2410                             }
    2411                         }
    2412                     }
    2413                     return;
    2414                 }
    2415             }
    2416             else
    2417             {
    2418                
    2419                 //
    2420                 // Lower triangular matrix
    2421                 //
    2422                 if( optype==0 )
    2423                 {
    2424                    
    2425                     //
    2426                     // A^(-1)*X
    2427                     //
    2428                     for(i=0; i<=m-1; i++)
    2429                     {
    2430                         for(j=0; j<=i-1; j++)
    2431                         {
    2432                             vr = a[i1+i,j1+j];
    2433                             for(i_=j2; i_<=j2+n-1;i_++)
    2434                             {
    2435                                 x[i2+i,i_] = x[i2+i,i_] - vr*x[i2+j,i_];
    2436                             }
    2437                         }
    2438                         if( isunit )
    2439                         {
    2440                             vd = 1;
    2441                         }
    2442                         else
    2443                         {
    2444                             vd = 1/a[i1+j,j1+j];
    2445                         }
    2446                         for(i_=j2; i_<=j2+n-1;i_++)
    2447                         {
    2448                             x[i2+i,i_] = vd*x[i2+i,i_];
    2449                         }
    2450                     }
    2451                     return;
    2452                 }
    2453                 if( optype==1 )
    2454                 {
    2455                    
    2456                     //
    2457                     // A^(-T)*X
    2458                     //
    2459                     for(i=m-1; i>=0; i--)
    2460                     {
    2461                         if( isunit )
    2462                         {
    2463                             vd = 1;
    2464                         }
    2465                         else
    2466                         {
    2467                             vd = 1/a[i1+i,j1+i];
    2468                         }
    2469                         for(i_=j2; i_<=j2+n-1;i_++)
    2470                         {
    2471                             x[i2+i,i_] = vd*x[i2+i,i_];
    2472                         }
    2473                         for(j=i-1; j>=0; j--)
    2474                         {
    2475                             vr = a[i1+i,j1+j];
    2476                             for(i_=j2; i_<=j2+n-1;i_++)
    2477                             {
    2478                                 x[i2+j,i_] = x[i2+j,i_] - vr*x[i2+i,i_];
    2479                             }
    2480                         }
    2481                     }
    2482                     return;
    2483                 }
    2484             }
    2485         }
    2486 
    2487 
    2488         /*************************************************************************
    2489         Level 2 subroutine
    2490         *************************************************************************/
    2491         private static void cmatrixsyrk2(int n,
    2492             int k,
    2493             double alpha,
    2494             ref AP.Complex[,] a,
    2495             int ia,
    2496             int ja,
    2497             int optypea,
    2498             double beta,
    2499             ref AP.Complex[,] c,
    2500             int ic,
    2501             int jc,
    2502             bool isupper)
    2503         {
    2504             int i = 0;
    2505             int j = 0;
    2506             int j1 = 0;
    2507             int j2 = 0;
    2508             AP.Complex v = 0;
    2509             int i_ = 0;
    2510             int i1_ = 0;
    2511 
    2512            
    2513             //
    2514             // Fast exit (nothing to be done)
    2515             //
    2516             if( ((double)(alpha)==(double)(0) | k==0) & (double)(beta)==(double)(1) )
    2517             {
    2518                 return;
    2519             }
    2520            
    2521             //
    2522             // Try to call fast SYRK
    2523             //
    2524             if( ablasf.cmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper) )
    2525             {
    2526                 return;
    2527             }
    2528            
    2529             //
    2530             // SYRK
    2531             //
    2532             if( optypea==0 )
    2533             {
    2534                
    2535                 //
    2536                 // C=alpha*A*A^H+beta*C
    2537                 //
    2538                 for(i=0; i<=n-1; i++)
    2539                 {
    2540                     if( isupper )
    2541                     {
    2542                         j1 = i;
    2543                         j2 = n-1;
    2544                     }
    2545                     else
    2546                     {
    2547                         j1 = 0;
    2548                         j2 = i;
    2549                     }
    2550                     for(j=j1; j<=j2; j++)
    2551                     {
    2552                         if( (double)(alpha)!=(double)(0) & k>0 )
    2553                         {
    2554                             v = 0.0;
    2555                             for(i_=ja; i_<=ja+k-1;i_++)
    2556                             {
    2557                                 v += a[ia+i,i_]*AP.Math.Conj(a[ia+j,i_]);
    2558                             }
    2559                         }
    2560                         else
    2561                         {
    2562                             v = 0;
    2563                         }
    2564                         if( (double)(beta)==(double)(0) )
    2565                         {
    2566                             c[ic+i,jc+j] = alpha*v;
    2567                         }
    2568                         else
    2569                         {
    2570                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    2571                         }
    2572                     }
    2573                 }
    2574                 return;
    2575             }
    2576             else
    2577             {
    2578                
    2579                 //
    2580                 // C=alpha*A^H*A+beta*C
    2581                 //
    2582                 for(i=0; i<=n-1; i++)
    2583                 {
    2584                     if( isupper )
    2585                     {
    2586                         j1 = i;
    2587                         j2 = n-1;
    2588                     }
    2589                     else
    2590                     {
    2591                         j1 = 0;
    2592                         j2 = i;
    2593                     }
    2594                     if( (double)(beta)==(double)(0) )
    2595                     {
    2596                         for(j=j1; j<=j2; j++)
    2597                         {
    2598                             c[ic+i,jc+j] = 0;
    2599                         }
    2600                     }
    2601                     else
    2602                     {
    2603                         for(i_=jc+j1; i_<=jc+j2;i_++)
    2604                         {
    2605                             c[ic+i,i_] = beta*c[ic+i,i_];
    2606                         }
    2607                     }
    2608                 }
    2609                 for(i=0; i<=k-1; i++)
    2610                 {
    2611                     for(j=0; j<=n-1; j++)
    2612                     {
    2613                         if( isupper )
    2614                         {
    2615                             j1 = j;
    2616                             j2 = n-1;
    2617                         }
    2618                         else
    2619                         {
    2620                             j1 = 0;
    2621                             j2 = j;
    2622                         }
    2623                         v = alpha*AP.Math.Conj(a[ia+i,ja+j]);
    2624                         i1_ = (ja+j1) - (jc+j1);
    2625                         for(i_=jc+j1; i_<=jc+j2;i_++)
    2626                         {
    2627                             c[ic+j,i_] = c[ic+j,i_] + v*a[ia+i,i_+i1_];
    2628                         }
    2629                     }
    2630                 }
    2631                 return;
    2632             }
    2633         }
    2634 
    2635 
    2636         /*************************************************************************
    2637         Level 2 subrotuine
    2638         *************************************************************************/
    2639         private static void rmatrixsyrk2(int n,
    2640             int k,
    2641             double alpha,
    2642             ref double[,] a,
    2643             int ia,
    2644             int ja,
    2645             int optypea,
    2646             double beta,
    2647             ref double[,] c,
    2648             int ic,
    2649             int jc,
    2650             bool isupper)
    2651         {
    2652             int i = 0;
    2653             int j = 0;
    2654             int j1 = 0;
    2655             int j2 = 0;
    2656             double v = 0;
    2657             int i_ = 0;
    2658             int i1_ = 0;
    2659 
    2660            
    2661             //
    2662             // Fast exit (nothing to be done)
    2663             //
    2664             if( ((double)(alpha)==(double)(0) | k==0) & (double)(beta)==(double)(1) )
    2665             {
    2666                 return;
    2667             }
    2668            
    2669             //
    2670             // Try to call fast SYRK
    2671             //
    2672             if( ablasf.rmatrixsyrkf(n, k, alpha, ref a, ia, ja, optypea, beta, ref c, ic, jc, isupper) )
    2673             {
    2674                 return;
    2675             }
    2676            
    2677             //
    2678             // SYRK
    2679             //
    2680             if( optypea==0 )
    2681             {
    2682                
    2683                 //
    2684                 // C=alpha*A*A^H+beta*C
    2685                 //
    2686                 for(i=0; i<=n-1; i++)
    2687                 {
    2688                     if( isupper )
    2689                     {
    2690                         j1 = i;
    2691                         j2 = n-1;
    2692                     }
    2693                     else
    2694                     {
    2695                         j1 = 0;
    2696                         j2 = i;
    2697                     }
    2698                     for(j=j1; j<=j2; j++)
    2699                     {
    2700                         if( (double)(alpha)!=(double)(0) & k>0 )
    2701                         {
    2702                             v = 0.0;
    2703                             for(i_=ja; i_<=ja+k-1;i_++)
    2704                             {
    2705                                 v += a[ia+i,i_]*a[ia+j,i_];
    2706                             }
    2707                         }
    2708                         else
    2709                         {
    2710                             v = 0;
    2711                         }
    2712                         if( (double)(beta)==(double)(0) )
    2713                         {
    2714                             c[ic+i,jc+j] = alpha*v;
    2715                         }
    2716                         else
    2717                         {
    2718                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    2719                         }
    2720                     }
    2721                 }
    2722                 return;
    2723             }
    2724             else
    2725             {
    2726                
    2727                 //
    2728                 // C=alpha*A^H*A+beta*C
    2729                 //
    2730                 for(i=0; i<=n-1; i++)
    2731                 {
    2732                     if( isupper )
    2733                     {
    2734                         j1 = i;
    2735                         j2 = n-1;
    2736                     }
    2737                     else
    2738                     {
    2739                         j1 = 0;
    2740                         j2 = i;
    2741                     }
    2742                     if( (double)(beta)==(double)(0) )
    2743                     {
    2744                         for(j=j1; j<=j2; j++)
    2745                         {
    2746                             c[ic+i,jc+j] = 0;
    2747                         }
    2748                     }
    2749                     else
    2750                     {
    2751                         for(i_=jc+j1; i_<=jc+j2;i_++)
    2752                         {
    2753                             c[ic+i,i_] = beta*c[ic+i,i_];
    2754                         }
    2755                     }
    2756                 }
    2757                 for(i=0; i<=k-1; i++)
    2758                 {
    2759                     for(j=0; j<=n-1; j++)
    2760                     {
    2761                         if( isupper )
    2762                         {
    2763                             j1 = j;
    2764                             j2 = n-1;
    2765                         }
    2766                         else
    2767                         {
    2768                             j1 = 0;
    2769                             j2 = j;
    2770                         }
    2771                         v = alpha*a[ia+i,ja+j];
    2772                         i1_ = (ja+j1) - (jc+j1);
    2773                         for(i_=jc+j1; i_<=jc+j2;i_++)
    2774                         {
    2775                             c[ic+j,i_] = c[ic+j,i_] + v*a[ia+i,i_+i1_];
    2776                         }
    2777                     }
    2778                 }
    2779                 return;
    2780             }
    2781         }
    2782 
    2783 
    2784         /*************************************************************************
    2785         GEMM kernel
    2786 
    2787           -- ALGLIB routine --
    2788              16.12.2009
    2789              Bochkanov Sergey
    2790         *************************************************************************/
    2791         private static void cmatrixgemmk(int m,
    2792             int n,
    2793             int k,
    2794             AP.Complex alpha,
    2795             ref AP.Complex[,] a,
    2796             int ia,
    2797             int ja,
    2798             int optypea,
    2799             ref AP.Complex[,] b,
    2800             int ib,
    2801             int jb,
    2802             int optypeb,
    2803             AP.Complex beta,
    2804             ref AP.Complex[,] c,
    2805             int ic,
    2806             int jc)
    2807         {
    2808             int i = 0;
    2809             int j = 0;
    2810             AP.Complex v = 0;
    2811             int i_ = 0;
    2812             int i1_ = 0;
    2813 
    2814            
    2815             //
    2816             // Special case
    2817             //
    2818             if( m*n==0 )
    2819             {
    2820                 return;
    2821             }
    2822            
    2823             //
    2824             // Try optimized code
    2825             //
    2826             if( ablasf.cmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc) )
    2827             {
    2828                 return;
    2829             }
    2830            
    2831             //
    2832             // Another special case
    2833             //
    2834             if( k==0 )
    2835             {
    2836                 if( beta!=0 )
    2837                 {
    2838                     for(i=0; i<=m-1; i++)
    2839                     {
    2840                         for(j=0; j<=n-1; j++)
    2841                         {
    2842                             c[ic+i,jc+j] = beta*c[ic+i,jc+j];
    2843                         }
    2844                     }
    2845                 }
    2846                 else
    2847                 {
    2848                     for(i=0; i<=m-1; i++)
    2849                     {
    2850                         for(j=0; j<=n-1; j++)
    2851                         {
    2852                             c[ic+i,jc+j] = 0;
    2853                         }
    2854                     }
    2855                 }
    2856                 return;
    2857             }
    2858            
    2859             //
    2860             // General case
    2861             //
    2862             if( optypea==0 & optypeb!=0 )
    2863             {
    2864                
    2865                 //
    2866                 // A*B'
    2867                 //
    2868                 for(i=0; i<=m-1; i++)
    2869                 {
    2870                     for(j=0; j<=n-1; j++)
    2871                     {
    2872                         if( k==0 | alpha==0 )
    2873                         {
    2874                             v = 0;
    2875                         }
    2876                         else
    2877                         {
    2878                             if( optypeb==1 )
    2879                             {
    2880                                 i1_ = (jb)-(ja);
    2881                                 v = 0.0;
    2882                                 for(i_=ja; i_<=ja+k-1;i_++)
    2883                                 {
    2884                                     v += a[ia+i,i_]*b[ib+j,i_+i1_];
    2885                                 }
    2886                             }
    2887                             else
    2888                             {
    2889                                 i1_ = (jb)-(ja);
    2890                                 v = 0.0;
    2891                                 for(i_=ja; i_<=ja+k-1;i_++)
    2892                                 {
    2893                                     v += a[ia+i,i_]*AP.Math.Conj(b[ib+j,i_+i1_]);
    2894                                 }
    2895                             }
    2896                         }
    2897                         if( beta==0 )
    2898                         {
    2899                             c[ic+i,jc+j] = alpha*v;
    2900                         }
    2901                         else
    2902                         {
    2903                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    2904                         }
    2905                     }
    2906                 }
    2907                 return;
    2908             }
    2909             if( optypea==0 & optypeb==0 )
    2910             {
    2911                
    2912                 //
    2913                 // A*B
    2914                 //
    2915                 for(i=0; i<=m-1; i++)
    2916                 {
    2917                     if( beta!=0 )
    2918                     {
    2919                         for(i_=jc; i_<=jc+n-1;i_++)
    2920                         {
    2921                             c[ic+i,i_] = beta*c[ic+i,i_];
    2922                         }
    2923                     }
    2924                     else
    2925                     {
    2926                         for(j=0; j<=n-1; j++)
    2927                         {
    2928                             c[ic+i,jc+j] = 0;
    2929                         }
    2930                     }
    2931                     if( alpha!=0 )
    2932                     {
    2933                         for(j=0; j<=k-1; j++)
    2934                         {
    2935                             v = alpha*a[ia+i,ja+j];
    2936                             i1_ = (jb) - (jc);
    2937                             for(i_=jc; i_<=jc+n-1;i_++)
    2938                             {
    2939                                 c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
    2940                             }
    2941                         }
    2942                     }
    2943                 }
    2944                 return;
    2945             }
    2946             if( optypea!=0 & optypeb!=0 )
    2947             {
    2948                
    2949                 //
    2950                 // A'*B'
    2951                 //
    2952                 for(i=0; i<=m-1; i++)
    2953                 {
    2954                     for(j=0; j<=n-1; j++)
    2955                     {
    2956                         if( alpha==0 )
    2957                         {
    2958                             v = 0;
    2959                         }
    2960                         else
    2961                         {
    2962                             if( optypea==1 )
    2963                             {
    2964                                 if( optypeb==1 )
    2965                                 {
    2966                                     i1_ = (jb)-(ia);
    2967                                     v = 0.0;
    2968                                     for(i_=ia; i_<=ia+k-1;i_++)
    2969                                     {
    2970                                         v += a[i_,ja+i]*b[ib+j,i_+i1_];
    2971                                     }
    2972                                 }
    2973                                 else
    2974                                 {
    2975                                     i1_ = (jb)-(ia);
    2976                                     v = 0.0;
    2977                                     for(i_=ia; i_<=ia+k-1;i_++)
    2978                                     {
    2979                                         v += a[i_,ja+i]*AP.Math.Conj(b[ib+j,i_+i1_]);
    2980                                     }
    2981                                 }
    2982                             }
    2983                             else
    2984                             {
    2985                                 if( optypeb==1 )
    2986                                 {
    2987                                     i1_ = (jb)-(ia);
    2988                                     v = 0.0;
    2989                                     for(i_=ia; i_<=ia+k-1;i_++)
    2990                                     {
    2991                                         v += AP.Math.Conj(a[i_,ja+i])*b[ib+j,i_+i1_];
    2992                                     }
    2993                                 }
    2994                                 else
    2995                                 {
    2996                                     i1_ = (jb)-(ia);
    2997                                     v = 0.0;
    2998                                     for(i_=ia; i_<=ia+k-1;i_++)
    2999                                     {
    3000                                         v += AP.Math.Conj(a[i_,ja+i])*AP.Math.Conj(b[ib+j,i_+i1_]);
    3001                                     }
    3002                                 }
    3003                             }
    3004                         }
    3005                         if( beta==0 )
    3006                         {
    3007                             c[ic+i,jc+j] = alpha*v;
    3008                         }
    3009                         else
    3010                         {
    3011                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    3012                         }
    3013                     }
    3014                 }
    3015                 return;
    3016             }
    3017             if( optypea!=0 & optypeb==0 )
    3018             {
    3019                
    3020                 //
    3021                 // A'*B
    3022                 //
    3023                 if( beta==0 )
    3024                 {
    3025                     for(i=0; i<=m-1; i++)
    3026                     {
    3027                         for(j=0; j<=n-1; j++)
    3028                         {
    3029                             c[ic+i,jc+j] = 0;
    3030                         }
    3031                     }
    3032                 }
    3033                 else
    3034                 {
    3035                     for(i=0; i<=m-1; i++)
    3036                     {
    3037                         for(i_=jc; i_<=jc+n-1;i_++)
    3038                         {
    3039                             c[ic+i,i_] = beta*c[ic+i,i_];
    3040                         }
    3041                     }
    3042                 }
    3043                 if( alpha!=0 )
    3044                 {
    3045                     for(j=0; j<=k-1; j++)
    3046                     {
    3047                         for(i=0; i<=m-1; i++)
    3048                         {
    3049                             if( optypea==1 )
    3050                             {
    3051                                 v = alpha*a[ia+j,ja+i];
    3052                             }
    3053                             else
    3054                             {
    3055                                 v = alpha*AP.Math.Conj(a[ia+j,ja+i]);
    3056                             }
    3057                             i1_ = (jb) - (jc);
    3058                             for(i_=jc; i_<=jc+n-1;i_++)
    3059                             {
    3060                                 c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
    3061                             }
    3062                         }
    3063                     }
    3064                 }
    3065                 return;
    3066             }
    3067         }
    3068 
    3069 
    3070         /*************************************************************************
    3071         GEMM kernel
    3072 
    3073           -- ALGLIB routine --
    3074              16.12.2009
    3075              Bochkanov Sergey
    3076         *************************************************************************/
    3077         private static void rmatrixgemmk(int m,
    3078             int n,
    3079             int k,
    3080             double alpha,
    3081             ref double[,] a,
    3082             int ia,
    3083             int ja,
    3084             int optypea,
    3085             ref double[,] b,
    3086             int ib,
    3087             int jb,
    3088             int optypeb,
    3089             double beta,
    3090             ref double[,] c,
    3091             int ic,
    3092             int jc)
    3093         {
    3094             int i = 0;
    3095             int j = 0;
    3096             double v = 0;
    3097             int i_ = 0;
    3098             int i1_ = 0;
    3099 
    3100            
    3101             //
    3102             // if matrix size is zero
    3103             //
    3104             if( m*n==0 )
    3105             {
    3106                 return;
    3107             }
    3108            
    3109             //
    3110             // Try optimized code
    3111             //
    3112             if( ablasf.rmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc) )
    3113             {
    3114                 return;
    3115             }
    3116            
    3117             //
    3118             // if K=0, then C=Beta*C
    3119             //
    3120             if( k==0 )
    3121             {
    3122                 if( (double)(beta)!=(double)(1) )
    3123                 {
    3124                     if( (double)(beta)!=(double)(0) )
    3125                     {
    3126                         for(i=0; i<=m-1; i++)
    3127                         {
    3128                             for(j=0; j<=n-1; j++)
    3129                             {
    3130                                 c[ic+i,jc+j] = beta*c[ic+i,jc+j];
    3131                             }
    3132                         }
    3133                     }
    3134                     else
    3135                     {
    3136                         for(i=0; i<=m-1; i++)
    3137                         {
    3138                             for(j=0; j<=n-1; j++)
    3139                             {
    3140                                 c[ic+i,jc+j] = 0;
    3141                             }
    3142                         }
    3143                     }
    3144                 }
    3145                 return;
    3146             }
    3147            
    3148             //
    3149             // General case
    3150             //
    3151             if( optypea==0 & optypeb!=0 )
    3152             {
    3153                
    3154                 //
    3155                 // A*B'
    3156                 //
    3157                 for(i=0; i<=m-1; i++)
    3158                 {
    3159                     for(j=0; j<=n-1; j++)
    3160                     {
    3161                         if( k==0 | (double)(alpha)==(double)(0) )
    3162                         {
    3163                             v = 0;
    3164                         }
    3165                         else
    3166                         {
    3167                             i1_ = (jb)-(ja);
    3168                             v = 0.0;
    3169                             for(i_=ja; i_<=ja+k-1;i_++)
    3170                             {
    3171                                 v += a[ia+i,i_]*b[ib+j,i_+i1_];
    3172                             }
    3173                         }
    3174                         if( (double)(beta)==(double)(0) )
    3175                         {
    3176                             c[ic+i,jc+j] = alpha*v;
    3177                         }
    3178                         else
    3179                         {
    3180                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    3181                         }
    3182                     }
    3183                 }
    3184                 return;
    3185             }
    3186             if( optypea==0 & optypeb==0 )
    3187             {
    3188                
    3189                 //
    3190                 // A*B
    3191                 //
    3192                 for(i=0; i<=m-1; i++)
    3193                 {
    3194                     if( (double)(beta)!=(double)(0) )
    3195                     {
    3196                         for(i_=jc; i_<=jc+n-1;i_++)
    3197                         {
    3198                             c[ic+i,i_] = beta*c[ic+i,i_];
    3199                         }
    3200                     }
    3201                     else
    3202                     {
    3203                         for(j=0; j<=n-1; j++)
    3204                         {
    3205                             c[ic+i,jc+j] = 0;
    3206                         }
    3207                     }
    3208                     if( (double)(alpha)!=(double)(0) )
    3209                     {
    3210                         for(j=0; j<=k-1; j++)
    3211                         {
    3212                             v = alpha*a[ia+i,ja+j];
    3213                             i1_ = (jb) - (jc);
    3214                             for(i_=jc; i_<=jc+n-1;i_++)
    3215                             {
    3216                                 c[ic+i,i_] = c[ic+i,i_] + v*b[ib+j,i_+i1_];
    3217                             }
    3218                         }
    3219                     }
    3220                 }
    3221                 return;
    3222             }
    3223             if( optypea!=0 & optypeb!=0 )
    3224             {
    3225                
    3226                 //
    3227                 // A'*B'
    3228                 //
    3229                 for(i=0; i<=m-1; i++)
    3230                 {
    3231                     for(j=0; j<=n-1; j++)
    3232                     {
    3233                         if( (double)(alpha)==(double)(0) )
    3234                         {
    3235                             v = 0;
    3236                         }
    3237                         else
    3238                         {
    3239                             i1_ = (jb)-(ia);
    3240                             v = 0.0;
    3241                             for(i_=ia; i_<=ia+k-1;i_++)
    3242                             {
    3243                                 v += a[i_,ja+i]*b[ib+j,i_+i1_];
    3244                             }
    3245                         }
    3246                         if( (double)(beta)==(double)(0) )
    3247                         {
    3248                             c[ic+i,jc+j] = alpha*v;
    3249                         }
    3250                         else
    3251                         {
    3252                             c[ic+i,jc+j] = beta*c[ic+i,jc+j]+alpha*v;
    3253                         }
    3254                     }
    3255                 }
    3256                 return;
    3257             }
    3258             if( optypea!=0 & optypeb==0 )
    3259             {
    3260                
    3261                 //
    3262                 // A'*B
    3263