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
base-lu.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 John W. Eaton
4 Copyright (C) 2009 VZLU Prague
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 "base-lu.h"
29 
30 template <class lu_type>
31 base_lu<lu_type>::base_lu (const lu_type& l, const lu_type& u,
32  const PermMatrix& p)
33  : a_fact (u), l_fact (l), ipvt (p.pvec ())
34 {
35  if (l.columns () != u.rows ())
36  (*current_liboctave_error_handler) ("lu: dimension mismatch");
37 }
38 
39 template <class lu_type>
40 bool
41 base_lu <lu_type> :: packed (void) const
42 {
43  return l_fact.dims () == dim_vector ();
44 }
45 
46 template <class lu_type>
47 void
48 base_lu <lu_type> :: unpack (void)
49 {
50  if (packed ())
51  {
52  l_fact = L ();
53  a_fact = U (); // FIXME: sub-optimal
54  ipvt = getp ();
55  }
56 }
57 
58 template <class lu_type>
59 lu_type
60 base_lu <lu_type> :: L (void) const
61 {
62  if (packed ())
63  {
64  octave_idx_type a_nr = a_fact.rows ();
65  octave_idx_type a_nc = a_fact.cols ();
66  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
67 
68  lu_type l (a_nr, mn, lu_elt_type (0.0));
69 
70  for (octave_idx_type i = 0; i < a_nr; i++)
71  {
72  if (i < a_nc)
73  l.xelem (i, i) = 1.0;
74 
75  for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
76  l.xelem (i, j) = a_fact.xelem (i, j);
77  }
78 
79  return l;
80  }
81  else
82  return l_fact;
83 }
84 
85 template <class lu_type>
86 lu_type
87 base_lu <lu_type> :: U (void) const
88 {
89  if (packed ())
90  {
91  octave_idx_type a_nr = a_fact.rows ();
92  octave_idx_type a_nc = a_fact.cols ();
93  octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
94 
95  lu_type u (mn, a_nc, lu_elt_type (0.0));
96 
97  for (octave_idx_type i = 0; i < mn; i++)
98  {
99  for (octave_idx_type j = i; j < a_nc; j++)
100  u.xelem (i, j) = a_fact.xelem (i, j);
101  }
102 
103  return u;
104  }
105  else
106  return a_fact;
107 }
108 
109 template <class lu_type>
110 lu_type
111 base_lu <lu_type> :: Y (void) const
112 {
113  if (! packed ())
114  (*current_liboctave_error_handler)
115  ("lu: Y () not implemented for unpacked form");
116  return a_fact;
117 }
118 
119 template <class lu_type>
121 base_lu <lu_type> :: getp (void) const
122 {
123  if (packed ())
124  {
125  octave_idx_type a_nr = a_fact.rows ();
126 
127  Array<octave_idx_type> pvt (dim_vector (a_nr, 1));
128 
129  for (octave_idx_type i = 0; i < a_nr; i++)
130  pvt.xelem (i) = i;
131 
132  for (octave_idx_type i = 0; i < ipvt.length (); i++)
133  {
134  octave_idx_type k = ipvt.xelem (i);
135 
136  if (k != i)
137  {
138  octave_idx_type tmp = pvt.xelem (k);
139  pvt.xelem (k) = pvt.xelem (i);
140  pvt.xelem (i) = tmp;
141  }
142  }
143 
144  return pvt;
145  }
146  else
147  return ipvt;
148 }
149 
150 template <class lu_type>
152 base_lu <lu_type> :: P (void) const
153 {
154  return PermMatrix (getp (), false);
155 }
156 
157 template <class lu_type>
159 base_lu <lu_type> :: P_vec (void) const
160 {
161  octave_idx_type a_nr = a_fact.rows ();
162 
163  ColumnVector p (a_nr);
164 
165  Array<octave_idx_type> pvt = getp ();
166 
167  for (octave_idx_type i = 0; i < a_nr; i++)
168  p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
169 
170  return p;
171 }
172 
173 template <class lu_type>
174 bool
176 {
177  bool retval = true;
178 
179  octave_idx_type k = std::min (a_fact.rows (), a_fact.columns ());
180 
181  for (octave_idx_type i = 0; i < k; i++)
182  {
183  if (a_fact(i, i) == lu_elt_type ())
184  {
185  retval = false;
186  break;
187  }
188  }
189 
190  return retval;
191 }