GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fCmplxLU.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2015 John W. Eaton
4 Copyright (C) 2009 VZLU Prague, a.s.
5 
6 This file is part of Octave.
7 
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include "fCmplxLU.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 #include "fCColVector.h"
33 
34 // Instantiate the base LU class for the types we need.
35 
36 #include "base-lu.h"
37 #include "base-lu.cc"
38 
39 template class base_lu <FloatComplexMatrix>;
40 
41 // Define the constructor for this particular derivation.
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (cgetrf, CGETRF) (const octave_idx_type&, const octave_idx_type&,
47  FloatComplex*, const octave_idx_type&,
48  octave_idx_type*, octave_idx_type&);
49 
50 #ifdef HAVE_QRUPDATE_LUU
51  F77_RET_T
52  F77_FUNC (clu1up, CLU1UP) (const octave_idx_type&, const octave_idx_type&,
53  FloatComplex *, const octave_idx_type&,
54  FloatComplex *, const octave_idx_type&,
55  FloatComplex *, FloatComplex *);
56 
57  F77_RET_T
58  F77_FUNC (clup1up, CLUP1UP) (const octave_idx_type&, const octave_idx_type&,
59  FloatComplex *, const octave_idx_type&,
60  FloatComplex *, const octave_idx_type&,
61  octave_idx_type *, const FloatComplex *,
62  const FloatComplex *, FloatComplex *);
63 #endif
64 }
65 
67 {
68  octave_idx_type a_nr = a.rows ();
69  octave_idx_type a_nc = a.cols ();
70  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
71 
72  ipvt.resize (dim_vector (mn, 1));
73  octave_idx_type *pipvt = ipvt.fortran_vec ();
74 
75  a_fact = a;
76  FloatComplex *tmp_data = a_fact.fortran_vec ();
77 
78  octave_idx_type info = 0;
79 
80  F77_XFCN (cgetrf, CGETRF, (a_nr, a_nc, tmp_data, a_nr, pipvt, info));
81 
82  for (octave_idx_type i = 0; i < mn; i++)
83  pipvt[i] -= 1;
84 }
85 
86 #ifdef HAVE_QRUPDATE_LUU
87 
89  const FloatComplexColumnVector& v)
90 {
91  if (packed ())
92  unpack ();
93 
96 
97  octave_idx_type m = l.rows ();
98  octave_idx_type n = r.columns ();
99  octave_idx_type k = l.columns ();
100 
101  if (u.length () == m && v.length () == n)
102  {
103  FloatComplexColumnVector utmp = u;
104  FloatComplexColumnVector vtmp = v;
105  F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (), m, r.fortran_vec (), k,
106  utmp.fortran_vec (), vtmp.fortran_vec ()));
107  }
108  else
109  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
110 }
111 
113  const FloatComplexMatrix& v)
114 {
115  if (packed ())
116  unpack ();
117 
120 
121  octave_idx_type m = l.rows ();
122  octave_idx_type n = r.columns ();
123  octave_idx_type k = l.columns ();
124 
125  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
126  {
127  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
128  {
129  FloatComplexColumnVector utmp = u.column (i);
130  FloatComplexColumnVector vtmp = v.column (i);
131  F77_XFCN (clu1up, CLU1UP, (m, n, l.fortran_vec (),
132  m, r.fortran_vec (), k,
133  utmp.fortran_vec (), vtmp.fortran_vec ()));
134  }
135  }
136  else
137  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
138 }
139 
141  const FloatComplexColumnVector& v)
142 {
143  if (packed ())
144  unpack ();
145 
148 
149  octave_idx_type m = l.rows ();
150  octave_idx_type n = r.columns ();
151  octave_idx_type k = l.columns ();
152 
153  if (u.length () == m && v.length () == n)
154  {
155  FloatComplexColumnVector utmp = u;
156  FloatComplexColumnVector vtmp = v;
158  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
159  F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (),
160  m, r.fortran_vec (), k,
161  ipvt.fortran_vec (),
162  utmp.data (), vtmp.data (), w));
163  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
164  }
165  else
166  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
167 }
168 
170  const FloatComplexMatrix& v)
171 {
172  if (packed ())
173  unpack ();
174 
177 
178  octave_idx_type m = l.rows ();
179  octave_idx_type n = r.columns ();
180  octave_idx_type k = l.columns ();
181 
182  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
183  {
185  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
186  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
187  {
188  FloatComplexColumnVector utmp = u.column (i);
189  FloatComplexColumnVector vtmp = v.column (i);
190  F77_XFCN (clup1up, CLUP1UP, (m, n, l.fortran_vec (),
191  m, r.fortran_vec (), k,
192  ipvt.fortran_vec (),
193  utmp.data (), vtmp.data (), w));
194  }
195  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
196  }
197  else
198  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
199 }
200 
201 #else
202 
205 {
206  (*current_liboctave_error_handler)
207  ("luupdate: not available in this version");
208 }
209 
211  const FloatComplexMatrix&)
212 {
213  (*current_liboctave_error_handler)
214  ("luupdate: not available in this version");
215 }
216 
219 {
220  (*current_liboctave_error_handler)
221  ("luupdate: not available in this version");
222 }
223 
225  const FloatComplexMatrix&)
226 {
227  (*current_liboctave_error_handler)
228  ("luupdate: not available in this version");
229 }
230 
231 #endif
FloatComplexLU(void)
Definition: fCmplxLU.h:36
FloatComplexMatrix l_fact
Definition: base-lu.h:82
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:51
octave_idx_type rows(void) const
Definition: Array.h:313
liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:1002
FloatComplexMatrix a_fact
Definition: base-lu.h:81
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
Definition: Array.h:479
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1033
#define F77_RET_T
Definition: f77-fcn.h:264
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
void update_piv(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
Definition: fCmplxLU.cc:140
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
std::complex< float > FloatComplex
Definition: oct-cmplx.h:30
void update(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
Definition: fCmplxLU.cc:88
const T * fortran_vec(void) const
Definition: Array.h:481
Array< octave_idx_type > ipvt
Definition: base-lu.h:84
octave_idx_type cols(void) const
Definition: Array.h:321
octave_idx_type columns(void) const
Definition: Array.h:322
F77_RET_T F77_FUNC(cgetrf, CGETRF)(const octave_idx_type &