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
CmplxLU.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 "CmplxLU.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 #include "CColVector.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 <ComplexMatrix>;
40 
41 // Define the constructor for this particular derivation.
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (zgetrf, ZGETRF) (const octave_idx_type&, const octave_idx_type&,
47  Complex*, const octave_idx_type&,
48  octave_idx_type*, octave_idx_type&);
49 
50 #ifdef HAVE_QRUPDATE_LUU
51  F77_RET_T
52  F77_FUNC (zlu1up, ZLU1UP) (const octave_idx_type&, const octave_idx_type&,
53  Complex *, const octave_idx_type&,
54  Complex *, const octave_idx_type&,
55  Complex *, Complex *);
56 
57  F77_RET_T
58  F77_FUNC (zlup1up, ZLUP1UP) (const octave_idx_type&, const octave_idx_type&,
59  Complex *, const octave_idx_type&,
60  Complex *, const octave_idx_type&,
61  octave_idx_type *, const Complex *,
62  const Complex *, Complex *);
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  Complex *tmp_data = a_fact.fortran_vec ();
77 
78  octave_idx_type info = 0;
79 
80  F77_XFCN (zgetrf, ZGETRF, (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 ComplexColumnVector& v)
90 {
91  if (packed ())
92  unpack ();
93 
94  ComplexMatrix& l = l_fact;
95  ComplexMatrix& r = a_fact;
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  ComplexColumnVector utmp = u;
104  ComplexColumnVector vtmp = v;
105  F77_XFCN (zlu1up, ZLU1UP, (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 {
114  if (packed ())
115  unpack ();
116 
117  ComplexMatrix& l = l_fact;
118  ComplexMatrix& r = a_fact;
119 
120  octave_idx_type m = l.rows ();
121  octave_idx_type n = r.columns ();
122  octave_idx_type k = l.columns ();
123 
124  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
125  {
126  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
127  {
128  ComplexColumnVector utmp = u.column (i);
129  ComplexColumnVector vtmp = v.column (i);
130  F77_XFCN (zlu1up, ZLU1UP, (m, n, l.fortran_vec (),
131  m, r.fortran_vec (), k,
132  utmp.fortran_vec (), vtmp.fortran_vec ()));
133  }
134  }
135  else
136  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
137 }
138 
140  const ComplexColumnVector& v)
141 {
142  if (packed ())
143  unpack ();
144 
145  ComplexMatrix& l = l_fact;
146  ComplexMatrix& r = a_fact;
147 
148  octave_idx_type m = l.rows ();
149  octave_idx_type n = r.columns ();
150  octave_idx_type k = l.columns ();
151 
152  if (u.length () == m && v.length () == n)
153  {
154  ComplexColumnVector utmp = u;
155  ComplexColumnVector vtmp = v;
157  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
158  F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (),
159  m, r.fortran_vec (), k,
160  ipvt.fortran_vec (),
161  utmp.data (), vtmp.data (), w));
162  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
163  }
164  else
165  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
166 }
167 
169 {
170  if (packed ())
171  unpack ();
172 
173  ComplexMatrix& l = l_fact;
174  ComplexMatrix& r = a_fact;
175 
176  octave_idx_type m = l.rows ();
177  octave_idx_type n = r.columns ();
178  octave_idx_type k = l.columns ();
179 
180  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
181  {
183  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
184  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
185  {
186  ComplexColumnVector utmp = u.column (i);
187  ComplexColumnVector vtmp = v.column (i);
188  F77_XFCN (zlup1up, ZLUP1UP, (m, n, l.fortran_vec (),
189  m, r.fortran_vec (), k,
190  ipvt.fortran_vec (),
191  utmp.data (), vtmp.data (), w));
192  }
193  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
194  }
195  else
196  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
197 }
198 
199 #else
200 
202 {
203  (*current_liboctave_error_handler)
204  ("luupdate: not available in this version");
205 }
206 
207 void ComplexLU::update (const ComplexMatrix&, const ComplexMatrix&)
208 {
209  (*current_liboctave_error_handler)
210  ("luupdate: not available in this version");
211 }
212 
214  const ComplexColumnVector&)
215 {
216  (*current_liboctave_error_handler)
217  ("luupdate: not available in this version");
218 }
219 
221 {
222  (*current_liboctave_error_handler)
223  ("luupdate: not available in this version");
224 }
225 
226 #endif
void update(const ComplexColumnVector &u, const ComplexColumnVector &v)
Definition: CmplxLU.cc:88
F77_RET_T F77_FUNC(zgetrf, ZGETRF)(const octave_idx_type &
ComplexLU(void)
Definition: CmplxLU.h:36
ComplexMatrix l_fact
Definition: base-lu.h:82
void update_piv(const ComplexColumnVector &u, const ComplexColumnVector &v)
Definition: CmplxLU.cc:139
#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
bool packed(void) const
ComplexMatrix 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
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
std::complex< double > Complex
Definition: oct-cmplx.h:29
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
ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:996
octave_idx_type columns(void) const
Definition: Array.h:322