Mercurial > hg > octave-lyh

annotate liboctave/base-lu.cc @ 8920:eb63fbe60fab

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update copyright notices
author John W. Eaton <jwe@octave.org>
date Sat, 07 Mar 2009 10:41:27 -0500
parents 445d27d79f4e
children 50db3c5175b5
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
1 /*
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
2
8920
eb63fbe60fab update copyright notices
John W. Eaton <jwe@octave.org>
parents: 8367
diff changeset
3 Copyright (C) 1996, 1997, 2002, 2003, 2004, 2005, 2007, 2008 John W. Eaton
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
4
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
6
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
10 option) any later version.
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
11
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
15 for more details.
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
16
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5307
diff changeset
19 <http://www.gnu.org/licenses/>.
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
20
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
21 */
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
22
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
23 #ifdef HAVE_CONFIG_H
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
24 #include <config.h>
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
25 #endif
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
26
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
27 #include "base-lu.h"
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
28
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
29 template <class lu_type>
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
30 lu_type
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
31 base_lu <lu_type> :: L (void) const
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
32 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
33 octave_idx_type a_nr = a_fact.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
34 octave_idx_type a_nc = a_fact.cols ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
35 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
36
4329
d53c33d93440 [project @ 2003-02-18 20:00:48 by jwe]
jwe
parents: 4192
diff changeset
37 lu_type l (a_nr, mn, lu_elt_type (0.0));
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
38
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
39 for (octave_idx_type i = 0; i < a_nr; i++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
40 {
4811
198f3712c692 [project @ 2004-03-03 19:24:33 by jwe]
jwe
parents: 4329
diff changeset
41 if (i < a_nc)
198f3712c692 [project @ 2004-03-03 19:24:33 by jwe]
jwe
parents: 4329
diff changeset
42 l.xelem (i, i) = 1.0;
198f3712c692 [project @ 2004-03-03 19:24:33 by jwe]
jwe
parents: 4329
diff changeset
43
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
44 for (octave_idx_type j = 0; j < (i < a_nc ? i : a_nc); j++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
45 l.xelem (i, j) = a_fact.xelem (i, j);
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
46 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
47
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
48 return l;
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
49 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
50
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
51 template <class lu_type>
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
52 lu_type
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
53 base_lu <lu_type> :: U (void) const
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
54 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
55 octave_idx_type a_nr = a_fact.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
56 octave_idx_type a_nc = a_fact.cols ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
57 octave_idx_type mn = (a_nr < a_nc ? a_nr : a_nc);
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
58
4329
d53c33d93440 [project @ 2003-02-18 20:00:48 by jwe]
jwe
parents: 4192
diff changeset
59 lu_type u (mn, a_nc, lu_elt_type (0.0));
d53c33d93440 [project @ 2003-02-18 20:00:48 by jwe]
jwe
parents: 4192
diff changeset
60
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
61 for (octave_idx_type i = 0; i < mn; i++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
62 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
63 for (octave_idx_type j = i; j < a_nc; j++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
64 u.xelem (i, j) = a_fact.xelem (i, j);
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
65 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
66
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
67 return u;
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
68 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
69
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
70 template <class lu_type>
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
71 Array<octave_idx_type>
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
72 base_lu <lu_type> :: getp (void) const
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
73 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
74 octave_idx_type a_nr = a_fact.rows ();
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
75
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
76 Array<octave_idx_type> pvt (a_nr);
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
77
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
78 for (octave_idx_type i = 0; i < a_nr; i++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
79 pvt.xelem (i) = i;
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
80
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
81 for (octave_idx_type i = 0; i < ipvt.length(); i++)
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
82 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
83 octave_idx_type k = ipvt.xelem (i);
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
84
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
85 if (k != i)
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
86 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
87 octave_idx_type tmp = pvt.xelem (k);
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
88 pvt.xelem (k) = pvt.xelem (i);
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
89 pvt.xelem (i) = tmp;
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
90 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
91 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
92
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
93 return pvt;
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
94 }
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
95
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
96 template <class lu_type>
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
97 PermMatrix
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
98 base_lu <lu_type> :: P (void) const
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
99 {
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
100 return PermMatrix (getp (), false);
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
101 }
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
102
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
103 template <class lu_type>
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
104 ColumnVector
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
105 base_lu <lu_type> :: P_vec (void) const
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
106 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
107 octave_idx_type a_nr = a_fact.rows ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
108
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
109 ColumnVector p (a_nr);
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
110
8367
445d27d79f4e support permutation matrix objects
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
111 Array<octave_idx_type> pvt = getp ();
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
112
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
113 for (octave_idx_type i = 0; i < a_nr; i++)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
114 p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
115
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
116 return p;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
117 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
118
1991
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
119 /*
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
120 ;;; Local Variables: ***
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
121 ;;; mode: C++ ***
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
122 ;;; End: ***
413f6a81868f [project @ 1996-03-03 00:40:53 by jwe]
jwe
parents:
diff changeset
123 */