Mercurial > hg > octave-nkf

comparison liboctave/dMatrix.cc @ 9662:0d3b248f4ab6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
further improve mixed real-complex division
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 23 Sep 2009 11:10:52 +0200
parents afcf852256d2
children 1dba57e9d08d
comparison
equal deleted inserted replaced
9661:afcf852256d2 9662:0d3b248f4ab6
1909 double rcon; 1909 double rcon;
1910 return solve (typ, b, info, rcon, 0); 1910 return solve (typ, b, info, rcon, 0);
1911 } 1911 }
1912 1912
1913 Matrix 1913 Matrix
1914 Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info) const
1915 {
1916 double rcon;
1917 return solve (typ, b, info, rcon, 0);
1918 }
1919
1920 Matrix
1914 Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info, 1921 Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
1915 double& rcon) const 1922 double& rcon) const
1916 { 1923 {
1917 return solve (typ, b, info, rcon, 0); 1924 return solve (typ, b, info, rcon, 0);
1918 } 1925 }
1954 } 1961 }
1955 1962
1956 ComplexMatrix 1963 ComplexMatrix
1957 Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const 1964 Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const
1958 { 1965 {
1959 ComplexMatrix tmp (*this); 1966 octave_idx_type info;
1960 return tmp.solve (typ, b); 1967 double rcon;
1968 return solve (typ, b, info, rcon, 0);
1961 } 1969 }
1962 1970
1963 ComplexMatrix 1971 ComplexMatrix
1964 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, 1972 Matrix::solve (MatrixType &typ, const ComplexMatrix& b,
1965 octave_idx_type& info) const 1973 octave_idx_type& info) const
1966 { 1974 {
1967 ComplexMatrix tmp (*this); 1975 double rcon;
1968 return tmp.solve (typ, b, info); 1976 return solve (typ, b, info, rcon, 0);
1969 } 1977 }
1970 1978
1971 ComplexMatrix 1979 ComplexMatrix
1972 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, 1980 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
1973 double& rcon) const 1981 double& rcon) const
1974 { 1982 {
1975 ComplexMatrix tmp (*this); 1983 return solve (typ, b, info, rcon, 0);
1976 return tmp.solve (typ, b, info, rcon); 1984 }
1985
1986 static Matrix
1987 stack_complex_matrix (const ComplexMatrix& cm)
1988 {
1989 octave_idx_type m = cm.rows (), n = cm.cols (), nel = m*n;
1990 Matrix retval (m, 2*n);
1991 const Complex *cmd = cm.data ();
1992 double *rd = retval.fortran_vec ();
1993 for (octave_idx_type i = 0; i < nel; i++)
1994 {
1995 rd[i] = std::real (cmd[i]);
1996 rd[nel+i] = std::imag (cmd[i]);
1997 }
1998 return retval;
1999 }
2000
2001 static ComplexMatrix
2002 unstack_complex_matrix (const Matrix& sm)
2003 {
2004 octave_idx_type m = sm.rows (), n = sm.cols () / 2, nel = m*n;
2005 ComplexMatrix retval (m, n);
2006 const double *smd = sm.data ();
2007 Complex *rd = retval.fortran_vec ();
2008 for (octave_idx_type i = 0; i < nel; i++)
2009 rd[i] = Complex (smd[i], smd[nel+i]);
2010 return retval;
1977 } 2011 }
1978 2012
1979 ComplexMatrix 2013 ComplexMatrix
1980 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info, 2014 Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
1981 double& rcon, solve_singularity_handler sing_handler, 2015 double& rcon, solve_singularity_handler sing_handler,
1982 bool singular_fallback, blas_trans_type transt) const 2016 bool singular_fallback, blas_trans_type transt) const
1983 { 2017 {
1984 ComplexMatrix tmp (*this); 2018 Matrix tmp = stack_complex_matrix (b);
1985 return tmp.solve (typ, b, info, rcon, sing_handler, singular_fallback, transt); 2019 tmp = solve (typ, tmp, info, rcon, sing_handler, singular_fallback, transt);
2020 return unstack_complex_matrix (tmp);
1986 } 2021 }
1987 2022
1988 ColumnVector 2023 ColumnVector
1989 Matrix::solve (MatrixType &typ, const ColumnVector& b) const 2024 Matrix::solve (MatrixType &typ, const ColumnVector& b) const
1990 { 2025 {