Free cookie consent management tool by TermsFeed Policy Generator

source: tags/3.3.2/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.5.0/ALGLIB-2.5.0/ablas.cs @ 13398

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

Sorted usings and removed unused usings in entire solution (#1094)
File size: 82.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2009-2010, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
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_];
1528                }
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_];
1548                }
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_]);
1572                }
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_];
1605                }
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_];
1625                }
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_]);
1649                }
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_];
1917                }
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_];
1937                }
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_];
1970                }
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_];
1990                }
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_];
2458                }
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_]);
2464                }
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                  }
2526                }
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                  }
2540                }
2541              }
2542            }
2543            if (beta == 0) {
2544              c[ic + i, jc + j] = alpha * v;
2545            } else {
2546              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
2547            }
2548          }
2549        }
2550        return;
2551      }
2552      if (optypea != 0 & optypeb == 0) {
2553
2554        //
2555        // A'*B
2556        //
2557        if (beta == 0) {
2558          for (i = 0; i <= m - 1; i++) {
2559            for (j = 0; j <= n - 1; j++) {
2560              c[ic + i, jc + j] = 0;
2561            }
2562          }
2563        } else {
2564          for (i = 0; i <= m - 1; i++) {
2565            for (i_ = jc; i_ <= jc + n - 1; i_++) {
2566              c[ic + i, i_] = beta * c[ic + i, i_];
2567            }
2568          }
2569        }
2570        if (alpha != 0) {
2571          for (j = 0; j <= k - 1; j++) {
2572            for (i = 0; i <= m - 1; i++) {
2573              if (optypea == 1) {
2574                v = alpha * a[ia + j, ja + i];
2575              } else {
2576                v = alpha * AP.Math.Conj(a[ia + j, ja + i]);
2577              }
2578              i1_ = (jb) - (jc);
2579              for (i_ = jc; i_ <= jc + n - 1; i_++) {
2580                c[ic + i, i_] = c[ic + i, i_] + v * b[ib + j, i_ + i1_];
2581              }
2582            }
2583          }
2584        }
2585        return;
2586      }
2587    }
2588
2589
2590    /*************************************************************************
2591    GEMM kernel
2592
2593      -- ALGLIB routine --
2594         16.12.2009
2595         Bochkanov Sergey
2596    *************************************************************************/
2597    private static void rmatrixgemmk(int m,
2598        int n,
2599        int k,
2600        double alpha,
2601        ref double[,] a,
2602        int ia,
2603        int ja,
2604        int optypea,
2605        ref double[,] b,
2606        int ib,
2607        int jb,
2608        int optypeb,
2609        double beta,
2610        ref double[,] c,
2611        int ic,
2612        int jc) {
2613      int i = 0;
2614      int j = 0;
2615      double v = 0;
2616      int i_ = 0;
2617      int i1_ = 0;
2618
2619
2620      //
2621      // if matrix size is zero
2622      //
2623      if (m * n == 0) {
2624        return;
2625      }
2626
2627      //
2628      // Try optimized code
2629      //
2630      if (ablasf.rmatrixgemmf(m, n, k, alpha, ref a, ia, ja, optypea, ref b, ib, jb, optypeb, beta, ref c, ic, jc)) {
2631        return;
2632      }
2633
2634      //
2635      // if K=0, then C=Beta*C
2636      //
2637      if (k == 0) {
2638        if ((double)(beta) != (double)(1)) {
2639          if ((double)(beta) != (double)(0)) {
2640            for (i = 0; i <= m - 1; i++) {
2641              for (j = 0; j <= n - 1; j++) {
2642                c[ic + i, jc + j] = beta * c[ic + i, jc + j];
2643              }
2644            }
2645          } else {
2646            for (i = 0; i <= m - 1; i++) {
2647              for (j = 0; j <= n - 1; j++) {
2648                c[ic + i, jc + j] = 0;
2649              }
2650            }
2651          }
2652        }
2653        return;
2654      }
2655
2656      //
2657      // General case
2658      //
2659      if (optypea == 0 & optypeb != 0) {
2660
2661        //
2662        // A*B'
2663        //
2664        for (i = 0; i <= m - 1; i++) {
2665          for (j = 0; j <= n - 1; j++) {
2666            if (k == 0 | (double)(alpha) == (double)(0)) {
2667              v = 0;
2668            } else {
2669              i1_ = (jb) - (ja);
2670              v = 0.0;
2671              for (i_ = ja; i_ <= ja + k - 1; i_++) {
2672                v += a[ia + i, i_] * b[ib + j, i_ + i1_];
2673              }
2674            }
2675            if ((double)(beta) == (double)(0)) {
2676              c[ic + i, jc + j] = alpha * v;
2677            } else {
2678              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
2679            }
2680          }
2681        }
2682        return;
2683      }
2684      if (optypea == 0 & optypeb == 0) {
2685
2686        //
2687        // A*B
2688        //
2689        for (i = 0; i <= m - 1; i++) {
2690          if ((double)(beta) != (double)(0)) {
2691            for (i_ = jc; i_ <= jc + n - 1; i_++) {
2692              c[ic + i, i_] = beta * c[ic + i, i_];
2693            }
2694          } else {
2695            for (j = 0; j <= n - 1; j++) {
2696              c[ic + i, jc + j] = 0;
2697            }
2698          }
2699          if ((double)(alpha) != (double)(0)) {
2700            for (j = 0; j <= k - 1; j++) {
2701              v = alpha * a[ia + i, ja + j];
2702              i1_ = (jb) - (jc);
2703              for (i_ = jc; i_ <= jc + n - 1; i_++) {
2704                c[ic + i, i_] = c[ic + i, i_] + v * b[ib + j, i_ + i1_];
2705              }
2706            }
2707          }
2708        }
2709        return;
2710      }
2711      if (optypea != 0 & optypeb != 0) {
2712
2713        //
2714        // A'*B'
2715        //
2716        for (i = 0; i <= m - 1; i++) {
2717          for (j = 0; j <= n - 1; j++) {
2718            if ((double)(alpha) == (double)(0)) {
2719              v = 0;
2720            } else {
2721              i1_ = (jb) - (ia);
2722              v = 0.0;
2723              for (i_ = ia; i_ <= ia + k - 1; i_++) {
2724                v += a[i_, ja + i] * b[ib + j, i_ + i1_];
2725              }
2726            }
2727            if ((double)(beta) == (double)(0)) {
2728              c[ic + i, jc + j] = alpha * v;
2729            } else {
2730              c[ic + i, jc + j] = beta * c[ic + i, jc + j] + alpha * v;
2731            }
2732          }
2733        }
2734        return;
2735      }
2736      if (optypea != 0 & optypeb == 0) {
2737
2738        //
2739        // A'*B
2740        //
2741        if ((double)(beta) == (double)(0)) {
2742          for (i = 0; i <= m - 1; i++) {
2743            for (j = 0; j <= n - 1; j++) {
2744              c[ic + i, jc + j] = 0;
2745            }
2746          }
2747        } else {
2748          for (i = 0; i <= m - 1; i++) {
2749            for (i_ = jc; i_ <= jc + n - 1; i_++) {
2750              c[ic + i, i_] = beta * c[ic + i, i_];
2751            }
2752          }
2753        }
2754        if ((double)(alpha) != (double)(0)) {
2755          for (j = 0; j <= k - 1; j++) {
2756            for (i = 0; i <= m - 1; i++) {
2757              v = alpha * a[ia + j, ja + i];
2758              i1_ = (jb) - (jc);
2759              for (i_ = jc; i_ <= jc + n - 1; i_++) {
2760                c[ic + i, i_] = c[ic + i, i_] + v * b[ib + j, i_ + i1_];
2761              }
2762            }
2763          }
2764        }
2765        return;
2766      }
2767    }
2768  }
2769}
Note: See TracBrowser for help on using the repository browser.