GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fCmplxCHOL.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2015 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <vector>
29 
30 #include "fMatrix.h"
31 #include "fRowVector.h"
32 #include "fCmplxCHOL.h"
33 #include "f77-fcn.h"
34 #include "lo-error.h"
35 #include "oct-locbuf.h"
36 #include "oct-norm.h"
37 #ifndef HAVE_QRUPDATE
38 #include "dbleQR.h"
39 #endif
40 
41 extern "C"
42 {
43  F77_RET_T
44  F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL,
46  const octave_idx_type&, octave_idx_type&
48  F77_RET_T
49  F77_FUNC (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL,
50  const octave_idx_type&, FloatComplex*,
51  const octave_idx_type&, octave_idx_type&
53 
54  F77_RET_T
55  F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL,
56  const octave_idx_type&, FloatComplex*,
57  const octave_idx_type&, const float&,
58  float&, FloatComplex*, float*, octave_idx_type&
60 #ifdef HAVE_QRUPDATE
61 
62  F77_RET_T
63  F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*,
64  const octave_idx_type&, FloatComplex*, float*);
65 
66  F77_RET_T
67  F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*,
68  const octave_idx_type&, FloatComplex*,
69  float*, octave_idx_type&);
70 
71  F77_RET_T
72  F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*,
73  const octave_idx_type&, const octave_idx_type&,
74  FloatComplex*, float*, octave_idx_type&);
75 
76  F77_RET_T
77  F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*,
78  const octave_idx_type&, const octave_idx_type&,
79  float*);
80 
81  F77_RET_T
82  F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*,
83  const octave_idx_type&, const octave_idx_type&,
84  const octave_idx_type&, FloatComplex*, float*);
85 #endif
86 }
87 
89 FloatComplexCHOL::init (const FloatComplexMatrix& a, bool calc_cond)
90 {
91  octave_idx_type a_nr = a.rows ();
92  octave_idx_type a_nc = a.cols ();
93 
94  if (a_nr != a_nc)
95  {
96  (*current_liboctave_error_handler)
97  ("FloatComplexCHOL requires square matrix");
98  return -1;
99  }
100 
101  octave_idx_type n = a_nc;
102  octave_idx_type info;
103 
104  chol_mat.clear (n, n);
105  for (octave_idx_type j = 0; j < n; j++)
106  {
107  for (octave_idx_type i = 0; i <= j; i++)
108  chol_mat.xelem (i, j) = a(i, j);
109  for (octave_idx_type i = j+1; i < n; i++)
110  chol_mat.xelem (i, j) = 0.0f;
111  }
113 
114  // Calculate the norm of the matrix, for later use.
115  float anorm = 0;
116  if (calc_cond)
117  anorm = xnorm (a, 1);
118 
119  F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
120  F77_CHAR_ARG_LEN (1)));
121 
122  xrcond = 0.0;
123  if (info > 0)
124  chol_mat.resize (info - 1, info - 1);
125  else if (calc_cond)
126  {
127  octave_idx_type cpocon_info = 0;
128 
129  // Now calculate the condition number for non-singular matrix.
130  Array<FloatComplex> z (dim_vector (2*n, 1));
131  FloatComplex *pz = z.fortran_vec ();
132  Array<float> rz (dim_vector (n, 1));
133  float *prz = rz.fortran_vec ();
134  F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
135  n, anorm, xrcond, pz, prz, cpocon_info
136  F77_CHAR_ARG_LEN (1)));
137 
138  if (cpocon_info != 0)
139  info = -1;
140  }
141 
142  return info;
143 }
144 
145 static FloatComplexMatrix
147 {
148  FloatComplexMatrix retval;
149 
150  octave_idx_type r_nr = r.rows ();
151  octave_idx_type r_nc = r.cols ();
152 
153  if (r_nr == r_nc)
154  {
155  octave_idx_type n = r_nc;
156  octave_idx_type info;
157 
158  FloatComplexMatrix tmp = r;
159 
160  F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
161  tmp.fortran_vec (), n, info
162  F77_CHAR_ARG_LEN (1)));
163 
164  // If someone thinks of a more graceful way of doing this (or
165  // faster for that matter :-)), please let me know!
166 
167  if (n > 1)
168  for (octave_idx_type j = 0; j < r_nc; j++)
169  for (octave_idx_type i = j+1; i < r_nr; i++)
170  tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
171 
172  retval = tmp;
173  }
174  else
175  (*current_liboctave_error_handler) ("chol2inv requires square matrix");
176 
177  return retval;
178 }
179 
180 // Compute the inverse of a matrix using the Cholesky factorization.
183 {
184  return chol2inv_internal (chol_mat);
185 }
186 
187 void
189 {
190  if (R.is_square ())
191  chol_mat = R;
192  else
193  (*current_liboctave_error_handler) ("CHOL requires square matrix");
194 }
195 
196 #ifdef HAVE_QRUPDATE
197 
198 void
200 {
202 
203  if (u.length () == n)
204  {
205  FloatComplexColumnVector utmp = u;
206 
207  OCTAVE_LOCAL_BUFFER (float, rw, n);
208 
209  F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
210  utmp.fortran_vec (), rw));
211  }
212  else
213  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
214 }
215 
218 {
219  octave_idx_type info = -1;
220 
222 
223  if (u.length () == n)
224  {
225  FloatComplexColumnVector utmp = u;
226 
227  OCTAVE_LOCAL_BUFFER (float, rw, n);
228 
229  F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
230  utmp.fortran_vec (), rw, info));
231  }
232  else
233  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
234 
235  return info;
236 }
237 
240  octave_idx_type j)
241 {
242  octave_idx_type info = -1;
243 
245 
246  if (u.length () != n + 1)
247  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
248  else if (j < 0 || j > n)
249  (*current_liboctave_error_handler) ("cholinsert: index out of range");
250  else
251  {
252  FloatComplexColumnVector utmp = u;
253 
254  OCTAVE_LOCAL_BUFFER (float, rw, n);
255 
256  chol_mat.resize (n+1, n+1);
257 
258  F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
259  j + 1, utmp.fortran_vec (), rw, info));
260  }
261 
262  return info;
263 }
264 
265 void
267 {
269 
270  if (j < 0 || j > n-1)
271  (*current_liboctave_error_handler) ("choldelete: index out of range");
272  else
273  {
274  OCTAVE_LOCAL_BUFFER (float, rw, n);
275 
276  F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
277  j + 1, rw));
278 
279  chol_mat.resize (n-1, n-1);
280  }
281 }
282 
283 void
285 {
287 
288  if (i < 0 || i > n-1 || j < 0 || j > n-1)
289  (*current_liboctave_error_handler) ("cholshift: index out of range");
290  else
291  {
293  OCTAVE_LOCAL_BUFFER (float, rw, n);
294 
295  F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
296  i + 1, j + 1, w, rw));
297  }
298 }
299 
300 #else
301 
302 void
304 {
305  warn_qrupdate_once ();
306 
308 
309  if (u.length () == n)
310  {
313  false);
314  }
315  else
316  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
317 }
318 
319 static bool
320 singular (const FloatComplexMatrix& a)
321 {
322  for (octave_idx_type i = 0; i < a.rows (); i++)
323  if (a(i,i) == 0.0f) return true;
324  return false;
325 }
326 
329 {
330  warn_qrupdate_once ();
331 
332  octave_idx_type info = -1;
333 
335 
336  if (u.length () == n)
337  {
338  if (singular (chol_mat))
339  info = 2;
340  else
341  {
342  info = init (chol_mat.hermitian () * chol_mat
343  - FloatComplexMatrix (u)
344  * FloatComplexMatrix (u).hermitian (),
345  false);
346  if (info) info = 1;
347  }
348  }
349  else
350  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
351 
352  return info;
353 }
354 
357  octave_idx_type j)
358 {
359  warn_qrupdate_once ();
360 
361  octave_idx_type info = -1;
362 
364 
365  if (u.length () != n + 1)
366  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
367  else if (j < 0 || j > n)
368  (*current_liboctave_error_handler) ("cholinsert: index out of range");
369  else
370  {
371  if (singular (chol_mat))
372  info = 2;
373  else if (u(j).imag () != 0.0)
374  info = 3;
375  else
376  {
378  FloatComplexMatrix a1 (n+1, n+1);
379  for (octave_idx_type k = 0; k < n+1; k++)
380  for (octave_idx_type l = 0; l < n+1; l++)
381  {
382  if (l == j)
383  a1(k, l) = u(k);
384  else if (k == j)
385  a1(k, l) = std::conj (u(l));
386  else
387  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
388  }
389  info = init (a1, false);
390  if (info) info = 1;
391  }
392  }
393 
394  return info;
395 }
396 
397 void
399 {
400  warn_qrupdate_once ();
401 
403 
404  if (j < 0 || j > n-1)
405  (*current_liboctave_error_handler) ("choldelete: index out of range");
406  else
407  {
409  a.delete_elements (1, idx_vector (j));
410  a.delete_elements (0, idx_vector (j));
411  init (a, false);
412  }
413 }
414 
415 void
417 {
418  warn_qrupdate_once ();
419 
421 
422  if (i < 0 || i > n-1 || j < 0 || j > n-1)
423  (*current_liboctave_error_handler) ("cholshift: index out of range");
424  else
425  {
428  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
429  if (i < j)
430  {
431  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
432  p(j) = i;
433  }
434  else if (j < i)
435  {
436  p(j) = i;
437  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
438  }
439 
440  init (a.index (idx_vector (p), idx_vector (p)), false);
441  }
442 }
443 
444 #endif
445 
448 {
449  return chol2inv_internal (r);
450 }
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:175
#define F77_CHAR_ARG_LEN(l)
Definition: f77-fcn.h:253
F77_RET_T F77_FUNC(cpotrf, CPOTRF)(F77_CONST_CHAR_ARG_DECL
octave_idx_type insert_sym(const FloatComplexColumnVector &u, octave_idx_type j)
Definition: fCmplxCHOL.cc:239
octave_idx_type init(const FloatComplexMatrix &a, bool calc_cond)
Definition: fCmplxCHOL.cc:89
FloatComplexMatrix hermitian(void) const
Definition: fCMatrix.h:154
octave_idx_type downdate(const FloatComplexColumnVector &u)
Definition: fCmplxCHOL.cc:217
void delete_elements(const idx_vector &i)
Deleting elements.
Definition: Array.cc:1393
static FloatComplexMatrix chol2inv_internal(const FloatComplexMatrix &r)
Definition: fCmplxCHOL.cc:146
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
FloatComplexMatrix inverse(void) const
Definition: fCmplxCHOL.cc:182
octave_idx_type rows(void) const
Definition: Array.h:313
#define F77_CONST_CHAR_ARG2(x, l)
Definition: f77-fcn.h:251
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:244
FloatComplexMatrix chol2inv(const FloatComplexMatrix &r)
Definition: fCmplxCHOL.cc:447
std::complex< double > w(std::complex< double > z, double relerr=0)
void shift_sym(octave_idx_type i, octave_idx_type j)
Definition: fCmplxCHOL.cc:284
bool is_square(void) const
Definition: Array.h:470
OCTAVE_API double xnorm(const ColumnVector &x, double p)
Definition: oct-norm.cc:536
#define F77_RET_T
Definition: f77-fcn.h:264
FloatComplexMatrix chol_mat
Definition: fCmplxCHOL.h:91
T & xelem(octave_idx_type n)
Definition: Array.h:353
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
void clear(void)
Definition: Array.cc:84
#define F77_CONST_CHAR_ARG_DECL
Definition: f77-fcn.h:255
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
void set(const FloatComplexMatrix &R)
Definition: fCmplxCHOL.cc:188
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
const T * fortran_vec(void) const
Definition: Array.h:481
octave_idx_type cols(void) const
Definition: Array.h:321
void update(const FloatComplexColumnVector &u)
Definition: fCmplxCHOL.cc:199
void delete_sym(octave_idx_type j)
Definition: fCmplxCHOL.cc:266
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716
F77_RET_T const octave_idx_type FloatComplex const octave_idx_type octave_idx_type & F77_CHAR_ARG_LEN_DECL
Definition: fCmplxCHOL.cc:45