fCmplxSVD.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1994-2012 John W. Eaton
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include "fCmplxSVD.h"
00028 #include "f77-fcn.h"
00029 #include "lo-error.h"
00030 #include "oct-locbuf.h"
00031 
00032 extern "C"
00033 {
00034   F77_RET_T
00035   F77_FUNC (cgesvd, CGESVD) (F77_CONST_CHAR_ARG_DECL,
00036                              F77_CONST_CHAR_ARG_DECL,
00037                              const octave_idx_type&, const octave_idx_type&,
00038                              FloatComplex*, const octave_idx_type&, float*,
00039                              FloatComplex*, const octave_idx_type&,
00040                              FloatComplex*, const octave_idx_type&,
00041                              FloatComplex*, const octave_idx_type&,
00042                              float*, octave_idx_type&
00043                              F77_CHAR_ARG_LEN_DECL
00044                              F77_CHAR_ARG_LEN_DECL);
00045 
00046   F77_RET_T
00047   F77_FUNC (cgesdd, CGESDD) (F77_CONST_CHAR_ARG_DECL,
00048                              const octave_idx_type&, const octave_idx_type&,
00049                              FloatComplex*, const octave_idx_type&, float*,
00050                              FloatComplex*, const octave_idx_type&,
00051                              FloatComplex*, const octave_idx_type&,
00052                              FloatComplex*, const octave_idx_type&,
00053                              float*, octave_idx_type *, octave_idx_type&
00054                              F77_CHAR_ARG_LEN_DECL);
00055 }
00056 
00057 FloatComplexMatrix
00058 FloatComplexSVD::left_singular_matrix (void) const
00059 {
00060   if (type_computed == SVD::sigma_only)
00061     {
00062       (*current_liboctave_error_handler)
00063         ("FloatComplexSVD: U not computed because type == SVD::sigma_only");
00064       return FloatComplexMatrix ();
00065     }
00066   else
00067     return left_sm;
00068 }
00069 
00070 FloatComplexMatrix
00071 FloatComplexSVD::right_singular_matrix (void) const
00072 {
00073   if (type_computed == SVD::sigma_only)
00074     {
00075       (*current_liboctave_error_handler)
00076         ("FloatComplexSVD: V not computed because type == SVD::sigma_only");
00077       return FloatComplexMatrix ();
00078     }
00079   else
00080     return right_sm;
00081 }
00082 
00083 octave_idx_type
00084 FloatComplexSVD::init (const FloatComplexMatrix& a, SVD::type svd_type,
00085                        SVD::driver svd_driver)
00086 {
00087   octave_idx_type info;
00088 
00089   octave_idx_type m = a.rows ();
00090   octave_idx_type n = a.cols ();
00091 
00092   FloatComplexMatrix atmp = a;
00093   FloatComplex *tmp_data = atmp.fortran_vec ();
00094 
00095   octave_idx_type min_mn = m < n ? m : n;
00096   octave_idx_type max_mn = m > n ? m : n;
00097 
00098   char jobu = 'A';
00099   char jobv = 'A';
00100 
00101   octave_idx_type ncol_u = m;
00102   octave_idx_type nrow_vt = n;
00103   octave_idx_type nrow_s = m;
00104   octave_idx_type ncol_s = n;
00105 
00106   switch (svd_type)
00107     {
00108     case SVD::economy:
00109       jobu = jobv = 'S';
00110       ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
00111       break;
00112 
00113     case SVD::sigma_only:
00114 
00115       // Note:  for this case, both jobu and jobv should be 'N', but
00116       // there seems to be a bug in dgesvd from Lapack V2.0.  To
00117       // demonstrate the bug, set both jobu and jobv to 'N' and find
00118       // the singular values of [eye(3), eye(3)].  The result is
00119       // [-sqrt(2), -sqrt(2), -sqrt(2)].
00120       //
00121       // For Lapack 3.0, this problem seems to be fixed.
00122 
00123       jobu = jobv = 'N';
00124       ncol_u = nrow_vt = 1;
00125       break;
00126 
00127     default:
00128       break;
00129     }
00130 
00131   type_computed = svd_type;
00132 
00133   if (! (jobu == 'N' || jobu == 'O'))
00134     left_sm.resize (m, ncol_u);
00135 
00136   FloatComplex *u = left_sm.fortran_vec ();
00137 
00138   sigma.resize (nrow_s, ncol_s);
00139   float *s_vec = sigma.fortran_vec ();
00140 
00141   if (! (jobv == 'N' || jobv == 'O'))
00142     right_sm.resize (nrow_vt, n);
00143 
00144   FloatComplex *vt = right_sm.fortran_vec ();
00145 
00146   // Query CGESVD for the correct dimension of WORK.
00147 
00148   octave_idx_type lwork = -1;
00149 
00150   Array<FloatComplex> work (dim_vector (1, 1));
00151 
00152   octave_idx_type one = 1;
00153   octave_idx_type m1 = std::max (m, one);
00154   octave_idx_type nrow_vt1 = std::max (nrow_vt, one);
00155 
00156   if (svd_driver == SVD::GESVD)
00157     {
00158       octave_idx_type lrwork = 5*max_mn;
00159       Array<float> rwork (dim_vector (lrwork, 1));
00160 
00161       F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00162                                  F77_CONST_CHAR_ARG2 (&jobv, 1),
00163                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
00164                                  nrow_vt1, work.fortran_vec (), lwork,
00165                                  rwork.fortran_vec (), info
00166                                  F77_CHAR_ARG_LEN (1)
00167                                  F77_CHAR_ARG_LEN (1)));
00168 
00169       lwork = static_cast<octave_idx_type> (work(0).real ());
00170       work.resize (dim_vector (lwork, 1));
00171 
00172       F77_XFCN (cgesvd, CGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00173                                  F77_CONST_CHAR_ARG2 (&jobv, 1),
00174                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
00175                                  nrow_vt1, work.fortran_vec (), lwork,
00176                                  rwork.fortran_vec (), info
00177                                  F77_CHAR_ARG_LEN (1)
00178                                  F77_CHAR_ARG_LEN (1)));
00179     }
00180   else if (svd_driver == SVD::GESDD)
00181     {
00182       assert (jobu == jobv);
00183       char jobz = jobu;
00184 
00185       octave_idx_type lrwork;
00186       if (jobz == 'N')
00187         lrwork = 5*min_mn;
00188       else
00189         lrwork = min_mn * std::max (5*min_mn+7, 2*max_mn+2*min_mn+1);
00190       Array<float> rwork (dim_vector (lrwork, 1));
00191 
00192       OCTAVE_LOCAL_BUFFER (octave_idx_type, iwork, 8*min_mn);
00193 
00194       F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
00195                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
00196                                  nrow_vt1, work.fortran_vec (), lwork,
00197                                  rwork.fortran_vec (), iwork, info
00198                                  F77_CHAR_ARG_LEN (1)));
00199 
00200       lwork = static_cast<octave_idx_type> (work(0).real ());
00201       work.resize (dim_vector (lwork, 1));
00202 
00203       F77_XFCN (cgesdd, CGESDD, (F77_CONST_CHAR_ARG2 (&jobz, 1),
00204                                  m, n, tmp_data, m1, s_vec, u, m1, vt,
00205                                  nrow_vt1, work.fortran_vec (), lwork,
00206                                  rwork.fortran_vec (), iwork, info
00207                                  F77_CHAR_ARG_LEN (1)));
00208     }
00209   else
00210     assert (0); // impossible
00211 
00212   if (! (jobv == 'N' || jobv == 'O'))
00213     right_sm = right_sm.hermitian ();
00214 
00215   return info;
00216 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines