Mercurial > octave-nkf

comparison liboctave/array/CMatrix.cc @ 20633:ffc6cdcd02c5 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix segfault when complex double matrix calls ZGETRF (bug #45577). * CMatrix.cc (finverse, determinant, rcond, fsolve): Calculate norm of matrix and if it is NaN, skip calling ZGETRF in LAPACK and set info to non-zero value to signal an error.
author Rik <rik@octave.org>
date Sat, 10 Oct 2015 16:46:00 -0700
parents 17d647821d61
children 1c5a86b7f838
comparison
equal deleted inserted replaced
20632:7890893a0e69 20633:ffc6cdcd02c5
1129 1129
1130 info = 0; 1130 info = 0;
1131 1131
1132 // Calculate the norm of the matrix, for later use. 1132 // Calculate the norm of the matrix, for later use.
1133 double anorm; 1133 double anorm;
1134 if (calc_cond) 1134 //if (calc_cond) // Must always calculate anorm for bug #45577
1135 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)) 1135 anorm = retval.abs ().sum ().row (static_cast<octave_idx_type>(0)).max ();
1136 .max (); 1136
1137 1137 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1138 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info)); 1138 if (xisnan (anorm))
1139 info = -1;
1140 else
1141 F77_XFCN (zgetrf, ZGETRF, (nc, nc, tmp_data, nr, pipvt, info));
1139 1142
1140 // Throw-away extra info LAPACK gives so as to not change output. 1143 // Throw-away extra info LAPACK gives so as to not change output.
1141 rcon = 0.0; 1144 rcon = 0.0;
1142 if (info != 0) 1145 if (info != 0)
1143 info = -1; 1146 info = -1;
1696 1699
1697 info = 0; 1700 info = 0;
1698 1701
1699 // Calculate the norm of the matrix, for later use. 1702 // Calculate the norm of the matrix, for later use.
1700 double anorm = 0; 1703 double anorm = 0;
1701 if (calc_cond) anorm = xnorm (*this, 1); 1704 //if (calc_cond) // Must always calculate anorm for bug #45577
1702 1705 anorm = xnorm (*this, 1);
1703 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1706
1707 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1708 if (xisnan (anorm))
1709 info = -1;
1710 else
1711 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1704 1712
1705 // Throw-away extra info LAPACK gives so as to not change output. 1713 // Throw-away extra info LAPACK gives so as to not change output.
1706 rcon = 0.0; 1714 rcon = 0.0;
1707 if (info != 0) 1715 if (info != 0)
1708 { 1716 {
1889 Array<Complex> z (dim_vector (2 * nc, 1)); 1897 Array<Complex> z (dim_vector (2 * nc, 1));
1890 Complex *pz = z.fortran_vec (); 1898 Complex *pz = z.fortran_vec ();
1891 Array<double> rz (dim_vector (2 * nc, 1)); 1899 Array<double> rz (dim_vector (2 * nc, 1));
1892 double *prz = rz.fortran_vec (); 1900 double *prz = rz.fortran_vec ();
1893 1901
1894 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 1902 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1903 if (xisnan (anorm))
1904 info = -1;
1905 else
1906 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1895 1907
1896 if (info != 0) 1908 if (info != 0)
1897 { 1909 {
1898 rcon = 0.0; 1910 rcon = 0.0;
1899 mattype.mark_as_rectangular (); 1911 mattype.mark_as_rectangular ();
2223 // Calculate the norm of the matrix, for later use. 2235 // Calculate the norm of the matrix, for later use.
2224 if (anorm < 0.) 2236 if (anorm < 0.)
2225 anorm = atmp.abs ().sum ().row (static_cast<octave_idx_type>(0)) 2237 anorm = atmp.abs ().sum ().row (static_cast<octave_idx_type>(0))
2226 .max (); 2238 .max ();
2227 2239
2228 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info)); 2240 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
2241 if (xisnan (anorm))
2242 info = -2;
2243 else
2244 F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
2229 2245
2230 // Throw-away extra info LAPACK gives so as to not change output. 2246 // Throw-away extra info LAPACK gives so as to not change output.
2231 rcon = 0.0; 2247 rcon = 0.0;
2232 if (info != 0) 2248 if (info != 0)
2233 { 2249 {