fCmplxCHOL.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1994-2012 John W. Eaton
00004 Copyright (C) 2008-2009 Jaroslav Hajek
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 <vector>
00029 
00030 #include "fMatrix.h"
00031 #include "fRowVector.h"
00032 #include "fCmplxCHOL.h"
00033 #include "f77-fcn.h"
00034 #include "lo-error.h"
00035 #include "oct-locbuf.h"
00036 #include "oct-norm.h"
00037 #ifndef HAVE_QRUPDATE
00038 #include "dbleQR.h"
00039 #endif
00040 
00041 extern "C"
00042 {
00043   F77_RET_T
00044   F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
00045                              const octave_idx_type&, FloatComplex*,
00046                              const octave_idx_type&, octave_idx_type&
00047                              F77_CHAR_ARG_LEN_DECL);
00048   F77_RET_T
00049   F77_FUNC (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL,
00050                              const octave_idx_type&, FloatComplex*,
00051                              const octave_idx_type&, octave_idx_type&
00052                              F77_CHAR_ARG_LEN_DECL);
00053 
00054   F77_RET_T
00055   F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL,
00056                              const octave_idx_type&, FloatComplex*,
00057                              const octave_idx_type&, const float&,
00058                              float&, FloatComplex*, float*, octave_idx_type&
00059                              F77_CHAR_ARG_LEN_DECL);
00060 #ifdef HAVE_QRUPDATE
00061 
00062   F77_RET_T
00063   F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*,
00064                              const octave_idx_type&, FloatComplex*, float*);
00065 
00066   F77_RET_T
00067   F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*,
00068                              const octave_idx_type&, FloatComplex*,
00069                              float*, octave_idx_type&);
00070 
00071   F77_RET_T
00072   F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*,
00073                              const octave_idx_type&, const octave_idx_type&,
00074                              FloatComplex*, float*, octave_idx_type&);
00075 
00076   F77_RET_T
00077   F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*,
00078                              const octave_idx_type&, const octave_idx_type&,
00079                              float*);
00080 
00081   F77_RET_T
00082   F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*,
00083                              const octave_idx_type&, const octave_idx_type&,
00084                              const octave_idx_type&, FloatComplex*, float*);
00085 #endif
00086 }
00087 
00088 octave_idx_type
00089 FloatComplexCHOL::init (const FloatComplexMatrix& a, bool calc_cond)
00090 {
00091   octave_idx_type a_nr = a.rows ();
00092   octave_idx_type a_nc = a.cols ();
00093 
00094   if (a_nr != a_nc)
00095     {
00096       (*current_liboctave_error_handler)
00097         ("FloatComplexCHOL requires square matrix");
00098       return -1;
00099     }
00100 
00101   octave_idx_type n = a_nc;
00102   octave_idx_type info;
00103 
00104   chol_mat.clear (n, n);
00105   for (octave_idx_type j = 0; j < n; j++)
00106     {
00107       for (octave_idx_type i = 0; i <= j; i++)
00108         chol_mat.xelem (i, j) = a(i, j);
00109       for (octave_idx_type i = j+1; i < n; i++)
00110         chol_mat.xelem (i, j) = 0.0f;
00111     }
00112   FloatComplex *h = chol_mat.fortran_vec ();
00113 
00114   // Calculate the norm of the matrix, for later use.
00115   float anorm = 0;
00116   if (calc_cond)
00117     anorm = xnorm (a, 1);
00118 
00119   F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
00120                              F77_CHAR_ARG_LEN (1)));
00121 
00122   xrcond = 0.0;
00123   if (info > 0)
00124     chol_mat.resize (info - 1, info - 1);
00125   else if (calc_cond)
00126     {
00127       octave_idx_type cpocon_info = 0;
00128 
00129       // Now calculate the condition number for non-singular matrix.
00130       Array<FloatComplex> z (dim_vector (2*n, 1));
00131       FloatComplex *pz = z.fortran_vec ();
00132       Array<float> rz (dim_vector (n, 1));
00133       float *prz = rz.fortran_vec ();
00134       F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
00135                                  n, anorm, xrcond, pz, prz, cpocon_info
00136                                  F77_CHAR_ARG_LEN (1)));
00137 
00138       if (cpocon_info != 0)
00139         info = -1;
00140     }
00141 
00142   return info;
00143 }
00144 
00145 static FloatComplexMatrix
00146 chol2inv_internal (const FloatComplexMatrix& r)
00147 {
00148   FloatComplexMatrix retval;
00149 
00150   octave_idx_type r_nr = r.rows ();
00151   octave_idx_type r_nc = r.cols ();
00152 
00153   if (r_nr == r_nc)
00154     {
00155       octave_idx_type n = r_nc;
00156       octave_idx_type info;
00157 
00158       FloatComplexMatrix tmp = r;
00159 
00160       F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
00161                                  tmp.fortran_vec (), n, info
00162                                  F77_CHAR_ARG_LEN (1)));
00163 
00164       // If someone thinks of a more graceful way of doing this (or
00165       // faster for that matter :-)), please let me know!
00166 
00167       if (n > 1)
00168         for (octave_idx_type j = 0; j < r_nc; j++)
00169           for (octave_idx_type i = j+1; i < r_nr; i++)
00170             tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
00171 
00172       retval = tmp;
00173     }
00174   else
00175     (*current_liboctave_error_handler) ("chol2inv requires square matrix");
00176 
00177   return retval;
00178 }
00179 
00180 // Compute the inverse of a matrix using the Cholesky factorization.
00181 FloatComplexMatrix
00182 FloatComplexCHOL::inverse (void) const
00183 {
00184   return chol2inv_internal (chol_mat);
00185 }
00186 
00187 void
00188 FloatComplexCHOL::set (const FloatComplexMatrix& R)
00189 {
00190   if (R.is_square ())
00191     chol_mat = R;
00192   else
00193     (*current_liboctave_error_handler) ("CHOL requires square matrix");
00194 }
00195 
00196 #ifdef HAVE_QRUPDATE
00197 
00198 void
00199 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
00200 {
00201   octave_idx_type n = chol_mat.rows ();
00202 
00203   if (u.length () == n)
00204     {
00205       FloatComplexColumnVector utmp = u;
00206 
00207       OCTAVE_LOCAL_BUFFER (float, rw, n);
00208 
00209       F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
00210                                  utmp.fortran_vec (), rw));
00211     }
00212   else
00213     (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
00214 }
00215 
00216 octave_idx_type
00217 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
00218 {
00219   octave_idx_type info = -1;
00220 
00221   octave_idx_type n = chol_mat.rows ();
00222 
00223   if (u.length () == n)
00224     {
00225       FloatComplexColumnVector utmp = u;
00226 
00227       OCTAVE_LOCAL_BUFFER (float, rw, n);
00228 
00229       F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
00230                                  utmp.fortran_vec (), rw, info));
00231     }
00232   else
00233     (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
00234 
00235   return info;
00236 }
00237 
00238 octave_idx_type
00239 FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j)
00240 {
00241   octave_idx_type info = -1;
00242 
00243   octave_idx_type n = chol_mat.rows ();
00244 
00245   if (u.length () != n + 1)
00246     (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
00247   else if (j < 0 || j > n)
00248     (*current_liboctave_error_handler) ("cholinsert: index out of range");
00249   else
00250     {
00251       FloatComplexColumnVector utmp = u;
00252 
00253       OCTAVE_LOCAL_BUFFER (float, rw, n);
00254 
00255       chol_mat.resize (n+1, n+1);
00256 
00257       F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
00258                                  j + 1, utmp.fortran_vec (), rw, info));
00259     }
00260 
00261   return info;
00262 }
00263 
00264 void
00265 FloatComplexCHOL::delete_sym (octave_idx_type j)
00266 {
00267   octave_idx_type n = chol_mat.rows ();
00268 
00269   if (j < 0 || j > n-1)
00270     (*current_liboctave_error_handler) ("choldelete: index out of range");
00271   else
00272     {
00273       OCTAVE_LOCAL_BUFFER (float, rw, n);
00274 
00275       F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
00276                                  j + 1, rw));
00277 
00278       chol_mat.resize (n-1, n-1);
00279     }
00280 }
00281 
00282 void
00283 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
00284 {
00285   octave_idx_type n = chol_mat.rows ();
00286 
00287   if (i < 0 || i > n-1 || j < 0 || j > n-1)
00288     (*current_liboctave_error_handler) ("cholshift: index out of range");
00289   else
00290     {
00291       OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
00292       OCTAVE_LOCAL_BUFFER (float, rw, n);
00293 
00294       F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
00295                                  i + 1, j + 1, w, rw));
00296     }
00297 }
00298 
00299 #else
00300 
00301 void
00302 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
00303 {
00304   warn_qrupdate_once ();
00305 
00306   octave_idx_type n = chol_mat.rows ();
00307 
00308   if (u.length () == n)
00309     {
00310       init (chol_mat.hermitian () * chol_mat
00311             + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
00312     }
00313   else
00314     (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
00315 }
00316 
00317 static bool
00318 singular (const FloatComplexMatrix& a)
00319 {
00320   for (octave_idx_type i = 0; i < a.rows (); i++)
00321     if (a(i,i) == 0.0f) return true;
00322   return false;
00323 }
00324 
00325 octave_idx_type
00326 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
00327 {
00328   warn_qrupdate_once ();
00329 
00330   octave_idx_type info = -1;
00331 
00332   octave_idx_type n = chol_mat.rows ();
00333 
00334   if (u.length () == n)
00335     {
00336       if (singular (chol_mat))
00337         info = 2;
00338       else
00339         {
00340           info = init (chol_mat.hermitian () * chol_mat
00341                        - FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
00342           if (info) info = 1;
00343         }
00344     }
00345   else
00346     (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
00347 
00348   return info;
00349 }
00350 
00351 octave_idx_type
00352 FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j)
00353 {
00354   warn_qrupdate_once ();
00355 
00356   octave_idx_type info = -1;
00357 
00358   octave_idx_type n = chol_mat.rows ();
00359 
00360   if (u.length () != n + 1)
00361     (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
00362   else if (j < 0 || j > n)
00363     (*current_liboctave_error_handler) ("cholinsert: index out of range");
00364   else
00365     {
00366       if (singular (chol_mat))
00367         info = 2;
00368       else if (u(j).imag () != 0.0)
00369         info = 3;
00370       else
00371         {
00372           FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
00373           FloatComplexMatrix a1 (n+1, n+1);
00374           for (octave_idx_type k = 0; k < n+1; k++)
00375             for (octave_idx_type l = 0; l < n+1; l++)
00376               {
00377                 if (l == j)
00378                   a1(k, l) = u(k);
00379                 else if (k == j)
00380                   a1(k, l) = std::conj (u(l));
00381                 else
00382                   a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
00383               }
00384           info = init (a1, false);
00385           if (info) info = 1;
00386         }
00387     }
00388 
00389   return info;
00390 }
00391 
00392 void
00393 FloatComplexCHOL::delete_sym (octave_idx_type j)
00394 {
00395   warn_qrupdate_once ();
00396 
00397   octave_idx_type n = chol_mat.rows ();
00398 
00399   if (j < 0 || j > n-1)
00400     (*current_liboctave_error_handler) ("choldelete: index out of range");
00401   else
00402     {
00403       FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
00404       a.delete_elements (1, idx_vector (j));
00405       a.delete_elements (0, idx_vector (j));
00406       init (a, false);
00407     }
00408 }
00409 
00410 void
00411 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
00412 {
00413   warn_qrupdate_once ();
00414 
00415   octave_idx_type n = chol_mat.rows ();
00416 
00417   if (i < 0 || i > n-1 || j < 0 || j > n-1)
00418     (*current_liboctave_error_handler) ("cholshift: index out of range");
00419   else
00420     {
00421       FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
00422       Array<octave_idx_type> p (dim_vector (n, 1));
00423       for (octave_idx_type k = 0; k < n; k++) p(k) = k;
00424       if (i < j)
00425         {
00426           for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
00427           p(j) = i;
00428         }
00429       else if (j < i)
00430         {
00431           p(j) = i;
00432           for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
00433         }
00434 
00435       init (a.index (idx_vector (p), idx_vector (p)), false);
00436     }
00437 }
00438 
00439 #endif
00440 
00441 FloatComplexMatrix
00442 chol2inv (const FloatComplexMatrix& r)
00443 {
00444   return chol2inv_internal (r);
00445 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines