Mercurial > hg > octave-nkf

comparison liboctave/fCmplxQR.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 141b3fb5cef7
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
40 40
41 extern "C" 41 extern "C"
42 { 42 {
43 F77_RET_T 43 F77_RET_T
44 F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, 44 F77_FUNC (cgeqrf, CGEQRF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*,
45 const octave_idx_type&, FloatComplex*, FloatComplex*, 45 const octave_idx_type&, FloatComplex*, FloatComplex*,
46 const octave_idx_type&, octave_idx_type&); 46 const octave_idx_type&, octave_idx_type&);
47 47
48 F77_RET_T 48 F77_RET_T
49 F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 49 F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
50 FloatComplex*, const octave_idx_type&, FloatComplex*, 50 FloatComplex*, const octave_idx_type&, FloatComplex*,
51 FloatComplex*, const octave_idx_type&, octave_idx_type&); 51 FloatComplex*, const octave_idx_type&, octave_idx_type&);
52 52
53 #ifdef HAVE_QRUPDATE 53 #ifdef HAVE_QRUPDATE
54 54
55 F77_RET_T 55 F77_RET_T
56 F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 56 F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
129 octave_idx_type info; 129 octave_idx_type info;
130 130
131 if (qr_type == qr_type_raw) 131 if (qr_type == qr_type_raw)
132 { 132 {
133 for (octave_idx_type j = 0; j < min_mn; j++) 133 for (octave_idx_type j = 0; j < min_mn; j++)
134 { 134 {
135 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1; 135 octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
136 for (octave_idx_type i = limit + 1; i < m; i++) 136 for (octave_idx_type i = limit + 1; i < m; i++)
137 afact.elem (i, j) *= tau[j]; 137 afact.elem (i, j) *= tau[j];
138 } 138 }
139 139
140 r = afact; 140 r = afact;
141 } 141 }
142 else 142 else
143 { 143 {
180 F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau, 180 F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
181 &clwork, -1, info)); 181 &clwork, -1, info));
182 182
183 // allocate buffer and do the job. 183 // allocate buffer and do the job.
184 octave_idx_type lwork = clwork.real (); 184 octave_idx_type lwork = clwork.real ();
185 lwork = std::max (lwork, static_cast<octave_idx_type> (1)); 185 lwork = std::max (lwork, static_cast<octave_idx_type> (1));
186 OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork); 186 OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
187 F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau, 187 F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
188 work, lwork, info)); 188 work, lwork, info));
189 } 189 }
190 } 190 }
298 } 298 }
299 299
300 OCTAVE_LOCAL_BUFFER (float, rw, kmax); 300 OCTAVE_LOCAL_BUFFER (float, rw, kmax);
301 for (volatile octave_idx_type i = 0; i < js.length (); i++) 301 for (volatile octave_idx_type i = 0; i < js.length (); i++)
302 { 302 {
303 octave_idx_type ii = i; 303 octave_idx_type ii = i;
304 F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii), 304 F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii),
305 q.fortran_vec (), q.rows (), 305 q.fortran_vec (), q.rows (),
306 r.fortran_vec (), r.rows (), js(ii) + 1, 306 r.fortran_vec (), r.rows (), js(ii) + 1,
307 u.column (jsi(i)).data (), rw)); 307 u.column (jsi(i)).data (), rw));
308 } 308 }
320 (*current_liboctave_error_handler) ("qrdelete: index out of range"); 320 (*current_liboctave_error_handler) ("qrdelete: index out of range");
321 else 321 else
322 { 322 {
323 OCTAVE_LOCAL_BUFFER (float, rw, k); 323 OCTAVE_LOCAL_BUFFER (float, rw, k);
324 F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), q.rows (), 324 F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
325 r.fortran_vec (), r.rows (), j + 1, rw)); 325 r.fortran_vec (), r.rows (), j + 1, rw));
326 326
327 if (k < m) 327 if (k < m)
328 { 328 {
329 q.resize (m, k-1); 329 q.resize (m, k-1);
330 r.resize (k-1, n-1); 330 r.resize (k-1, n-1);
357 else if (nj > 0) 357 else if (nj > 0)
358 { 358 {
359 OCTAVE_LOCAL_BUFFER (float, rw, k); 359 OCTAVE_LOCAL_BUFFER (float, rw, k);
360 for (volatile octave_idx_type i = 0; i < js.length (); i++) 360 for (volatile octave_idx_type i = 0; i < js.length (); i++)
361 { 361 {
362 octave_idx_type ii = i; 362 octave_idx_type ii = i;
363 F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii, 363 F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii,
364 q.fortran_vec (), q.rows (), 364 q.fortran_vec (), q.rows (),
365 r.fortran_vec (), r.rows (), js(ii) + 1, rw)); 365 r.fortran_vec (), r.rows (), js(ii) + 1, rw));
366 } 366 }
367 if (k < m) 367 if (k < m)
393 q.resize (m + 1, m + 1); 393 q.resize (m + 1, m + 1);
394 r.resize (m + 1, n); 394 r.resize (m + 1, n);
395 FloatComplexRowVector utmp = u; 395 FloatComplexRowVector utmp = u;
396 OCTAVE_LOCAL_BUFFER (float, rw, k); 396 OCTAVE_LOCAL_BUFFER (float, rw, k);
397 F77_XFCN (cqrinr, CQRINR, (m, n, q.fortran_vec (), q.rows (), 397 F77_XFCN (cqrinr, CQRINR, (m, n, q.fortran_vec (), q.rows (),
398 r.fortran_vec (), r.rows (), 398 r.fortran_vec (), r.rows (),
399 j + 1, utmp.fortran_vec (), rw)); 399 j + 1, utmp.fortran_vec (), rw));
400 400
401 } 401 }
402 } 402 }
403 403
414 else 414 else
415 { 415 {
416 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m); 416 OCTAVE_LOCAL_BUFFER (FloatComplex, w, m);
417 OCTAVE_LOCAL_BUFFER (float, rw, m); 417 OCTAVE_LOCAL_BUFFER (float, rw, m);
418 F77_XFCN (cqrder, CQRDER, (m, n, q.fortran_vec (), q.rows (), 418 F77_XFCN (cqrder, CQRDER, (m, n, q.fortran_vec (), q.rows (),
419 r.fortran_vec (), r.rows (), j + 1, 419 r.fortran_vec (), r.rows (), j + 1,
420 w, rw)); 420 w, rw));
421 421
422 q.resize (m - 1, m - 1); 422 q.resize (m - 1, m - 1);
423 r.resize (m - 1, n); 423 r.resize (m - 1, n);
424 } 424 }