Mercurial > hg > octave-lyh
diff liboctave/fCmplxQRP.cc @ 8597:c86718093c1b
improve & fix QR classes
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Tue, 27 Jan 2009 12:40:06 +0100 |
parents | e3c9102431a9 |
children | 20dfb885f877 |
line wrap: on
line diff
--- a/liboctave/fCmplxQRP.cc +++ b/liboctave/fCmplxQRP.cc @@ -30,6 +30,7 @@ #include "fCmplxQRP.h" #include "f77-fcn.h" #include "lo-error.h" +#include "oct-locbuf.h" extern "C" { @@ -37,11 +38,6 @@ F77_FUNC (cgeqp3, CGEQP3) (const octave_idx_type&, const octave_idx_type&, FloatComplex*, const octave_idx_type&, octave_idx_type*, FloatComplex*, FloatComplex*, const octave_idx_type&, float*, octave_idx_type&); - - F77_RET_T - F77_FUNC (cungqr, CUNGQR) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, - FloatComplex*, const octave_idx_type&, FloatComplex*, - FloatComplex*, const octave_idx_type&, octave_idx_type&); } // It would be best to share some of this code with FloatComplexQR class... @@ -60,44 +56,34 @@ octave_idx_type m = a.rows (); octave_idx_type n = a.cols (); - if (m == 0 || n == 0) - { - (*current_liboctave_error_handler) - ("FloatComplexQR must have non-empty matrix"); - return; - } - octave_idx_type min_mn = m < n ? m : n; - Array<FloatComplex> tau (min_mn); - FloatComplex *ptau = tau.fortran_vec (); + OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn); octave_idx_type info = 0; - FloatComplexMatrix A_fact = a; - if (m > n && qr_type != QR::economy) - A_fact.resize (m, m, 0.0); - - FloatComplex *tmp_data = A_fact.fortran_vec (); - - Array<float> rwork (2*n); - float *prwork = rwork.fortran_vec (); + FloatComplexMatrix afact = a; + if (m > n && qr_type == QR::std) + afact.resize (m, m); MArray<octave_idx_type> jpvt (n, 0); - octave_idx_type *pjpvt = jpvt.fortran_vec (); - FloatComplex rlwork = 0; - // Workspace query... - F77_XFCN (cgeqp3, CGEQP3, (m, n, tmp_data, m, pjpvt, ptau, &rlwork, - -1, prwork, info)); + if (m > 0) + { + OCTAVE_LOCAL_BUFFER (float, rwork, 2*n); - octave_idx_type lwork = rlwork.real (); - Array<FloatComplex> work (lwork); - FloatComplex *pwork = work.fortran_vec (); + // workspace query. + FloatComplex clwork; + F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), + tau, &clwork, -1, rwork, info)); - // Code to enforce a certain permutation could go here... - - F77_XFCN (cgeqp3, CGEQP3, (m, n, tmp_data, m, pjpvt, ptau, pwork, - lwork, prwork, info)); + // allocate buffer and do the job. + octave_idx_type lwork = clwork.real (); lwork = std::max (lwork, 1); + OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork); + F77_XFCN (cgeqp3, CGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (), + tau, work, lwork, rwork, info)); + } + else + for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1; // Form Permutation matrix (if economy is requested, return the // indices only!) @@ -105,25 +91,8 @@ jpvt -= 1; p = PermMatrix (jpvt, true); - octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m; - if (qr_type == QR::economy && m > n) - r.resize (n, n, 0.0); - else - r.resize (m, n, 0.0); - - for (octave_idx_type j = 0; j < n; j++) - { - octave_idx_type limit = j < min_mn-1 ? j : min_mn-1; - for (octave_idx_type i = 0; i <= limit; i++) - r.elem (i, j) = A_fact.elem (i, j); - } - - F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau, - pwork, lwork, info)); - - q = A_fact; - q.resize (m, n2); + form (n, afact, tau, qr_type); } FloatColumnVector