Mercurial > octave-nkf

comparison libinterp/corefcn/lu.cc @ 18425:2a45b6b87bee stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
correct numerical errors in sparse LU factorization (bug #41116). * lu.cc: modified to apply pivots as warranted to L and U. * sparse-base-lu.cc: compute correct matrix size for single-output case.
author Michael C. Grant <mcg@cvxr.com>
date Mon, 03 Feb 2014 10:54:27 -0500
parents 175b392e91fe
children 807e21ebddd5
comparison
equal deleted inserted replaced
18418:1ad77b3e6bef 18425:2a45b6b87bee
204 if (arg_is_empty < 0) 204 if (arg_is_empty < 0)
205 return retval; 205 return retval;
206 else if (arg_is_empty > 0) 206 else if (arg_is_empty > 0)
207 return octave_value_list (5, SparseMatrix ()); 207 return octave_value_list (5, SparseMatrix ());
208 208
209 ColumnVector Qinit;
210 if (nargout < 4)
211 {
212 Qinit.resize (nc);
213 for (octave_idx_type i = 0; i < nc; i++)
214 Qinit (i) = i;
215 }
216
217 if (arg.is_real_type ()) 209 if (arg.is_real_type ())
218 { 210 {
211
219 SparseMatrix m = arg.sparse_matrix_value (); 212 SparseMatrix m = arg.sparse_matrix_value ();
220 213
221 switch (nargout) 214 if ( nargout < 4 )
222 { 215 {
223 case 0: 216
224 case 1: 217 ColumnVector Qinit;
225 case 2: 218 Qinit.resize (nc);
226 { 219 for (octave_idx_type i = 0; i < nc; i++)
227 SparseLU fact (m, Qinit, thres, false, true); 220 Qinit (i) = i;
228 221 SparseLU fact (m, Qinit, thres, false, true);
229 if (nargout < 2) 222
223 if ( nargout < 2 )
230 retval(0) = fact.Y (); 224 retval(0) = fact.Y ();
231 else 225 else
232 { 226 {
233 PermMatrix P = fact.Pr_mat (); 227
234 SparseMatrix L = P.transpose () * fact.L (); 228 retval(1)
235 retval(1) = octave_value (fact.U (), 229 = octave_value (
236 MatrixType (MatrixType::Upper)); 230 fact.U () * fact.Pc_mat ().transpose (),
237 231 MatrixType (MatrixType::Permuted_Upper,
238 retval(0) 232 nc, fact.col_perm ()));
239 = octave_value (L, MatrixType (MatrixType::Permuted_Lower, 233
240 nr, fact.row_perm ())); 234 PermMatrix P = fact.Pr_mat ();
241 } 235 SparseMatrix L = fact.L ();
242 } 236 if ( nargout < 3 )
243 break; 237 retval(0)
244 238 = octave_value ( P.transpose () * L,
245 case 3: 239 MatrixType (MatrixType::Permuted_Lower,
246 { 240 nr, fact.row_perm ()));
247 SparseLU fact (m, Qinit, thres, false, true); 241 else
248 242 {
249 if (vecout) 243 retval(0) = L;
250 retval(2) = fact.Pr_vec (); 244 if ( vecout )
251 else 245 retval(2) = fact.Pr_vec();
252 retval(2) = fact.Pr_mat (); 246 else
253 247 retval(2) = P;
254 retval(1) = octave_value (fact.U (), 248 }
255 MatrixType (MatrixType::Upper)); 249
256 retval(0) = octave_value (fact.L (), 250 }
257 MatrixType (MatrixType::Lower)); 251
258 } 252 }
259 break; 253 else
260 254 {
261 case 4: 255
262 default:
263 {
264 SparseLU fact (m, thres, scale); 256 SparseLU fact (m, thres, scale);
265 257
266 if (scale) 258 if (scale)
267 retval(4) = fact.R (); 259 retval(4) = fact.R ();
268 260
278 } 270 }
279 retval(1) = octave_value (fact.U (), 271 retval(1) = octave_value (fact.U (),
280 MatrixType (MatrixType::Upper)); 272 MatrixType (MatrixType::Upper));
281 retval(0) = octave_value (fact.L (), 273 retval(0) = octave_value (fact.L (),
282 MatrixType (MatrixType::Lower)); 274 MatrixType (MatrixType::Lower));
283 } 275 }
284 break; 276
285 }
286 } 277 }
287 else if (arg.is_complex_type ()) 278 else if (arg.is_complex_type ())
288 { 279 {
289 SparseComplexMatrix m = arg.sparse_complex_matrix_value (); 280 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
290 281
291 switch (nargout) 282 if ( nargout < 4 )
292 { 283 {
293 case 0: 284
294 case 1: 285 ColumnVector Qinit;
295 case 2: 286 Qinit.resize (nc);
296 { 287 for (octave_idx_type i = 0; i < nc; i++)
297 SparseComplexLU fact (m, Qinit, thres, false, true); 288 Qinit (i) = i;
298 289 SparseComplexLU fact (m, Qinit, thres, false, true);
299 if (nargout < 2) 290
291 if ( nargout < 2 )
292
300 retval(0) = fact.Y (); 293 retval(0) = fact.Y ();
301 else 294
302 { 295 else
303 PermMatrix P = fact.Pr_mat (); 296 {
304 SparseComplexMatrix L = P.transpose () * fact.L (); 297
305 retval(1) = octave_value (fact.U (), 298 retval(1)
306 MatrixType (MatrixType::Upper)); 299 = octave_value (
307 300 fact.U () * fact.Pc_mat ().transpose (),
308 retval(0) 301 MatrixType (MatrixType::Permuted_Upper,
309 = octave_value (L, MatrixType (MatrixType::Permuted_Lower, 302 nc, fact.col_perm ()));
310 nr, fact.row_perm ())); 303
311 } 304 PermMatrix P = fact.Pr_mat ();
312 } 305 SparseComplexMatrix L = fact.L ();
313 break; 306 if ( nargout < 3 )
314 307 retval(0)
315 case 3: 308 = octave_value ( P.transpose () * L,
316 { 309 MatrixType (MatrixType::Permuted_Lower,
317 SparseComplexLU fact (m, Qinit, thres, false, true); 310 nr, fact.row_perm ()));
318 311 else
319 if (vecout) 312 {
320 retval(2) = fact.Pr_vec (); 313 retval(0) = L;
321 else 314 if ( vecout )
322 retval(2) = fact.Pr_mat (); 315 retval(2) = fact.Pr_vec();
323 316 else
324 retval(1) = octave_value (fact.U (), 317 retval(2) = P;
325 MatrixType (MatrixType::Upper)); 318 }
326 retval(0) = octave_value (fact.L (), 319
327 MatrixType (MatrixType::Lower)); 320 }
328 } 321
329 break; 322 }
330 323 else
331 case 4: 324 {
332 default: 325
333 {
334 SparseComplexLU fact (m, thres, scale); 326 SparseComplexLU fact (m, thres, scale);
335 327
336 if (scale) 328 if (scale)
337 retval(4) = fact.R (); 329 retval(4) = fact.R ();
338 330
348 } 340 }
349 retval(1) = octave_value (fact.U (), 341 retval(1) = octave_value (fact.U (),
350 MatrixType (MatrixType::Upper)); 342 MatrixType (MatrixType::Upper));
351 retval(0) = octave_value (fact.L (), 343 retval(0) = octave_value (fact.L (),
352 MatrixType (MatrixType::Lower)); 344 MatrixType (MatrixType::Lower));
353 } 345 }
354 break; 346
355 }
356 } 347 }
357 else 348 else
358 gripe_wrong_type_arg ("lu", arg); 349 gripe_wrong_type_arg ("lu", arg);
359 } 350 }
360 else 351 else
580 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single"))); 571 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
581 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single"))); 572 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
582 573
583 %!error lu () 574 %!error lu ()
584 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2) 575 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
576
577 %!test
578 %! Bi = [1 2 3 4 5 2 3 6 7 8 4 5 7 8 9];
579 %! Bj = [1 3 4 5 6 7 8 9 11 12 13 14 15 16 17];
580 %! Bv = [1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1];
581 %! B = sparse (Bi,Bj,Bv);
582 %! [L,U] = lu (B);
583 %! assert (nnz (B - L*U) == 0);
584 %! [L,U,P] = lu(B);
585 %! assert (nnz (B - P'*L*U) == 0);
586 %! [L,U,P,Q] = lu (B);
587 %! assert (nnz (B - P'*L*U*Q') == 0);
588
585 */ 589 */
586 590
587 static 591 static
588 bool check_lu_dims (const octave_value& l, const octave_value& u, 592 bool check_lu_dims (const octave_value& l, const octave_value& u,
589 const octave_value& p) 593 const octave_value& p)