GNU Octave  3.8.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-2013 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, vtmp = v;
103  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
104  m, r.fortran_vec (), k,
105  utmp.fortran_vec (), vtmp.fortran_vec ()));
106  }
107  else
108  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
109 }
110 
111 void FloatLU::update (const FloatMatrix& u, const FloatMatrix& v)
112 {
113  if (packed ())
114  unpack ();
115 
116  FloatMatrix& l = l_fact;
117  FloatMatrix& r = a_fact;
118 
119  octave_idx_type m = l.rows ();
120  octave_idx_type n = r.columns ();
121  octave_idx_type k = l.columns ();
122 
123  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
124  {
125  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
126  {
127  FloatColumnVector utmp = u.column (i), vtmp = v.column (i);
128  F77_XFCN (slu1up, SLU1UP, (m, n, l.fortran_vec (),
129  m, r.fortran_vec (), k,
130  utmp.fortran_vec (), vtmp.fortran_vec ()));
131  }
132  }
133  else
134  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
135 }
136 
138  const FloatColumnVector& v)
139 {
140  if (packed ())
141  unpack ();
142 
143  FloatMatrix& l = l_fact;
144  FloatMatrix& r = a_fact;
145 
146  octave_idx_type m = l.rows ();
147  octave_idx_type n = r.columns ();
148  octave_idx_type k = l.columns ();
149 
150  if (u.length () == m && v.length () == n)
151  {
152  FloatColumnVector utmp = u, vtmp = v;
153  OCTAVE_LOCAL_BUFFER (float, w, m);
154  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
155  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
156  m, r.fortran_vec (), k,
157  ipvt.fortran_vec (),
158  utmp.data (), vtmp.data (), w));
159  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
160  }
161  else
162  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
163 }
164 
165 void FloatLU::update_piv (const FloatMatrix& u, const FloatMatrix& v)
166 {
167  if (packed ())
168  unpack ();
169 
170  FloatMatrix& l = l_fact;
171  FloatMatrix& r = a_fact;
172 
173  octave_idx_type m = l.rows ();
174  octave_idx_type n = r.columns ();
175  octave_idx_type k = l.columns ();
176 
177  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
178  {
179  OCTAVE_LOCAL_BUFFER (float, w, m);
180  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
181  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
182  {
183  FloatColumnVector utmp = u.column (i), vtmp = v.column (i);
184  F77_XFCN (slup1up, SLUP1UP, (m, n, l.fortran_vec (),
185  m, r.fortran_vec (), k,
186  ipvt.fortran_vec (),
187  utmp.data (), vtmp.data (), w));
188  }
189  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
190  }
191  else
192  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
193 }
194 
195 #else
196 
198 {
199  (*current_liboctave_error_handler)
200  ("luupdate: not available in this version");
201 }
202 
203 void FloatLU::update (const FloatMatrix&, const FloatMatrix&)
204 {
205  (*current_liboctave_error_handler)
206  ("luupdate: not available in this version");
207 }
208 
210 {
211  (*current_liboctave_error_handler)
212  ("luupdate: not available in this version");
213 }
214 
215 void FloatLU::update_piv (const FloatMatrix&, const FloatMatrix&)
216 {
217  (*current_liboctave_error_handler)
218  ("luupdate: not available in this version");
219 }
220 
221 #endif