Mercurial > hg > octave-lyh

annotate liboctave/fCmplxCHOL.cc @ 11188:4cb1522e4d0f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use function handle as input to cellfun, rather than quoted function name or anonymous function wrapper.
author Rik <octave@nomad.inbox5.com>
date Wed, 03 Nov 2010 17:20:56 -0700
parents 12884915a8e4
children 8a5e980da6aa
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
1 /*
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
2
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
4 John W. Eaton
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
5 Copyright (C) 2008, 2009 Jaroslav Hajek
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
6
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
7 This file is part of Octave.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
8
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
9 Octave is free software; you can redistribute it and/or modify it
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
10 under the terms of the GNU General Public License as published by the
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
11 Free Software Foundation; either version 3 of the License, or (at your
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
12 option) any later version.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
13
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
14 Octave is distributed in the hope that it will be useful, but WITHOUT
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
17 for more details.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
18
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
19 You should have received a copy of the GNU General Public License
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
20 along with Octave; see the file COPYING. If not, see
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
21 <http://www.gnu.org/licenses/>.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
22
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
23 */
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
24
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
25 #ifdef HAVE_CONFIG_H
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
26 #include <config.h>
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
27 #endif
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
28
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
29 #include <vector>
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
30
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
31 #include "fMatrix.h"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
32 #include "fRowVector.h"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
33 #include "fCmplxCHOL.h"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
34 #include "f77-fcn.h"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
35 #include "lo-error.h"
8377
25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents: 7789
diff changeset
36 #include "oct-locbuf.h"
9862
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
37 #include "oct-norm.h"
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
38 #ifndef HAVE_QRUPDATE
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
39 #include "dbleQR.h"
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
40 #endif
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
41
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 extern "C"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
44 F77_RET_T
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
45 F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
46 FloatComplex*, const octave_idx_type&, octave_idx_type&
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
47 F77_CHAR_ARG_LEN_DECL);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
48 F77_RET_T
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 F77_FUNC (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
50 FloatComplex*, const octave_idx_type&, octave_idx_type&
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
51 F77_CHAR_ARG_LEN_DECL);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
52
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 F77_RET_T
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
54 F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
55 FloatComplex*, const octave_idx_type&, const float&,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
56 float&, FloatComplex*, float*,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
57 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
58 #ifdef HAVE_QRUPDATE
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
59
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
60 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
61 F77_FUNC (cch1up, CCH1UP) (const octave_idx_type&, FloatComplex*, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
62 FloatComplex*, float*);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
63
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
64 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
65 F77_FUNC (cch1dn, CCH1DN) (const octave_idx_type&, FloatComplex*, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
66 FloatComplex*, float*, octave_idx_type&);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
67
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
68 F77_RET_T
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
69 F77_FUNC (cchinx, CCHINX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
70 const octave_idx_type&, FloatComplex*, float*,
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
71 octave_idx_type&);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
72
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
73 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
74 F77_FUNC (cchdex, CCHDEX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
75 const octave_idx_type&, float*);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
76
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
77 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
78 F77_FUNC (cchshx, CCHSHX) (const octave_idx_type&, FloatComplex*, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
79 const octave_idx_type&, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
80 FloatComplex*, float*);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
81 #endif
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
82 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
83
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 octave_idx_type
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
85 FloatComplexCHOL::init (const FloatComplexMatrix& a, bool calc_cond)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 octave_idx_type a_nr = a.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 octave_idx_type a_nc = a.cols ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
89
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 if (a_nr != a_nc)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
91 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
92 (*current_liboctave_error_handler)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
93 ("FloatComplexCHOL requires square matrix");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 return -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
95 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
96
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 octave_idx_type n = a_nc;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
98 octave_idx_type info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
99
9862
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
100 chol_mat.clear (n, n);
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
101 for (octave_idx_type j = 0; j < n; j++)
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
102 {
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
103 for (octave_idx_type i = 0; i <= j; i++)
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
104 chol_mat.xelem (i, j) = a(i, j);
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
105 for (octave_idx_type i = j+1; i < n; i++)
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
106 chol_mat.xelem (i, j) = 0.0f;
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
107 }
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 FloatComplex *h = chol_mat.fortran_vec ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
109
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 // Calculate the norm of the matrix, for later use.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 float anorm = 0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 if (calc_cond)
9862
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
113 anorm = xnorm (a, 1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
114
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
116 F77_CHAR_ARG_LEN (1)));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
117
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 xrcond = 0.0;
9862
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
119 if (info > 0)
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
120 chol_mat.resize (info - 1, info - 1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 else if (calc_cond)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
122 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
123 octave_idx_type cpocon_info = 0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
124
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 // Now calculate the condition number for non-singular matrix.
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
126 Array<FloatComplex> z (2*n, 1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 FloatComplex *pz = z.fortran_vec ();
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
128 Array<float> rz (n, 1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
129 float *prz = rz.fortran_vec ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
130 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
131 n, anorm, xrcond, pz, prz, cpocon_info
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
132 F77_CHAR_ARG_LEN (1)));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
133
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 if (cpocon_info != 0)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
135 info = -1;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
136 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
137
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
139 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
140
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 static FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 chol2inv_internal (const FloatComplexMatrix& r)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
144 FloatComplexMatrix retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
145
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 octave_idx_type r_nr = r.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
147 octave_idx_type r_nc = r.cols ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
148
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 if (r_nr == r_nc)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 octave_idx_type n = r_nc;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
152 octave_idx_type info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
153
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
154 FloatComplexMatrix tmp = r;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
155
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
157 tmp.fortran_vec (), n, info
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
158 F77_CHAR_ARG_LEN (1)));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
159
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
160 // If someone thinks of a more graceful way of doing this (or
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
161 // faster for that matter :-)), please let me know!
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
162
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 if (n > 1)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
164 for (octave_idx_type j = 0; j < r_nc; j++)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
165 for (octave_idx_type i = j+1; i < r_nr; i++)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
166 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
167
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 retval = tmp;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
171 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
172
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
174 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
175
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
176 // Compute the inverse of a matrix using the Cholesky factorization.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 FloatComplexCHOL::inverse (void) const
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
180 return chol2inv_internal (chol_mat);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
181 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
182
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 FloatComplexCHOL::set (const FloatComplexMatrix& R)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 if (R.is_square ())
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 chol_mat = R;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 (*current_liboctave_error_handler) ("CHOL requires square matrix");
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
190 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
191
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
192 #ifdef HAVE_QRUPDATE
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
193
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
194 void
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
195 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
197 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
198
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 if (u.length () == n)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
200 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
201 FloatComplexColumnVector utmp = u;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
202
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
203 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
204
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
205 F77_XFCN (cch1up, CCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
206 utmp.fortran_vec (), rw));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
208 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
209 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
210 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
211
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
212 octave_idx_type
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
213 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
215 octave_idx_type info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
216
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
217 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
218
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 if (u.length () == n)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
220 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
221 FloatComplexColumnVector utmp = u;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
222
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
223 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
224
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
225 F77_XFCN (cch1dn, CCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
226 utmp.fortran_vec (), rw, info));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
228 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
229 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
230
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
231 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
232 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
233
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
234 octave_idx_type
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
235 FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
236 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
237 octave_idx_type info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
238
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
239 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
240
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
241 if (u.length () != n + 1)
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
242 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
243 else if (j < 0 || j > n)
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
244 (*current_liboctave_error_handler) ("cholinsert: index out of range");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
246 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
247 FloatComplexColumnVector utmp = u;
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
248
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
249 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
250
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
251 chol_mat.resize (n+1, n+1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
252
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
253 F77_XFCN (cchinx, CCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
254 j + 1, utmp.fortran_vec (), rw, info));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
255 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
256
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
258 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
259
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
260 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
261 FloatComplexCHOL::delete_sym (octave_idx_type j)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
263 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
264
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
265 if (j < 0 || j > n-1)
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
266 (*current_liboctave_error_handler) ("choldelete: index out of range");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
268 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
269 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
270
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
271 F77_XFCN (cchdex, CCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
272 j + 1, rw));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
273
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
274 chol_mat.resize (n-1, n-1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
275 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
276 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
277
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
278 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
279 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
280 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
281 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
282
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
283 if (i < 0 || i > n-1 || j < 0 || j > n-1)
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
284 (*current_liboctave_error_handler) ("cholshift: index out of range");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
285 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
286 {
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
287 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
288 OCTAVE_LOCAL_BUFFER (float, rw, n);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
289
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
290 F77_XFCN (cchshx, CCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
291 i + 1, j + 1, w, rw));
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
292 }
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
293 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
294
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
295 #else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
296
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
297 void
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
298 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
299 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
300 warn_qrupdate_once ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
301
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
302 octave_idx_type n = chol_mat.rows ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
303
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
304 if (u.length () == n)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
305 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
306 init (chol_mat.hermitian () * chol_mat
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
307 + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
308 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
309 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
310 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
311 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
312
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
313 static bool
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
314 singular (const FloatComplexMatrix& a)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
315 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
316 for (octave_idx_type i = 0; i < a.rows (); i++)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
317 if (a(i,i) == 0.0f) return true;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
318 return false;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
319 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
320
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
321 octave_idx_type
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
322 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
323 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
324 warn_qrupdate_once ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
325
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
326 octave_idx_type info = -1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
327
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
328 octave_idx_type n = chol_mat.rows ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
329
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
330 if (u.length () == n)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
331 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
332 if (singular (chol_mat))
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
333 info = 2;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
334 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
335 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
336 info = init (chol_mat.hermitian () * chol_mat
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
337 - FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (), false);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
338 if (info) info = 1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
339 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
340 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
341 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
342 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
343
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
344 return info;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
345 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
346
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
347 octave_idx_type
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
348 FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u, octave_idx_type j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
349 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
350 warn_qrupdate_once ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
351
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
352 octave_idx_type info = -1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
353
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
354 octave_idx_type n = chol_mat.rows ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
355
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
356 if (u.length () != n + 1)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
357 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
358 else if (j < 0 || j > n)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
359 (*current_liboctave_error_handler) ("cholinsert: index out of range");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
360 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
361 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
362 if (singular (chol_mat))
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
363 info = 2;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
364 else if (u(j).imag () != 0.0)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
365 info = 3;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
366 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
367 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
368 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
369 FloatComplexMatrix a1 (n+1, n+1);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
370 for (octave_idx_type k = 0; k < n+1; k++)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
371 for (octave_idx_type l = 0; l < n+1; l++)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
372 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
373 if (l == j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
374 a1(k, l) = u(k);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
375 else if (k == j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
376 a1(k, l) = std::conj (u(l));
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
377 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
378 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
379 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
380 info = init (a1, false);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
381 if (info) info = 1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
382 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
383 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
384
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
385 return info;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
386 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
387
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
388 void
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
389 FloatComplexCHOL::delete_sym (octave_idx_type j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
390 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
391 warn_qrupdate_once ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
392
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
393 octave_idx_type n = chol_mat.rows ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
394
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
395 if (j < 0 || j > n-1)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
396 (*current_liboctave_error_handler) ("choldelete: index out of range");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
397 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
398 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
399 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
400 a.delete_elements (1, idx_vector (j));
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
401 a.delete_elements (0, idx_vector (j));
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
402 init (a, false);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
403 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
404 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
405
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
406 void
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
407 FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
408 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
409 warn_qrupdate_once ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
410
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
411 octave_idx_type n = chol_mat.rows ();
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
412
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
413 if (i < 0 || i > n-1 || j < 0 || j > n-1)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
414 (*current_liboctave_error_handler) ("cholshift: index out of range");
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
415 else
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
416 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
417 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
418 Array<octave_idx_type> p (n, 1);
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
419 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
420 if (i < j)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
421 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
422 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
423 p(j) = i;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
424 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
425 else if (j < i)
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
426 {
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
427 p(j) = i;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
428 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
429 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
430
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
431 init (a.index (idx_vector (p), idx_vector (p)), false);
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
432 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
433 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
434
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
435 #endif
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
436
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
437 FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
438 chol2inv (const FloatComplexMatrix& r)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
439 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
440 return chol2inv_internal (r);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
441 }