GNU Octave  3.8.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
CmplxCHOL.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2013 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 "dMatrix.h"
31 #include "dRowVector.h"
32 #include "CmplxCHOL.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 (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL,
46  const octave_idx_type&, octave_idx_type&
48  F77_RET_T
49  F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL,
50  const octave_idx_type&, Complex*,
51  const octave_idx_type&, octave_idx_type&
53 
54  F77_RET_T
55  F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL,
56  const octave_idx_type&, Complex*,
57  const octave_idx_type&, const double&,
58  double&, Complex*, double*, octave_idx_type&
60 #ifdef HAVE_QRUPDATE
61 
62  F77_RET_T
63  F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*,
64  const octave_idx_type&, Complex*, double*);
65 
66  F77_RET_T
67  F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*,
68  const octave_idx_type&, Complex*, double*,
69  octave_idx_type&);
70 
71  F77_RET_T
72  F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*,
73  const octave_idx_type&, const octave_idx_type&,
74  Complex*, double*, octave_idx_type&);
75 
76  F77_RET_T
77  F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*,
78  const octave_idx_type&, const octave_idx_type&,
79  double*);
80 
81  F77_RET_T
82  F77_FUNC (zchshx, ZCHSHX) (const octave_idx_type&, Complex*,
83  const octave_idx_type&, const octave_idx_type&,
84  const octave_idx_type&, Complex*, double*);
85 #endif
86 }
87 
89 ComplexCHOL::init (const ComplexMatrix& 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  ("ComplexCHOL 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.0;
111  }
112  Complex *h = chol_mat.fortran_vec ();
113 
114  // Calculate the norm of the matrix, for later use.
115  double anorm = 0;
116  if (calc_cond)
117  anorm = xnorm (a, 1);
118 
119  F77_XFCN (zpotrf, ZPOTRF, (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 zpocon_info = 0;
128 
129  // Now calculate the condition number for non-singular matrix.
130  Array<Complex> z (dim_vector (2*n, 1));
131  Complex *pz = z.fortran_vec ();
132  Array<double> rz (dim_vector (n, 1));
133  double *prz = rz.fortran_vec ();
134  F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
135  n, anorm, xrcond, pz, prz, zpocon_info
136  F77_CHAR_ARG_LEN (1)));
137 
138  if (zpocon_info != 0)
139  info = -1;
140  }
141 
142  return info;
143 }
144 
145 static ComplexMatrix
147 {
148  ComplexMatrix 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  ComplexMatrix tmp = r;
159 
160  F77_XFCN (zpotri, ZPOTRI, (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  ComplexColumnVector utmp = u;
206 
207  OCTAVE_LOCAL_BUFFER (double, rw, n);
208 
209  F77_XFCN (zch1up, ZCH1UP, (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  ComplexColumnVector utmp = u;
226 
227  OCTAVE_LOCAL_BUFFER (double, rw, n);
228 
229  F77_XFCN (zch1dn, ZCH1DN, (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 {
241  octave_idx_type info = -1;
242 
244 
245  if (u.length () != n + 1)
246  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
247  else if (j < 0 || j > n)
248  (*current_liboctave_error_handler) ("cholinsert: index out of range");
249  else
250  {
251  ComplexColumnVector utmp = u;
252 
253  OCTAVE_LOCAL_BUFFER (double, rw, n);
254 
255  chol_mat.resize (n+1, n+1);
256 
257  F77_XFCN (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
258  j + 1, utmp.fortran_vec (), rw, info));
259  }
260 
261  return info;
262 }
263 
264 void
266 {
268 
269  if (j < 0 || j > n-1)
270  (*current_liboctave_error_handler) ("choldelete: index out of range");
271  else
272  {
273  OCTAVE_LOCAL_BUFFER (double, rw, n);
274 
275  F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
276  j + 1, rw));
277 
278  chol_mat.resize (n-1, n-1);
279  }
280 }
281 
282 void
284 {
286 
287  if (i < 0 || i > n-1 || j < 0 || j > n-1)
288  (*current_liboctave_error_handler) ("cholshift: index out of range");
289  else
290  {
292  OCTAVE_LOCAL_BUFFER (double, rw, n);
293 
294  F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
295  i + 1, j + 1, w, rw));
296  }
297 }
298 
299 #else
300 
301 void
303 {
304  warn_qrupdate_once ();
305 
307 
308  if (u.length () == n)
309  {
311  + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
312  }
313  else
314  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
315 }
316 
317 static bool
318 singular (const ComplexMatrix& a)
319 {
320  for (octave_idx_type i = 0; i < a.rows (); i++)
321  if (a(i,i) == 0.0) return true;
322  return false;
323 }
324 
327 {
328  warn_qrupdate_once ();
329 
330  octave_idx_type info = -1;
331 
333 
334  if (u.length () == n)
335  {
336  if (singular (chol_mat))
337  info = 2;
338  else
339  {
340  info = init (chol_mat.hermitian () * chol_mat
341  - ComplexMatrix (u) * ComplexMatrix (u).hermitian (),
342  false);
343  if (info) info = 1;
344  }
345  }
346  else
347  (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
348 
349  return info;
350 }
351 
354 {
355  warn_qrupdate_once ();
356 
357  octave_idx_type info = -1;
358 
360 
361  if (u.length () != n + 1)
362  (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
363  else if (j < 0 || j > n)
364  (*current_liboctave_error_handler) ("cholinsert: index out of range");
365  else
366  {
367  if (singular (chol_mat))
368  info = 2;
369  else if (u(j).imag () != 0.0)
370  info = 3;
371  else
372  {
374  ComplexMatrix a1 (n+1, n+1);
375  for (octave_idx_type k = 0; k < n+1; k++)
376  for (octave_idx_type l = 0; l < n+1; l++)
377  {
378  if (l == j)
379  a1(k, l) = u(k);
380  else if (k == j)
381  a1(k, l) = std::conj (u(l));
382  else
383  a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
384  }
385  info = init (a1, false);
386  if (info) info = 1;
387  }
388  }
389 
390  return info;
391 }
392 
393 void
395 {
396  warn_qrupdate_once ();
397 
399 
400  if (j < 0 || j > n-1)
401  (*current_liboctave_error_handler) ("choldelete: index out of range");
402  else
403  {
405  a.delete_elements (1, idx_vector (j));
406  a.delete_elements (0, idx_vector (j));
407  init (a, false);
408  }
409 }
410 
411 void
413 {
414  warn_qrupdate_once ();
415 
417 
418  if (i < 0 || i > n-1 || j < 0 || j > n-1)
419  (*current_liboctave_error_handler) ("cholshift: index out of range");
420  else
421  {
424  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
425  if (i < j)
426  {
427  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
428  p(j) = i;
429  }
430  else if (j < i)
431  {
432  p(j) = i;
433  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
434  }
435 
436  init (a.index (idx_vector (p), idx_vector (p)), false);
437  }
438 }
439 
440 #endif
441 
444 {
445  return chol2inv_internal (r);
446 }