Mercurial > octave-nkf

comparison liboctave/dMatrix.cc @ 7079:6d3e53a2f963

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-30 19:26:32 by jwe]
author jwe
date Tue, 30 Oct 2007 19:26:33 +0000
parents 0bade2dc44a1
children d07cb867891b
comparison
equal deleted inserted replaced
7078:405cf85b435c 7079:6d3e53a2f963
2102 octave_idx_type lwork = -1; 2102 octave_idx_type lwork = -1;
2103 2103
2104 Array<double> work (1); 2104 Array<double> work (1);
2105 2105
2106 // FIXME: Can SMLSIZ be other than 25? 2106 // FIXME: Can SMLSIZ be other than 25?
2107 octave_idx_type liwork = 3 * minmn * 25 + 11 * minmn; 2107 octave_idx_type smlsiz = 25;
2108
2109 // We compute the size of iwork because DGELSD in older versions
2110 // of LAPACK does not return it on a query call.
2111 #if defined (HAVE_LOG2)
2112 double tmp = log2 (minmn) / static_cast<double> (smlsiz+1) + 1;
2113 #else
2114 double tmp = log (minmn) / static_cast<double> (smlsiz+1) / log (2) + 1;
2115 #endif
2116 octave_idx_type nlvl = static_cast<int> (tmp);
2117 if (nlvl < 0)
2118 nlvl = 0;
2119
2120 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2121 if (liwork < 1)
2122 liwork = 1;
2108 Array<octave_idx_type> iwork (liwork); 2123 Array<octave_idx_type> iwork (liwork);
2109 octave_idx_type* piwork = iwork.fortran_vec (); 2124 octave_idx_type* piwork = iwork.fortran_vec ();
2110 2125
2111 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2126 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2112 ps, rcond, rank, work.fortran_vec (), 2127 ps, rcond, rank, work.fortran_vec (),
2135 ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank); 2150 ("dgelsd: rank deficient %dx%d matrix, rank = %d", m, n, rank);
2136 if (s.elem (0) == 0.0) 2151 if (s.elem (0) == 0.0)
2137 rcond = 0.0; 2152 rcond = 0.0;
2138 else 2153 else
2139 rcond = s.elem (minmn - 1) / s.elem (0); 2154 rcond = s.elem (minmn - 1) / s.elem (0);
2155
2156 retval.resize (n, nrhs);
2140 } 2157 }
2141 } 2158 }
2142 } 2159 }
2143 2160
2144 return retval; 2161 return retval;
2248 octave_idx_type lwork = -1; 2265 octave_idx_type lwork = -1;
2249 2266
2250 Array<double> work (1); 2267 Array<double> work (1);
2251 2268
2252 // FIXME: Can SMLSIZ be other than 25? 2269 // FIXME: Can SMLSIZ be other than 25?
2253 octave_idx_type liwork = 3 * minmn * 25 + 11 * minmn; 2270 octave_idx_type smlsiz = 25;
2271
2272 // We compute the size of iwork because DGELSD in older versions
2273 // of LAPACK does not return it on a query call.
2274 #if defined (HAVE_LOG2)
2275 double tmp = log2 (minmn) / static_cast<double> (smlsiz+1) + 1;
2276 #else
2277 double tmp = log (minmn) / static_cast<double> (smlsiz+1) / log (2) + 1;
2278 #endif
2279 octave_idx_type nlvl = static_cast<int> (tmp);
2280 if (nlvl < 0)
2281 nlvl = 0;
2282
2283 octave_idx_type liwork = 3 * minmn * nlvl + 11 * minmn;
2284 if (liwork < 1)
2285 liwork = 1;
2254 Array<octave_idx_type> iwork (liwork); 2286 Array<octave_idx_type> iwork (liwork);
2255 octave_idx_type* piwork = iwork.fortran_vec (); 2287 octave_idx_type* piwork = iwork.fortran_vec ();
2256 2288
2257 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2289 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2258 ps, rcond, rank, work.fortran_vec (), 2290 ps, rcond, rank, work.fortran_vec (),
2282 if (s.elem (0) == 0.0) 2314 if (s.elem (0) == 0.0)
2283 rcond = 0.0; 2315 rcond = 0.0;
2284 else 2316 else
2285 rcond = s.elem (minmn - 1) / s.elem (0); 2317 rcond = s.elem (minmn - 1) / s.elem (0);
2286 } 2318 }
2319
2320 retval.resize (n, nrhs);
2287 } 2321 }
2288 } 2322 }
2289 2323
2290 return retval; 2324 return retval;
2291 } 2325 }