annotate liboctave/sparse-dmsolve.cc @ 7334:274d9642ac10 release-3-0-0

[project @ 2007-12-22 01:42:45 by jwe]
author jwe
date Sat, 22 Dec 2007 01:42:46 +0000
parents a1dbe9d80eee
children b166043585a8
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
1 /*
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
2
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
3 Copyright (C) 2006, 2007 David Bateman
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
4
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6924
diff changeset
5 This file is part of Octave.
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6924
diff changeset
6
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
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: 6924
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: 6924
diff changeset
10 option) any later version.
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
11
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
15 for more details.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
16
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
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: 6924
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: 6924
diff changeset
19 <http://www.gnu.org/licenses/>.
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
20
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
21 */
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
22
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
23 #ifdef HAVE_CONFIG_H
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
24 #include <config.h>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
25 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
26
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
27 #include <vector>
5794
1138ced03f14 [project @ 2006-05-08 20:23:04 by jwe]
jwe
parents: 5792
diff changeset
28
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
29 #include "MArray2.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
30 #include "MSparse.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
31 #include "SparseQR.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
32 #include "SparseCmplxQR.h"
5785
6b9cec830d72 [project @ 2006-05-03 19:32:46 by dbateman]
dbateman
parents: 5775
diff changeset
33 #include "MatrixType.h"
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
34 #include "oct-sort.h"
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
35
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
36 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
37 static MSparse<T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
38 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
39 const octave_idx_type *Q, octave_idx_type rst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
40 octave_idx_type rend, octave_idx_type cst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
41 octave_idx_type cend, octave_idx_type maxnz = -1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
42 bool lazy = false)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
43 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
44 octave_idx_type nz = (rend - rst) * (cend - cst);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
45 maxnz = (maxnz < 0 ? A.nnz () : maxnz);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
46 MSparse<T> B (rend - rst, cend - cst, (nz < maxnz ? nz : maxnz));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
47 // Some sparse functions can support lazy indexing (where elements
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
48 // in the row are in no particular order), even though octave in
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
49 // general can't. For those functions that can using it is a big
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
50 // win here in terms of speed.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
51 if (lazy)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
52 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
53 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
54 for (octave_idx_type j = cst ; j < cend ; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
55 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
56 octave_idx_type qq = (Q ? Q [j] : j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
57 B.xcidx (j - cst) = nz;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
58 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
59 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
60 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
61 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
62 if (r >= rst && r < rend)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
63 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
64 B.xdata (nz) = A.data (p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
65 B.xridx (nz++) = r - rst ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
66 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
67 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
68 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
69 B.xcidx (cend - cst) = nz ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
70 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
71 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
72 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
73 OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
74 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
75 octave_idx_type *ri = B.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
76 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
77 for (octave_idx_type j = cst ; j < cend ; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
78 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
79 octave_idx_type qq = (Q ? Q [j] : j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
80 B.xcidx (j - cst) = nz;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
81 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
82 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
83 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
84 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
85 if (r >= rst && r < rend)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
86 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
87 X [r-rst] = A.data (p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
88 B.xridx (nz++) = r - rst ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
89 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
90 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
91 sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
92 for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
93 B.xdata (p) = X [B.xridx (p)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
94 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
95 B.xcidx (cend - cst) = nz ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
96 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
97
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
98 return B;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
99 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
100
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
101 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
102 static MSparse<double>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
103 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
104 const octave_idx_type *Q, octave_idx_type rst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
105 octave_idx_type rend, octave_idx_type cst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
106 octave_idx_type cend, octave_idx_type maxnz,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
107 bool lazy);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
108
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
109 static MSparse<Complex>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
110 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
111 const octave_idx_type *Q, octave_idx_type rst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
112 octave_idx_type rend, octave_idx_type cst,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
113 octave_idx_type cend, octave_idx_type maxnz,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
114 bool lazy);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
115 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
116
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
117 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
118 static MArray2<T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
119 dmsolve_extract (const MArray2<T> &m, const octave_idx_type *,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
120 const octave_idx_type *, octave_idx_type r1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
121 octave_idx_type r2, octave_idx_type c1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
122 octave_idx_type c2)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
123 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
124 r2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
125 c2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
126 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
127 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
128
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
129 octave_idx_type new_r = r2 - r1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
130 octave_idx_type new_c = c2 - c1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
131
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
132 MArray2<T> result (new_r, new_c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
133
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
134 for (octave_idx_type j = 0; j < new_c; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
135 for (octave_idx_type i = 0; i < new_r; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
136 result.xelem (i, j) = m.elem (r1+i, c1+j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
137
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
138 return result;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
139 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
140
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
141 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
142 static MArray2<double>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
143 dmsolve_extract (const MArray2<double> &m, const octave_idx_type *,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
144 const octave_idx_type *, octave_idx_type r1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
145 octave_idx_type r2, octave_idx_type c1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
146 octave_idx_type c2)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
147
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
148 static MArray2<Complex>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
149 dmsolve_extract (const MArray2<Complex> &m, const octave_idx_type *,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
150 const octave_idx_type *, octave_idx_type r1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
151 octave_idx_type r2, octave_idx_type c1,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
152 octave_idx_type c2)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
153 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
154
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
155 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
156 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
157 dmsolve_insert (MArray2<T> &a, const MArray2<T> &b, const octave_idx_type *Q,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
158 octave_idx_type r, octave_idx_type c)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
159 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
160 T *ax = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
161 const T *bx = b.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
162 octave_idx_type anr = a.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
163 octave_idx_type nr = b.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
164 octave_idx_type nc = b.cols();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
165 for (octave_idx_type j = 0; j < nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
166 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
167 octave_idx_type aoff = (c + j) * anr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
168 octave_idx_type boff = j * nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
169 for (octave_idx_type i = 0; i < nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
170 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
171 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
172 ax [Q [r + i] + aoff] = bx [i + boff];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
173 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
174 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
175 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
176
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
177 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
178 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
179 dmsolve_insert (MArray2<double> &a, const MArray2<double> &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
180 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
181
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
182 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
183 dmsolve_insert (MArray2<Complex> &a, const MArray2<Complex> &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
184 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
185 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
186
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
187 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
188 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
189 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
190 octave_idx_type r, octave_idx_type c)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
191 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
192 octave_idx_type b_rows = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
193 octave_idx_type b_cols = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
194 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
195 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
196
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
197 OCTAVE_LOCAL_BUFFER (octave_idx_type, Qinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
198 for (octave_idx_type i = 0; i < nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
199 Qinv [Q [i]] = i;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
200
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
201 // First count the number of elements in the final array
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
202 octave_idx_type nel = a.xcidx(c) + b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
203
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
204 if (c + b_cols < nc)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
205 nel += a.xcidx(nc) - a.xcidx(c + b_cols);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
206
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
207 for (octave_idx_type i = c; i < c + b_cols; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
208 for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
209 if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
210 nel++;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
211
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
212 OCTAVE_LOCAL_BUFFER (T, X, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
213 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
214 MSparse<T> tmp (a);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
215 a = MSparse<T> (nr, nc, nel);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
216 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
217
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
218 for (octave_idx_type i = 0; i < tmp.cidx(c); i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
219 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
220 a.xdata(i) = tmp.xdata(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
221 a.xridx(i) = tmp.xridx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
222 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
223 for (octave_idx_type i = 0; i < c + 1; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
224 a.xcidx(i) = tmp.xcidx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
225
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
226 octave_idx_type ii = a.xcidx(c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
227
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
228 for (octave_idx_type i = c; i < c + b_cols; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
229 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
230 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
231
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
232 for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
233 if (Qinv [tmp.xridx(j)] < r || Qinv [tmp.xridx(j)] >= r + b_rows)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
234 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
235 X [tmp.xridx(j)] = tmp.xdata(j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
236 a.xridx(ii++) = tmp.xridx(j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
237 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
238
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
239 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
240
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
241 for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
242 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
243 X [Q [r + b.ridx(j)]] = b.data(j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
244 a.xridx(ii++) = Q [r + b.ridx(j)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
245 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
246
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
247 sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
248 for (octave_idx_type p = a.xcidx (i); p < ii; p++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
249 a.xdata (p) = X [a.xridx (p)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
250 a.xcidx(i+1) = ii;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
251 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
252
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
253 for (octave_idx_type i = c + b_cols; i < nc; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
254 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
255 for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
256 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
257 a.xdata(ii) = tmp.xdata(j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
258 a.xridx(ii++) = tmp.xridx(j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
259 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
260 a.xcidx(i+1) = ii;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
261 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
262 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
263
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
264 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
265 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
266 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
267 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
268
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
269 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
270 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
271 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
272 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
273
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
274 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
275 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
276 dmsolve_permute (MArray2<RT> &a, const MArray2<T>& b, const octave_idx_type *p)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
277 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
278 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
279 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
280 const T *Bx = b.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
281 a.resize(b_nr, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
282 RT *Btx = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
283 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
284 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
285 octave_idx_type off = j * b_nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
286 for (octave_idx_type i = 0; i < b_nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
287 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
288 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
289 Btx [p [i] + off] = Bx [ i + off];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
290 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
291 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
292 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
293
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
294 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
295 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
296 dmsolve_permute (MArray2<double> &a, const MArray2<double>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
297 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
298
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
299 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
300 dmsolve_permute (MArray2<Complex> &a, const MArray2<double>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
301 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
302
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
303 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
304 dmsolve_permute (MArray2<Complex> &a, const MArray2<Complex>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
305 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
306 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
307
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
308 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
309 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
310 dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
311 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
312 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
313 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
314 octave_idx_type b_nz = b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
315 octave_idx_type nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
316 a = MSparse<RT> (b_nr, b_nc, b_nz);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
317 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
318 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
319 OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
320 a.xcidx(0) = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
321 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
322 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
323 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
324 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
325 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
326 octave_idx_type r = p [b.ridx (i)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
327 X [r] = b.data (i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
328 a.xridx(nz++) = p [b.ridx (i)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
329 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
330 sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
331 for (octave_idx_type i = a.cidx (j); i < nz; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
332 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
333 OCTAVE_QUIT;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
334 a.xdata (i) = X [a.xridx (i)];
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
335 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
336 a.xcidx(j+1) = nz;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
337 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
338 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
339
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
340 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
341 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
342 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
343 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
344
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
345 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
346 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
347 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
348
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
349 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
350 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
351 const octave_idx_type *p);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
352 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
353
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
354 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
355 solve_singularity_warning (double)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
356 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
357 // Dummy singularity handler so that LU solver doesn't flag
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
358 // an error for numerically rank defficient matrices
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
359 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
360
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
361 template <class RT, class ST, class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
362 RT
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
363 dmsolve (const ST &a, const T &b, octave_idx_type &info)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
364 {
5684
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
365 #ifdef HAVE_CXSPARSE
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
366 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
367 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
368 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
369 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
370 RT retval;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
371
6924
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman]
dbateman
parents: 5797
diff changeset
372 if (nr < 0 || nc < 0 || nr != b_nr)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
373 (*current_liboctave_error_handler)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
374 ("matrix dimension mismatch in solution of minimum norm problem");
6924
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman]
dbateman
parents: 5797
diff changeset
375 else if (nr == 0 || nc == 0 || b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman]
dbateman
parents: 5797
diff changeset
376 retval = RT (nc, b_nc, 0.0);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
377 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
378 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
379 octave_idx_type nnz_remaining = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
380 CXSPARSE_DNAME () csm;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
381 csm.m = nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
382 csm.n = nc;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
383 csm.x = NULL;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
384 csm.nz = -1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
385 csm.nzmax = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
386 // Cast away const on A, with full knowledge that CSparse won't touch it.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
387 // Prevents the methods below making a copy of the data.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
388 csm.p = const_cast<octave_idx_type *>(a.cidx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
389 csm.i = const_cast<octave_idx_type *>(a.ridx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
390
5792
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
391 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
392 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm, 0);
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
393 octave_idx_type *p = dm->p;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
394 octave_idx_type *q = dm->q;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
395 #else
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
396 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
397 octave_idx_type *p = dm->P;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
398 octave_idx_type *q = dm->Q;
5792
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
399 #endif
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
400 OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
401 for (octave_idx_type i = 0; i < nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
402 pinv [p [i]] = i;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
403 RT btmp;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
404 dmsolve_permute (btmp, b, pinv);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
405 info = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
406 retval.resize (nc, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
407
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
408 // Leading over-determined block
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
409 if (dm->rr [2] < nr && dm->cc [3] < nc)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
410 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
411 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
412 nnz_remaining, true);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
413 nnz_remaining -= m.nnz();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
414 RT mtmp =
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
415 qrsolve (m, dmsolve_extract (btmp, NULL, NULL, dm->rr[2], b_nr, 0,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
416 b_nc), info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
417 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
418 if (dm->rr [2] > 0 && !info)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
419 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
420 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
421 dm->cc [3], nc, nnz_remaining, true);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
422 nnz_remaining -= m.nnz();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
423 RT ctmp = dmsolve_extract (btmp, NULL, NULL, 0,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
424 dm->rr[2], 0, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
425 btmp.insert (ctmp - m * mtmp, 0, 0);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
426 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
427 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
428
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
429 // Structurally non-singular blocks
5775
ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe]
jwe
parents: 5697
diff changeset
430 // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
431 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
432 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
433 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
434 dm->cc [2], dm->cc [3], nnz_remaining, false);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
435 nnz_remaining -= m.nnz();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
436 RT btmp2 = dmsolve_extract (btmp, NULL, NULL, dm->rr [1], dm->rr [2],
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
437 0, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
438 double rcond = 0.0;
5785
6b9cec830d72 [project @ 2006-05-03 19:32:46 by dbateman]
dbateman
parents: 5775
diff changeset
439 MatrixType mtyp (MatrixType::Full);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
440 RT mtmp = m.solve (mtyp, btmp2, info, rcond,
5697
2fe20065a545 [project @ 2006-03-21 16:01:46 by dbateman]
dbateman
parents: 5684
diff changeset
441 solve_singularity_warning, false);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
442 if (info != 0)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
443 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
444 info = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
445 mtmp = qrsolve (m, btmp2, info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
446 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
447
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
448 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
449 if (dm->rr [1] > 0 && !info)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
450 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
451 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
452 dm->cc [3], nnz_remaining, true);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
453 nnz_remaining -= m.nnz();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
454 RT ctmp = dmsolve_extract (btmp, NULL, NULL, 0,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
455 dm->rr[1], 0, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
456 btmp.insert (ctmp - m * mtmp, 0, 0);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
457 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
458 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
459
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
460 // Trailing under-determined block
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
461 if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
462 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
463 ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
464 dm->cc [2], nnz_remaining, true);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
465 RT mtmp =
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
466 qrsolve (m, dmsolve_extract(btmp, NULL, NULL, 0, dm->rr [1] , 0,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
467 b_nc), info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
468 dmsolve_insert (retval, mtmp, q, 0, 0);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
469 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
470
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
471 CXSPARSE_DNAME (_dfree) (dm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
472 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
473 return retval;
5684
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
474 #else
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
475 return RT ();
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
476 #endif
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
477 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
478
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
479 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
480 extern Matrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
481 dmsolve (const SparseMatrix &a, const Matrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
482 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
483
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
484 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
485 dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
486 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
487
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
488 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
489 dmsolve (const SparseComplexMatrix &a, const Matrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
490 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
491
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
492 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
493 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
494 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
495
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
496 extern SparseMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
497 dmsolve (const SparseMatrix &a, const SparseMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
498 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
499
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
500 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
501 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
502 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
503
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
504 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
505 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
506 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
507
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
508 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
509 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
510 octave_idx_type &info);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
511 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
512
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
513 /*
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
514 ;;; Local Variables: ***
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
515 ;;; mode: C++ ***
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
516 ;;; End: ***
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
517 */