Mercurial > octave-nkf

annotate liboctave/fCmplxCHOL.cc @ 8547:d66c9b6e506a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
imported patch qrupdate.diff
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 20 Jan 2009 21:16:42 +0100
parents 25bc2d31e1bf
children a6edd5c23cb5
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"
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
37
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
38 extern "C"
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
39 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
40 F77_RET_T
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
41 F77_FUNC (cpotrf, CPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
42 FloatComplex*, const octave_idx_type&, octave_idx_type&
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
43 F77_CHAR_ARG_LEN_DECL);
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 (cpotri, CPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
46 FloatComplex*, const octave_idx_type&, octave_idx_type&
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
47 F77_CHAR_ARG_LEN_DECL);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
48
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
49 F77_RET_T
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
50 F77_FUNC (cpocon, CPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
51 FloatComplex*, const octave_idx_type&, const float&,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
52 float&, FloatComplex*, float*,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
53 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
54 #ifdef HAVE_QRUPDATE
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
55
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
56 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
57 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
58 FloatComplex*, float*);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
59
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 (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
62 FloatComplex*, float*, octave_idx_type&);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
63
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
64 F77_RET_T
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
65 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
66 const octave_idx_type&, FloatComplex*, float*,
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
67 octave_idx_type&);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
68
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
69 F77_RET_T
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
70 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
71 const octave_idx_type&, float*);
7789
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 (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
75 const octave_idx_type&, const octave_idx_type&,
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
76 FloatComplex*, float*);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
77 #endif
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
78 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
79
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
80 octave_idx_type
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
81 FloatComplexCHOL::init (const FloatComplexMatrix& a, bool calc_cond)
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 octave_idx_type a_nr = a.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
84 octave_idx_type a_nc = a.cols ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
85
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
86 if (a_nr != a_nc)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
87 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
88 (*current_liboctave_error_handler)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
89 ("FloatComplexCHOL requires square matrix");
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
90 return -1;
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
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
93 octave_idx_type n = a_nc;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
94 octave_idx_type info;
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 chol_mat = a;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
97 FloatComplex *h = chol_mat.fortran_vec ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
98
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
99 // Calculate the norm of the matrix, for later use.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
100 float anorm = 0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
101 if (calc_cond)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
102 anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
103
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
104 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
105 F77_CHAR_ARG_LEN (1)));
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
106
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
107 xrcond = 0.0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
108 if (info != 0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
109 info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
110 else if (calc_cond)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
111 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
112 octave_idx_type cpocon_info = 0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
113
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
114 // Now calculate the condition number for non-singular matrix.
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
115 Array<FloatComplex> z (2*n);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
116 FloatComplex *pz = z.fortran_vec ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
117 Array<float> rz (n);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
118 float *prz = rz.fortran_vec ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
119 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
120 n, anorm, xrcond, pz, prz, cpocon_info
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
121 F77_CHAR_ARG_LEN (1)));
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 if (cpocon_info != 0)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
124 info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
125 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
126 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
127 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
128 // 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
129 // faster for that matter :-)), please let me know!
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
130
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
131 if (n > 1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
132 for (octave_idx_type j = 0; j < a_nc; j++)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
133 for (octave_idx_type i = j+1; i < a_nr; i++)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
134 chol_mat.xelem (i, j) = 0.0;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
135 }
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 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
138 }
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 static FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
141 chol2inv_internal (const FloatComplexMatrix& r)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
142 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
143 FloatComplexMatrix retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
144
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
145 octave_idx_type r_nr = r.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
146 octave_idx_type r_nc = r.cols ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
147
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
148 if (r_nr == r_nc)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
149 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
150 octave_idx_type n = r_nc;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
151 octave_idx_type info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
152
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
153 FloatComplexMatrix tmp = r;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
154
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
155 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
156 tmp.fortran_vec (), n, info
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
157 F77_CHAR_ARG_LEN (1)));
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
158
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
159 // 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
160 // faster for that matter :-)), please let me know!
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
161
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
162 if (n > 1)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
163 for (octave_idx_type j = 0; j < r_nc; j++)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
164 for (octave_idx_type i = j+1; i < r_nr; i++)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
165 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
166
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
167 retval = tmp;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
168 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
169 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
170 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
171
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
172 return retval;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
173 }
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 // 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
176 FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
177 FloatComplexCHOL::inverse (void) const
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
178 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
179 return chol2inv_internal (chol_mat);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
180 }
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 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
183 FloatComplexCHOL::set (const FloatComplexMatrix& R)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
184 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
185 if (R.is_square ())
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
186 chol_mat = R;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
187 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
188 (*current_liboctave_error_handler) ("CHOL requires square matrix");
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
189 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
190
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
191 #ifdef HAVE_QRUPDATE
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
192
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
193 void
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
194 FloatComplexCHOL::update (const FloatComplexColumnVector& u)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
195 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
196 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
197
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
198 if (u.length () == n)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
199 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
200 FloatComplexColumnVector utmp = u;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
201
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
202 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
203
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
204 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
205 utmp.fortran_vec (), rw));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
206 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
207 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
208 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
209 }
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 octave_idx_type
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
212 FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
213 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
214 octave_idx_type info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
215
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
216 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
217
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
218 if (u.length () == n)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
219 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
220 FloatComplexColumnVector utmp = u;
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
221
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
222 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
223
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
224 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
225 utmp.fortran_vec (), rw, info));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
226 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
227 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
228 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
229
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
230 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
231 }
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 octave_idx_type
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
234 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
235 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
236 octave_idx_type info = -1;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
237
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
238 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
239
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
240 if (u.length () != n + 1)
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
241 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
242 else if (j < 0 || j > n)
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
243 (*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
244 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
245 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
246 FloatComplexColumnVector utmp = u;
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
247
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
248 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
249
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
250 chol_mat.resize (n+1, n+1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
251
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
252 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
253 j + 1, utmp.fortran_vec (), rw, info));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
254 }
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 return info;
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
257 }
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 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
260 FloatComplexCHOL::delete_sym (octave_idx_type j)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
261 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
262 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
263
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
264 if (j < 0 || j > n-1)
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
265 (*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
266 else
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
267 {
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
268 OCTAVE_LOCAL_BUFFER (float, rw, n);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
269
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
270 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
271 j + 1, rw));
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
272
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
273 chol_mat.resize (n-1, n-1);
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
274 }
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 void
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
278 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
279 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
280 octave_idx_type n = chol_mat.rows ();
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
281
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
282 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
283 (*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
284 else
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
285 {
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
286 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
287 OCTAVE_LOCAL_BUFFER (float, rw, n);
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
288
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
289 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
290 i + 1, j + 1, w, rw));
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
291 }
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
292 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
293
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
294 #endif
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
295
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
296 FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
297 chol2inv (const FloatComplexMatrix& r)
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
298 {
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
299 return chol2inv_internal (r);
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
300 }
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
301
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
302 /*
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
303 ;;; Local Variables: ***
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
304 ;;; mode: C++ ***
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
305 ;;; End: ***
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
diff changeset
306 */