Mercurial > hg > octave-lyh

comparison liboctave/CMatrix.cc @ 7476:e9f10b4c05cf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix workspace size calculation for xGELSD
author Jason Riedy
date Tue, 12 Feb 2008 21:04:27 -0500
parents a7a987b229b7
children 8b22207ef9ca
comparison
equal deleted inserted replaced
7475:aa5208636bea 7476:e9f10b4c05cf
61 61
62 // Fortran functions we call. 62 // Fortran functions we call.
63 63
64 extern "C" 64 extern "C"
65 { 65 {
66 octave_idx_type
67 F77_FUNC (ilaenv, ILAENV) (const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
68 F77_CONST_CHAR_ARG_DECL,
69 const octave_idx_type&, const octave_idx_type&,
70 const octave_idx_type&, const octave_idx_type&
71 F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
72
66 F77_RET_T 73 F77_RET_T
67 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL, 74 F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
68 const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type&, 75 const octave_idx_type&, Complex*, const octave_idx_type&, octave_idx_type&,
69 octave_idx_type&, double*, octave_idx_type& 76 octave_idx_type&, double*, octave_idx_type&
70 F77_CHAR_ARG_LEN_DECL); 77 F77_CHAR_ARG_LEN_DECL);
2490 // Ask ZGELSD what the dimension of WORK should be. 2497 // Ask ZGELSD what the dimension of WORK should be.
2491 octave_idx_type lwork = -1; 2498 octave_idx_type lwork = -1;
2492 2499
2493 Array<Complex> work (1); 2500 Array<Complex> work (1);
2494 2501
2495 // FIXME: Can SMLSIZ be other than 25? 2502 const octave_idx_type smlsiz
2496 octave_idx_type smlsiz = 25; 2503 = F77_FUNC (ilaenv, ILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6),
2504 F77_CONST_CHAR_ARG2 (" ", 1),
2505 0, 0, 0, 0
2506 F77_CHAR_ARG_LEN (6)
2507 F77_CHAR_ARG_LEN (1));
2497 2508
2498 // We compute the size of rwork and iwork because ZGELSD in 2509 // We compute the size of rwork and iwork because ZGELSD in
2499 // older versions of LAPACK does not return them on a query 2510 // older versions of LAPACK does not return them on a query
2500 // call. 2511 // call.
2501 double dminmn = static_cast<double> (minmn); 2512 double dminmn = static_cast<double> (minmn);
2523 octave_idx_type* piwork = iwork.fortran_vec (); 2534 octave_idx_type* piwork = iwork.fortran_vec ();
2524 2535
2525 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2536 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2526 ps, rcond, rank, work.fortran_vec (), 2537 ps, rcond, rank, work.fortran_vec (),
2527 lwork, prwork, piwork, info)); 2538 lwork, prwork, piwork, info));
2539
2540 // The workspace query is broken in at least LAPACK 3.0.0
2541 // through 3.1.1 when n > m. The obtuse formula below
2542 // should provide sufficient workspace for DGELSD to operate
2543 // efficiently.
2544 if (n > m)
2545 {
2546 octave_idx_type addend = m;
2547
2548 if (2*m-4 > addend)
2549 addend = 2*m-4;
2550
2551 if (nrhs > addend)
2552 addend = nrhs;
2553
2554 if (n-3*m > addend)
2555 addend = n-3*m;
2556
2557 const octave_idx_type lworkaround = 4*m + m*m + addend;
2558
2559 if (std::real (work(0)) < lworkaround)
2560 work(0) = lworkaround;
2561 }
2528 2562
2529 if (f77_exception_encountered) 2563 if (f77_exception_encountered)
2530 (*current_liboctave_error_handler) 2564 (*current_liboctave_error_handler)
2531 ("unrecoverable error in zgelsd"); 2565 ("unrecoverable error in zgelsd");
2532 else 2566 else