CmplxQRP.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1994-2012 John W. Eaton
00004 Copyright (C) 2009 VZLU Prague
00005 
00006 This file is part of Octave.
00007 
00008 Octave is free software; you can redistribute it and/or modify it
00009 under the terms of the GNU General Public License as published by the
00010 Free Software Foundation; either version 3 of the License, or (at your
00011 option) any later version.
00012 
00013 Octave is distributed in the hope that it will be useful, but WITHOUT
00014 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00016 for more details.
00017 
00018 You should have received a copy of the GNU General Public License
00019 along with Octave; see the file COPYING.  If not, see
00020 <http://www.gnu.org/licenses/>.
00021 
00022 */
00023 
00024 #ifdef HAVE_CONFIG_H
00025 #include <config.h>
00026 #endif
00027 
00028 #include <cassert>
00029 
00030 #include "CmplxQRP.h"
00031 #include "f77-fcn.h"
00032 #include "lo-error.h"
00033 #include "oct-locbuf.h"
00034 
00035 extern "C"
00036 {
00037   F77_RET_T
00038   F77_FUNC (zgeqp3, ZGEQP3) (const octave_idx_type&, const octave_idx_type&,
00039                              Complex*, const octave_idx_type&,
00040                              octave_idx_type*, Complex*, Complex*,
00041                              const octave_idx_type&, double*,
00042                              octave_idx_type&);
00043 }
00044 
00045 // It would be best to share some of this code with ComplexQR class...
00046 
00047 ComplexQRP::ComplexQRP (const ComplexMatrix& a, qr_type_t qr_type)
00048   : ComplexQR (), p ()
00049 {
00050   init (a, qr_type);
00051 }
00052 
00053 void
00054 ComplexQRP::init (const ComplexMatrix& a, qr_type_t qr_type)
00055 {
00056   assert (qr_type != qr_type_raw);
00057 
00058   octave_idx_type m = a.rows ();
00059   octave_idx_type n = a.cols ();
00060 
00061   octave_idx_type min_mn = m < n ? m : n;
00062   OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
00063 
00064   octave_idx_type info = 0;
00065 
00066   ComplexMatrix afact = a;
00067   if (m > n && qr_type == qr_type_std)
00068     afact.resize (m, m);
00069 
00070   MArray<octave_idx_type> jpvt (dim_vector (n, 1), 0);
00071 
00072   if (m > 0)
00073     {
00074       OCTAVE_LOCAL_BUFFER (double, rwork, 2*n);
00075 
00076       // workspace query.
00077       Complex clwork;
00078       F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
00079                                  tau, &clwork, -1, rwork, info));
00080 
00081       // allocate buffer and do the job.
00082       octave_idx_type lwork = clwork.real ();
00083       lwork = std::max (lwork, static_cast<octave_idx_type> (1));
00084       OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
00085       F77_XFCN (zgeqp3, ZGEQP3, (m, n, afact.fortran_vec (), m, jpvt.fortran_vec (),
00086                                  tau, work, lwork, rwork, info));
00087     }
00088   else
00089     for (octave_idx_type i = 0; i < n; i++) jpvt(i) = i+1;
00090 
00091   // Form Permutation matrix (if economy is requested, return the
00092   // indices only!)
00093 
00094   jpvt -= static_cast<octave_idx_type> (1);
00095   p = PermMatrix (jpvt, true);
00096 
00097 
00098   form (n, afact, tau, qr_type);
00099 }
00100 
00101 RowVector
00102 ComplexQRP::Pvec (void) const
00103 {
00104   Array<double> pa (p.pvec ());
00105   RowVector pv (MArray<double> (pa) + 1.0);
00106   return pv;
00107 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines