This subroutine calculates op(A^-1)*X where: * X is MxN general matrix * A is MxM upper/lower triangular/unitriangular matrix * "op" may be identity transformation, transposition, conjugate transposition Multiplication result replaces X. Cache-oblivious algorithm is used. INPUT PARAMETERS N - matrix size, N>=0 M - matrix size, N>=0 A - matrix, actial matrix is stored
| 862 | Bochkanov Sergey |
| 863 | *************************************************************************/ |
| 864 | void cmatrixlefttrsm(int m, |
| 865 | int n, |
| 866 | const ap::complex_2d_array& a, |
| 867 | int i1, |
| 868 | int j1, |
| 869 | bool isupper, |
| 870 | bool isunit, |
| 871 | int optype, |
| 872 | ap::complex_2d_array& x, |
| 873 | int i2, |
| 874 | int j2) |
| 875 | { |
| 876 | int s1; |
| 877 | int s2; |
| 878 | int bs; |
| 879 | |
| 880 | bs = ablascomplexblocksize(a); |
| 881 | if( m<=bs&&n<=bs ) |
| 882 | { |
| 883 | cmatrixlefttrsm2(m, n, a, i1, j1, isupper, isunit, optype, x, i2, j2); |
| 884 | return; |
| 885 | } |
| 886 | if( n>=m ) |
| 887 | { |
| 888 | |
| 889 | // |
| 890 | // Split X: op(A)^-1*X = op(A)^-1*(X1 X2) |
| 891 | // |
| 892 | ablascomplexsplitlength(x, n, s1, s2); |
| 893 | cmatrixlefttrsm(m, s1, a, i1, j1, isupper, isunit, optype, x, i2, j2); |
| 894 | cmatrixlefttrsm(m, s2, a, i1, j1, isupper, isunit, optype, x, i2, j2+s1); |
| 895 | } |
| 896 | else |
| 897 | { |
| 898 | |
| 899 | // |
| 900 | // Split A |
| 901 | // |
| 902 | ablascomplexsplitlength(a, m, s1, s2); |
| 903 | if( isupper&&optype==0 ) |
| 904 | { |
| 905 | |
| 906 | // |
| 907 | // (A1 A12)-1 ( X1 ) |
| 908 | // A^-1*X* = ( ) *( ) |
| 909 | // ( A2) ( X2 ) |
| 910 | // |
| 911 | cmatrixlefttrsm(s2, n, a, i1+s1, j1+s1, isupper, isunit, optype, x, i2+s1, j2); |
| 912 | cmatrixgemm(s1, n, s2, -1.0, a, i1, j1+s1, 0, x, i2+s1, j2, 0, 1.0, x, i2, j2); |
| 913 | cmatrixlefttrsm(s1, n, a, i1, j1, isupper, isunit, optype, x, i2, j2); |
| 914 | return; |
| 915 | } |
| 916 | if( isupper&&optype!=0 ) |
| 917 | { |
| 918 | |
| 919 | // |
| 920 | // (A1' )-1 ( X1 ) |
| 921 | // A^-1*X = ( ) *( ) |
nothing calls this directly
no test coverage detected