Mercurial > hg > octave-lyh

comparison liboctave/dMatrix.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
57 57
58 // Fortran functions we call. 58 // Fortran functions we call.
59 59
60 extern "C" 60 extern "C"
61 { 61 {
62 octave_idx_type
63 F77_FUNC (ilaenv, ILAENV) (const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
64 F77_CONST_CHAR_ARG_DECL,
65 const octave_idx_type&, const octave_idx_type&,
66 const octave_idx_type&, const octave_idx_type&
67 F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);
68
62 F77_RET_T 69 F77_RET_T
63 F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL, 70 F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
64 const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type&, 71 const octave_idx_type&, double*, const octave_idx_type&, octave_idx_type&,
65 octave_idx_type&, double*, octave_idx_type& 72 octave_idx_type&, double*, octave_idx_type&
66 F77_CHAR_ARG_LEN_DECL); 73 F77_CHAR_ARG_LEN_DECL);
2101 // Ask DGELSD what the dimension of WORK should be. 2108 // Ask DGELSD what the dimension of WORK should be.
2102 octave_idx_type lwork = -1; 2109 octave_idx_type lwork = -1;
2103 2110
2104 Array<double> work (1); 2111 Array<double> work (1);
2105 2112
2106 // FIXME: Can SMLSIZ be other than 25? 2113 const octave_idx_type smlsiz
2107 octave_idx_type smlsiz = 25; 2114 = F77_FUNC (ilaenv, ILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2115 F77_CONST_CHAR_ARG2 (" ", 1),
2116 0, 0, 0, 0
2117 F77_CHAR_ARG_LEN (6)
2118 F77_CHAR_ARG_LEN (1));
2108 2119
2109 // We compute the size of iwork because DGELSD in older versions 2120 // We compute the size of iwork because DGELSD in older versions
2110 // of LAPACK does not return it on a query call. 2121 // of LAPACK does not return it on a query call.
2111 double dminmn = static_cast<double> (minmn); 2122 double dminmn = static_cast<double> (minmn);
2112 double dsmlsizp1 = static_cast<double> (smlsiz+1); 2123 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2126 octave_idx_type* piwork = iwork.fortran_vec (); 2137 octave_idx_type* piwork = iwork.fortran_vec ();
2127 2138
2128 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn, 2139 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2129 ps, rcond, rank, work.fortran_vec (), 2140 ps, rcond, rank, work.fortran_vec (),
2130 lwork, piwork, info)); 2141 lwork, piwork, info));
2142
2143 // The workspace query is broken in at least LAPACK 3.0.0
2144 // through 3.1.1 when n > m. The obtuse formula below
2145 // should provide sufficient workspace for DGELSD to operate
2146 // efficiently.
2147 if (n > m)
2148 {
2149 const octave_idx_type wlalsd
2150 = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
2151
2152 octave_idx_type addend = m;
2153
2154 if (2*m-4 > addend)
2155 addend = 2*m-4;
2156
2157 if (nrhs > addend)
2158 addend = nrhs;
2159
2160 if (n-3*m > addend)
2161 addend = n-3*m;
2162
2163 if (wlalsd > addend)
2164 addend = wlalsd;
2165
2166 const octave_idx_type lworkaround = 4*m + m*m + addend;
2167
2168 if (work(0) < lworkaround)
2169 work(0) = lworkaround;
2170 }
2131 2171
2132 if (f77_exception_encountered) 2172 if (f77_exception_encountered)
2133 (*current_liboctave_error_handler) 2173 (*current_liboctave_error_handler)
2134 ("unrecoverable error in dgelsd"); 2174 ("unrecoverable error in dgelsd");
2135 else 2175 else