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
floatLU.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 "floatLU.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 #include "fColVector.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 <FloatMatrix>;
40 
41 // Define the constructor for this particular derivation.
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (sgetrf, SGETRF) (const octave_idx_type&, const octave_idx_type&,
47  float*, const octave_idx_type&, octave_idx_type*,
48  octave_idx_type&);
49 
50 #ifdef HAVE_QRUPDATE_LUU
51  F77_RET_T
52  F77_FUNC (slu1up, SLU1UP) (const octave_idx_type&, const octave_idx_type&,
53  float *, const octave_idx_type&,
54  float *, const octave_idx_type&,
55  float *, float *);
56 
57  F77_RET_T
58  F77_FUNC (slup1up, SLUP1UP) (const octave_idx_type&, const octave_idx_type&,
59  float *, const octave_idx_type&,
60  float *, const octave_idx_type&,
61  octave_idx_type *, const float *,
62  const float *, float *);
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  float *tmp_data = a_fact.fortran_vec ();
77 
78  octave_idx_type info = 0;
79 
80  F77_XFCN (sgetrf, SGETRF, (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 {
90  if (packed ())
91  unpack ();
92 
93  FloatMatrix& l = l_fact;
94  FloatMatrix& r = a_fact;
95 
96  octave_idx_type m = l.rows ();
97  octave_idx_type n = r.columns ();
98  octave_idx_type k = l.columns ();
99 
100  if (u.length () == m && v.length () == n)
101  {
102  FloatColumnVector utmp = u;
103  FloatColumnVector vtmp = v;
104  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
105  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 
112 void FloatLU::update (const FloatMatrix& u, const FloatMatrix& v)
113 {
114  if (packed ())
115  unpack ();
116 
117  FloatMatrix& l = l_fact;
118  FloatMatrix& 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  FloatColumnVector utmp = u.column (i);
129  FloatColumnVector vtmp = v.column (i);
130  F77_XFCN (slu1up, SLU1UP, (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 FloatColumnVector& v)
141 {
142  if (packed ())
143  unpack ();
144 
145  FloatMatrix& l = l_fact;
146  FloatMatrix& 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  FloatColumnVector utmp = u;
155  FloatColumnVector vtmp = v;
156  OCTAVE_LOCAL_BUFFER (float, w, m);
157  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
158  F77_XFCN (slup1up, SLUP1UP, (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 
168 void FloatLU::update_piv (const FloatMatrix& u, const FloatMatrix& v)
169 {
170  if (packed ())
171  unpack ();
172 
173  FloatMatrix& l = l_fact;
174  FloatMatrix& 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  {
182  OCTAVE_LOCAL_BUFFER (float, w, m);
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  FloatColumnVector utmp = u.column (i);
187  FloatColumnVector vtmp = v.column (i);
188  F77_XFCN (slup1up, SLUP1UP, (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 FloatLU::update (const FloatMatrix&, const FloatMatrix&)
208 {
209  (*current_liboctave_error_handler)
210  ("luupdate: not available in this version");
211 }
212 
214 {
215  (*current_liboctave_error_handler)
216  ("luupdate: not available in this version");
217 }
218 
219 void FloatLU::update_piv (const FloatMatrix&, const FloatMatrix&)
220 {
221  (*current_liboctave_error_handler)
222  ("luupdate: not available in this version");
223 }
224 
225 #endif
FloatMatrix 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
bool packed(void) const
FloatLU(void)
Definition: floatLU.h:36
FloatMatrix 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
void update_piv(const FloatColumnVector &u, const FloatColumnVector &v)
Definition: floatLU.cc:139
F77_RET_T F77_FUNC(sgetrf, SGETRF)(const octave_idx_type &
#define F77_RET_T
Definition: f77-fcn.h:264
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:652
void update(const FloatColumnVector &u, const FloatColumnVector &v)
Definition: floatLU.cc:88
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
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