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
dbleLU.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 "dbleLU.h"
29 #include "f77-fcn.h"
30 #include "lo-error.h"
31 #include "oct-locbuf.h"
32 #include "dColVector.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 <Matrix>;
40 
41 // Define the constructor for this particular derivation.
42 
43 extern "C"
44 {
45  F77_RET_T
46  F77_FUNC (dgetrf, DGETRF) (const octave_idx_type&, const octave_idx_type&,
47  double*, const octave_idx_type&,
48  octave_idx_type*, octave_idx_type&);
49 
50 #ifdef HAVE_QRUPDATE_LUU
51  F77_RET_T
52  F77_FUNC (dlu1up, DLU1UP) (const octave_idx_type&, const octave_idx_type&,
53  double *, const octave_idx_type&,
54  double *, const octave_idx_type&,
55  double *, double *);
56 
57  F77_RET_T
58  F77_FUNC (dlup1up, DLUP1UP) (const octave_idx_type&, const octave_idx_type&,
59  double *, const octave_idx_type&,
60  double *, const octave_idx_type&,
61  octave_idx_type *, const double *,
62  const double *, double *);
63 #endif
64 }
65 
66 LU::LU (const Matrix& a)
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  double *tmp_data = a_fact.fortran_vec ();
77 
78  octave_idx_type info = 0;
79 
80  F77_XFCN (dgetrf, DGETRF, (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 
88 void LU::update (const ColumnVector& u, const ColumnVector& v)
89 {
90  if (packed ())
91  unpack ();
92 
93  Matrix& l = l_fact;
94  Matrix& 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  ColumnVector utmp = u;
103  ColumnVector vtmp = v;
104  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (), 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 LU::update (const Matrix& u, const Matrix& v)
112 {
113  if (packed ())
114  unpack ();
115 
116  Matrix& l = l_fact;
117  Matrix& 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  ColumnVector utmp = u.column (i);
128  ColumnVector vtmp = v.column (i);
129  F77_XFCN (dlu1up, DLU1UP, (m, n, l.fortran_vec (),
130  m, r.fortran_vec (), k,
131  utmp.fortran_vec (), vtmp.fortran_vec ()));
132  }
133  }
134  else
135  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
136 }
137 
138 void LU::update_piv (const ColumnVector& u, const ColumnVector& v)
139 {
140  if (packed ())
141  unpack ();
142 
143  Matrix& l = l_fact;
144  Matrix& 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  ColumnVector utmp = u;
153  ColumnVector vtmp = v;
154  OCTAVE_LOCAL_BUFFER (double, w, m);
155  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
156  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
157  m, r.fortran_vec (), k,
158  ipvt.fortran_vec (),
159  utmp.data (), vtmp.data (), w));
160  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
161  }
162  else
163  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
164 }
165 
166 void LU::update_piv (const Matrix& u, const Matrix& v)
167 {
168  if (packed ())
169  unpack ();
170 
171  Matrix& l = l_fact;
172  Matrix& r = a_fact;
173 
174  octave_idx_type m = l.rows ();
175  octave_idx_type n = r.columns ();
176  octave_idx_type k = l.columns ();
177 
178  if (u.rows () == m && v.rows () == n && u.cols () == v.cols ())
179  {
180  OCTAVE_LOCAL_BUFFER (double, w, m);
181  for (octave_idx_type i = 0; i < m; i++) ipvt(i) += 1; // increment
182  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
183  {
184  ColumnVector utmp = u.column (i);
185  ColumnVector vtmp = v.column (i);
186  F77_XFCN (dlup1up, DLUP1UP, (m, n, l.fortran_vec (),
187  m, r.fortran_vec (), k,
188  ipvt.fortran_vec (),
189  utmp.data (), vtmp.data (), w));
190  }
191  for (octave_idx_type i = 0; i < m; i++) ipvt(i) -= 1; // decrement
192  }
193  else
194  (*current_liboctave_error_handler) ("luupdate: dimensions mismatch");
195 }
196 
197 #else
198 
199 void LU::update (const ColumnVector&, const ColumnVector&)
200 {
201  (*current_liboctave_error_handler)
202  ("luupdate: not available in this version");
203 }
204 
205 void LU::update (const Matrix&, const Matrix&)
206 {
207  (*current_liboctave_error_handler)
208  ("luupdate: not available in this version");
209 }
210 
211 void LU::update_piv (const ColumnVector&, const ColumnVector&)
212 {
213  (*current_liboctave_error_handler)
214  ("luupdate: not available in this version");
215 }
216 
217 void LU::update_piv (const Matrix&, const Matrix&)
218 {
219  (*current_liboctave_error_handler)
220  ("luupdate: not available in this version");
221 }
222 
223 #endif
Matrix l_fact
Definition: base-lu.h:82
void update(const ColumnVector &u, const ColumnVector &v)
Definition: dbleLU.cc:88
#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
Matrix 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
Definition: dMatrix.h:35
LU(void)
Definition: dbleLU.h:35
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
F77_RET_T F77_FUNC(dgetrf, DGETRF)(const octave_idx_type &
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:645
void update_piv(const ColumnVector &u, const ColumnVector &v)
Definition: dbleLU.cc:138
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
void unpack(void)
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