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
sparse-base-lu.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2013 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
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 "sparse-base-lu.h"
29 
30 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
31 lu_type
33 {
34  octave_idx_type nr = Lfact.rows ();
35  octave_idx_type nc = Ufact.rows ();
36  octave_idx_type rcmin = (nr > nc ? nr : nc);
37 
38  lu_type Yout (nr, nc, Lfact.nnz () + Ufact.nnz ());
39  octave_idx_type ii = 0;
40  Yout.xcidx (0) = 0;
41 
42  for (octave_idx_type j = 0; j < nc; j++)
43  {
44  for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx (j + 1); i++)
45  {
46  Yout.xridx (ii) = Ufact.ridx (i);
47  Yout.xdata (ii++) = Ufact.data (i);
48  }
49  if (j < rcmin)
50  {
51  // Note the +1 skips the 1.0 on the diagonal
52  for (octave_idx_type i = Lfact.cidx (j) + 1;
53  i < Lfact.cidx (j +1); i++)
54  {
55  Yout.xridx (ii) = Lfact.ridx (i);
56  Yout.xdata (ii++) = Lfact.data (i);
57  }
58  }
59  Yout.xcidx (j + 1) = ii;
60  }
61 
62  return Yout;
63 }
64 
65 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
66 p_type
68 {
69 
70  octave_idx_type nr = Lfact.rows ();
71 
72  p_type Pout (nr, nr, nr);
73 
74  for (octave_idx_type i = 0; i < nr; i++)
75  {
76  Pout.cidx (i) = i;
77  Pout.ridx (P (i)) = i;
78  Pout.data (i) = 1;
79  }
80  Pout.cidx (nr) = nr;
81 
82  return Pout;
83 }
84 
85 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
88 {
89 
90  octave_idx_type nr = Lfact.rows ();
91 
92  ColumnVector Pout (nr);
93 
94  for (octave_idx_type i = 0; i < nr; i++)
95  Pout.xelem (i) = static_cast<double> (P(i) + 1);
96 
97  return Pout;
98 }
99 
100 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
103 {
104  return PermMatrix (P, false);
105 }
106 
107 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
108 p_type
110 {
111  octave_idx_type nc = Ufact.cols ();
112 
113  p_type Pout (nc, nc, nc);
114 
115  for (octave_idx_type i = 0; i < nc; i++)
116  {
117  Pout.cidx (i) = i;
118  Pout.ridx (i) = Q (i);
119  Pout.data (i) = 1;
120  }
121  Pout.cidx (nc) = nc;
122 
123  return Pout;
124 }
125 
126 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
129 {
130 
131  octave_idx_type nc = Ufact.cols ();
132 
133  ColumnVector Pout (nc);
134 
135  for (octave_idx_type i = 0; i < nc; i++)
136  Pout.xelem (i) = static_cast<double> (Q(i) + 1);
137 
138  return Pout;
139 }
140 
141 template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
144 {
145  return PermMatrix (Q, true);
146 }