Mercurial > hg > octave-lyh

comparison liboctave/fMatrix.cc @ 9661:afcf852256d2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
optimize / and '\ for triangular matrices
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 23 Sep 2009 10:00:16 +0200
parents 3429c956de6f
children 0d3b248f4ab6
comparison
equal deleted inserted replaced
9660:0256e187d13b 9661:afcf852256d2
1520 } 1520 }
1521 1521
1522 FloatMatrix 1522 FloatMatrix
1523 FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, 1523 FloatMatrix::utsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info,
1524 float& rcon, solve_singularity_handler sing_handler, 1524 float& rcon, solve_singularity_handler sing_handler,
1525 bool calc_cond) const 1525 bool calc_cond, blas_trans_type transt) const
1526 { 1526 {
1527 FloatMatrix retval; 1527 FloatMatrix retval;
1528 1528
1529 octave_idx_type nr = rows (); 1529 octave_idx_type nr = rows ();
1530 octave_idx_type nc = cols (); 1530 octave_idx_type nc = cols ();
1596 { 1596 {
1597 retval = b; 1597 retval = b;
1598 float *result = retval.fortran_vec (); 1598 float *result = retval.fortran_vec ();
1599 1599
1600 char uplo = 'U'; 1600 char uplo = 'U';
1601 char trans = 'N'; 1601 char trans = get_blas_char (transt);
1602 char dia = 'N'; 1602 char dia = 'N';
1603 1603
1604 F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), 1604 F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1605 F77_CONST_CHAR_ARG2 (&trans, 1), 1605 F77_CONST_CHAR_ARG2 (&trans, 1),
1606 F77_CONST_CHAR_ARG2 (&dia, 1), 1606 F77_CONST_CHAR_ARG2 (&dia, 1),
1620 } 1620 }
1621 1621
1622 FloatMatrix 1622 FloatMatrix
1623 FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, 1623 FloatMatrix::ltsolve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info,
1624 float& rcon, solve_singularity_handler sing_handler, 1624 float& rcon, solve_singularity_handler sing_handler,
1625 bool calc_cond) const 1625 bool calc_cond, blas_trans_type transt) const
1626 { 1626 {
1627 FloatMatrix retval; 1627 FloatMatrix retval;
1628 1628
1629 octave_idx_type nr = rows (); 1629 octave_idx_type nr = rows ();
1630 octave_idx_type nc = cols (); 1630 octave_idx_type nc = cols ();
1696 { 1696 {
1697 retval = b; 1697 retval = b;
1698 float *result = retval.fortran_vec (); 1698 float *result = retval.fortran_vec ();
1699 1699
1700 char uplo = 'L'; 1700 char uplo = 'L';
1701 char trans = 'N'; 1701 char trans = get_blas_char (transt);
1702 char dia = 'N'; 1702 char dia = 'N';
1703 1703
1704 F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1), 1704 F77_XFCN (strtrs, STRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1705 F77_CONST_CHAR_ARG2 (&trans, 1), 1705 F77_CONST_CHAR_ARG2 (&trans, 1),
1706 F77_CONST_CHAR_ARG2 (&dia, 1), 1706 F77_CONST_CHAR_ARG2 (&dia, 1),
1917 } 1917 }
1918 1918
1919 FloatMatrix 1919 FloatMatrix
1920 FloatMatrix::solve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info, 1920 FloatMatrix::solve (MatrixType &mattype, const FloatMatrix& b, octave_idx_type& info,
1921 float& rcon, solve_singularity_handler sing_handler, 1921 float& rcon, solve_singularity_handler sing_handler,
1922 bool singular_fallback) const 1922 bool singular_fallback, blas_trans_type transt) const
1923 { 1923 {
1924 FloatMatrix retval; 1924 FloatMatrix retval;
1925 int typ = mattype.type (); 1925 int typ = mattype.type ();
1926 1926
1927 if (typ == MatrixType::Unknown) 1927 if (typ == MatrixType::Unknown)
1928 typ = mattype.type (*this); 1928 typ = mattype.type (*this);
1929 1929
1930 // Only calculate the condition number for LU/Cholesky 1930 // Only calculate the condition number for LU/Cholesky
1931 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) 1931 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1932 retval = utsolve (mattype, b, info, rcon, sing_handler, false); 1932 retval = utsolve (mattype, b, info, rcon, sing_handler, false, transt);
1933 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) 1933 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1934 retval = ltsolve (mattype, b, info, rcon, sing_handler, false); 1934 retval = ltsolve (mattype, b, info, rcon, sing_handler, false, transt);
1935 else if (transt == blas_trans || transt == blas_conj_trans)
1936 return transpose ().solve (mattype, b, info, rcon, sing_handler, singular_fallback);
1935 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 1937 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1936 retval = fsolve (mattype, b, info, rcon, sing_handler, true); 1938 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
1937 else if (typ != MatrixType::Rectangular) 1939 else if (typ != MatrixType::Rectangular)
1938 { 1940 {
1939 (*current_liboctave_error_handler) ("unknown matrix type"); 1941 (*current_liboctave_error_handler) ("unknown matrix type");
1974 } 1976 }
1975 1977
1976 FloatComplexMatrix 1978 FloatComplexMatrix
1977 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info, 1979 FloatMatrix::solve (MatrixType &typ, const FloatComplexMatrix& b, octave_idx_type& info,
1978 float& rcon, solve_singularity_handler sing_handler, 1980 float& rcon, solve_singularity_handler sing_handler,
1979 bool singular_fallback) const 1981 bool singular_fallback, blas_trans_type transt) const
1980 { 1982 {
1981 FloatComplexMatrix tmp (*this); 1983 FloatComplexMatrix tmp (*this);
1982 return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback); 1984 return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt);
1983 } 1985 }
1984 1986
1985 FloatColumnVector 1987 FloatColumnVector
1986 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const 1988 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b) const
1987 { 1989 {
2004 return solve (typ, b, info, rcon, 0); 2006 return solve (typ, b, info, rcon, 0);
2005 } 2007 }
2006 2008
2007 FloatColumnVector 2009 FloatColumnVector
2008 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info, 2010 FloatMatrix::solve (MatrixType &typ, const FloatColumnVector& b, octave_idx_type& info,
2009 float& rcon, solve_singularity_handler sing_handler) const 2011 float& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const
2010 { 2012 {
2011 FloatMatrix tmp (b); 2013 FloatMatrix tmp (b);
2012 return solve (typ, tmp, info, rcon, sing_handler).column(static_cast<octave_idx_type> (0)); 2014 return solve (typ, tmp, info, rcon, sing_handler, transt).column(static_cast<octave_idx_type> (0));
2013 } 2015 }
2014 2016
2015 FloatComplexColumnVector 2017 FloatComplexColumnVector
2016 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b) const 2018 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b) const
2017 { 2019 {
2036 } 2038 }
2037 2039
2038 FloatComplexColumnVector 2040 FloatComplexColumnVector
2039 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b, 2041 FloatMatrix::solve (MatrixType &typ, const FloatComplexColumnVector& b,
2040 octave_idx_type& info, float& rcon, 2042 octave_idx_type& info, float& rcon,
2041 solve_singularity_handler sing_handler) const 2043 solve_singularity_handler sing_handler, blas_trans_type transt) const
2042 { 2044 {
2043 FloatComplexMatrix tmp (*this); 2045 FloatComplexMatrix tmp (*this);
2044 return tmp.solve(typ, b, info, rcon, sing_handler); 2046 return tmp.solve(typ, b, info, rcon, sing_handler, transt);
2045 } 2047 }
2046 2048
2047 FloatMatrix 2049 FloatMatrix
2048 FloatMatrix::solve (const FloatMatrix& b) const 2050 FloatMatrix::solve (const FloatMatrix& b) const
2049 { 2051 {
2065 return solve (b, info, rcon, 0); 2067 return solve (b, info, rcon, 0);
2066 } 2068 }
2067 2069
2068 FloatMatrix 2070 FloatMatrix
2069 FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info, 2071 FloatMatrix::solve (const FloatMatrix& b, octave_idx_type& info,
2070 float& rcon, solve_singularity_handler sing_handler) const 2072 float& rcon, solve_singularity_handler sing_handler, blas_trans_type transt) const
2071 { 2073 {
2072 MatrixType mattype (*this); 2074 MatrixType mattype (*this);
2073 return solve (mattype, b, info, rcon, sing_handler); 2075 return solve (mattype, b, info, rcon, sing_handler, true, transt);
2074 } 2076 }
2075 2077
2076 FloatComplexMatrix 2078 FloatComplexMatrix
2077 FloatMatrix::solve (const FloatComplexMatrix& b) const 2079 FloatMatrix::solve (const FloatComplexMatrix& b) const
2078 { 2080 {
2094 return tmp.solve (b, info, rcon); 2096 return tmp.solve (b, info, rcon);
2095 } 2097 }
2096 2098
2097 FloatComplexMatrix 2099 FloatComplexMatrix
2098 FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon, 2100 FloatMatrix::solve (const FloatComplexMatrix& b, octave_idx_type& info, float& rcon,
2099 solve_singularity_handler sing_handler) const 2101 solve_singularity_handler sing_handler, blas_trans_type transt) const
2100 { 2102 {
2101 FloatComplexMatrix tmp (*this); 2103 FloatComplexMatrix tmp (*this);
2102 return tmp.solve (b, info, rcon, sing_handler); 2104 return tmp.solve (b, info, rcon, sing_handler, transt);
2103 } 2105 }
2104 2106
2105 FloatColumnVector 2107 FloatColumnVector
2106 FloatMatrix::solve (const FloatColumnVector& b) const 2108 FloatMatrix::solve (const FloatColumnVector& b) const
2107 { 2109 {
2122 return solve (b, info, rcon, 0); 2124 return solve (b, info, rcon, 0);
2123 } 2125 }
2124 2126
2125 FloatColumnVector 2127 FloatColumnVector
2126 FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon, 2128 FloatMatrix::solve (const FloatColumnVector& b, octave_idx_type& info, float& rcon,
2127 solve_singularity_handler sing_handler) const 2129 solve_singularity_handler sing_handler, blas_trans_type transt) const
2128 { 2130 {
2129 MatrixType mattype (*this); 2131 MatrixType mattype (*this);
2130 return solve (mattype, b, info, rcon, sing_handler); 2132 return solve (mattype, b, info, rcon, sing_handler, transt);
2131 } 2133 }
2132 2134
2133 FloatComplexColumnVector 2135 FloatComplexColumnVector
2134 FloatMatrix::solve (const FloatComplexColumnVector& b) const 2136 FloatMatrix::solve (const FloatComplexColumnVector& b) const
2135 { 2137 {
2151 return tmp.solve (b, info, rcon); 2153 return tmp.solve (b, info, rcon);
2152 } 2154 }
2153 2155
2154 FloatComplexColumnVector 2156 FloatComplexColumnVector
2155 FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon, 2157 FloatMatrix::solve (const FloatComplexColumnVector& b, octave_idx_type& info, float& rcon,
2156 solve_singularity_handler sing_handler) const 2158 solve_singularity_handler sing_handler, blas_trans_type transt) const
2157 { 2159 {
2158 FloatComplexMatrix tmp (*this); 2160 FloatComplexMatrix tmp (*this);
2159 return tmp.solve (b, info, rcon, sing_handler); 2161 return tmp.solve (b, info, rcon, sing_handler, transt);
2160 } 2162 }
2161 2163
2162 FloatMatrix 2164 FloatMatrix
2163 FloatMatrix::lssolve (const FloatMatrix& b) const 2165 FloatMatrix::lssolve (const FloatMatrix& b) const
2164 { 2166 {