Mercurial > hg > octave-nkf

annotate liboctave/sparse-dmsolve.cc @ 11574:a83bad07f7e3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
attempt better backward compatibility for Array resize functions
author John W. Eaton <jwe@octave.org>
date Thu, 20 Jan 2011 00:12:05 -0500
parents 57632dea2446
children 12df7854fa7c
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
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 10350
diff changeset
3 Copyright (C) 2006-2011 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
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
29 #include "MArray.h"
5683
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"
8377
25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents: 7520
diff changeset
35 #include "oct-locbuf.h"
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
36
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
37 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
38 static MSparse<T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
39 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
40 const octave_idx_type *Q, octave_idx_type rst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
41 octave_idx_type rend, octave_idx_type cst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
42 octave_idx_type cend, octave_idx_type maxnz = -1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
43 bool lazy = false)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
44 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
45 octave_idx_type nz = (rend - rst) * (cend - cst);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
46 maxnz = (maxnz < 0 ? A.nnz () : maxnz);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
47 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
48 // Some sparse functions can support lazy indexing (where elements
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
49 // 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
50 // 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
51 // win here in terms of speed.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
52 if (lazy)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
53 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
54 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
55 for (octave_idx_type j = cst ; j < cend ; j++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
56 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
57 octave_idx_type qq = (Q ? Q [j] : j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
58 B.xcidx (j - cst) = nz;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
59 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
60 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
61 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
62 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
63 if (r >= rst && r < rend)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
64 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
65 B.xdata (nz) = A.data (p);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
66 B.xridx (nz++) = r - rst ;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
67 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
68 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
69 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
70 B.xcidx (cend - cst) = nz ;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
71 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
72 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
73 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
74 OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
75 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
76 octave_idx_type *ri = B.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
77 nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
78 for (octave_idx_type j = cst ; j < cend ; j++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
79 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
80 octave_idx_type qq = (Q ? Q [j] : j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
81 B.xcidx (j - cst) = nz;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
82 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
83 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
84 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
85 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
86 if (r >= rst && r < rend)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
87 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
88 X [r-rst] = A.data (p);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
89 B.xridx (nz++) = r - rst ;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
90 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
91 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
92 sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
93 for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
94 B.xdata (p) = X [B.xridx (p)];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
95 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
96 B.xcidx (cend - cst) = nz ;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
99 return B;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
102 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
103 static MSparse<double>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
104 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
105 const octave_idx_type *Q, octave_idx_type rst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
106 octave_idx_type rend, octave_idx_type cst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
107 octave_idx_type cend, octave_idx_type maxnz,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
108 bool lazy);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
109
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
110 static MSparse<Complex>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
111 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
112 const octave_idx_type *Q, octave_idx_type rst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
113 octave_idx_type rend, octave_idx_type cst,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
114 octave_idx_type cend, octave_idx_type maxnz,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
115 bool lazy);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
116 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
117
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
118 template <class T>
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
119 static MArray<T>
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
120 dmsolve_extract (const MArray<T> &m, const octave_idx_type *,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
121 const octave_idx_type *, octave_idx_type r1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
122 octave_idx_type r2, octave_idx_type c1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
123 octave_idx_type c2)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
124 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
125 r2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
126 c2 -= 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
127 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
128 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
129
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
130 octave_idx_type new_r = r2 - r1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
131 octave_idx_type new_c = c2 - c1 + 1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
132
11570
57632dea2446 attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
133 MArray<T> result (dim_vector (new_r, new_c));
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
134
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
135 for (octave_idx_type j = 0; j < new_c; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
136 for (octave_idx_type i = 0; i < new_r; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
137 result.xelem (i, j) = m.elem (r1+i, c1+j);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
138
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
139 return result;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
142 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
143 static MArray<double>
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
144 dmsolve_extract (const MArray<double> &m, const octave_idx_type *,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
145 const octave_idx_type *, octave_idx_type r1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
146 octave_idx_type r2, octave_idx_type c1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
147 octave_idx_type c2)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
148
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
149 static MArray<Complex>
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
150 dmsolve_extract (const MArray<Complex> &m, const octave_idx_type *,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
151 const octave_idx_type *, octave_idx_type r1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
152 octave_idx_type r2, octave_idx_type c1,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
153 octave_idx_type c2)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
154 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
155
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
156 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
157 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
158 dmsolve_insert (MArray<T> &a, const MArray<T> &b, const octave_idx_type *Q,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
159 octave_idx_type r, octave_idx_type c)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
160 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
161 T *ax = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
162 const T *bx = b.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
163 octave_idx_type anr = a.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
164 octave_idx_type nr = b.rows();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
165 octave_idx_type nc = b.cols();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
166 for (octave_idx_type j = 0; j < nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
167 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
168 octave_idx_type aoff = (c + j) * anr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
169 octave_idx_type boff = j * nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
170 for (octave_idx_type i = 0; i < nr; i++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
171 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
172 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
173 ax [Q [r + i] + aoff] = bx [i + boff];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
174 }
5683
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
178 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
179 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
180 dmsolve_insert (MArray<double> &a, const MArray<double> &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
181 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
182
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
183 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
184 dmsolve_insert (MArray<Complex> &a, const MArray<Complex> &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
185 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
186 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
187
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
188 template <class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
189 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
190 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
191 octave_idx_type r, octave_idx_type c)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
192 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
193 octave_idx_type b_rows = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
194 octave_idx_type b_cols = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
195 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
196 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
197
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
198 OCTAVE_LOCAL_BUFFER (octave_idx_type, Qinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
199 for (octave_idx_type i = 0; i < nr; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
200 Qinv [Q [i]] = i;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
201
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
202 // First count the number of elements in the final array
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
203 octave_idx_type nel = a.xcidx(c) + b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
204
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
205 if (c + b_cols < nc)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
206 nel += a.xcidx(nc) - a.xcidx(c + b_cols);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
207
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
208 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
209 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
210 if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
211 nel++;
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
212
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
213 OCTAVE_LOCAL_BUFFER (T, X, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
214 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
215 MSparse<T> tmp (a);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
216 a = MSparse<T> (nr, nc, nel);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
217 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
218
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
219 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
220 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
221 a.xdata(i) = tmp.xdata(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
222 a.xridx(i) = tmp.xridx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
223 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
224 for (octave_idx_type i = 0; i < c + 1; i++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
225 a.xcidx(i) = tmp.xcidx(i);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
226
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
227 octave_idx_type ii = a.xcidx(c);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
228
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
229 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
230 {
10142
829e69ec3110 make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents: 8920
diff changeset
231 octave_quit ();
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
232
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
233 for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
234 if (Qinv [tmp.xridx(j)] < r || Qinv [tmp.xridx(j)] >= r + b_rows)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
235 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
236 X [tmp.xridx(j)] = tmp.xdata(j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
237 a.xridx(ii++) = tmp.xridx(j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
238 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
239
10142
829e69ec3110 make OCTAVE_QUIT a function
Jaroslav Hajek <highegg@gmail.com>
parents: 8920
diff changeset
240 octave_quit ();
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
241
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
242 for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
243 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
244 X [Q [r + b.ridx(j)]] = b.data(j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
245 a.xridx(ii++) = Q [r + b.ridx(j)];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
246 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
247
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
248 sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
249 for (octave_idx_type p = a.xcidx (i); p < ii; p++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
250 a.xdata (p) = X [a.xridx (p)];
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
251 a.xcidx(i+1) = ii;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
254 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
255 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
256 for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
257 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
258 a.xdata(ii) = tmp.xdata(j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
259 a.xridx(ii++) = tmp.xridx(j);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
260 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
261 a.xcidx(i+1) = ii;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
265 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
266 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
267 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
268 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
269
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
270 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
271 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
272 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
273 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
274
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
275 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
276 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
277 dmsolve_permute (MArray<RT> &a, const MArray<T>& b, const octave_idx_type *p)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
278 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
279 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
280 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
281 const T *Bx = b.fortran_vec();
11574
a83bad07f7e3 attempt better backward compatibility for Array resize functions
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
282 a.resize (dim_vector (b_nr, b_nc));
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
283 RT *Btx = a.fortran_vec();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
284 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
285 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
286 octave_idx_type off = j * b_nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
287 for (octave_idx_type i = 0; i < b_nr; i++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
288 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
289 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
290 Btx [p [i] + off] = Bx [ i + off];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
291 }
5683
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
295 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
296 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
297 dmsolve_permute (MArray<double> &a, const MArray<double>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
298 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
299
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
300 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
301 dmsolve_permute (MArray<Complex> &a, const MArray<double>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
302 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
303
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
304 static void
10350
12884915a8e4 merge MArray classes & improve Array interface
Jaroslav Hajek <highegg@gmail.com>
parents: 10314
diff changeset
305 dmsolve_permute (MArray<Complex> &a, const MArray<Complex>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
306 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
307 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
308
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
309 template <class T, class RT>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
310 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
311 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
312 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
313 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
314 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
315 octave_idx_type b_nz = b.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
316 octave_idx_type nz = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
317 a = MSparse<RT> (b_nr, b_nc, b_nz);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
318 octave_sort<octave_idx_type> sort;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
319 octave_idx_type *ri = a.xridx();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
320 OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
321 a.xcidx(0) = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
322 for (octave_idx_type j = 0; j < b_nc; j++)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
323 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
324 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
325 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
326 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
327 octave_idx_type r = p [b.ridx (i)];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
328 X [r] = b.data (i);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
329 a.xridx(nz++) = p [b.ridx (i)];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
330 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
331 sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
332 for (octave_idx_type i = a.cidx (j); i < nz; i++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
333 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
334 octave_quit ();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
335 a.xdata (i) = X [a.xridx (i)];
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
336 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
337 a.xcidx(j+1) = nz;
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
341 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
342 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
343 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
344 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
345
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
346 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
347 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
348 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
349
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
350 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
351 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
352 const octave_idx_type *p);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
353 #endif
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
354
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
355 static void
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
356 solve_singularity_warning (double)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
357 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
358 // Dummy singularity handler so that LU solver doesn't flag
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
359 // an error for numerically rank defficient matrices
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
362 template <class RT, class ST, class T>
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
363 RT
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
364 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
365 {
5684
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
366 #ifdef HAVE_CXSPARSE
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
367 octave_idx_type nr = a.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
368 octave_idx_type nc = a.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
369 octave_idx_type b_nr = b.rows ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
370 octave_idx_type b_nc = b.cols ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
371 RT retval;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
372
6924
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman]
dbateman
parents: 5797
diff changeset
373 if (nr < 0 || nc < 0 || nr != b_nr)
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
374 (*current_liboctave_error_handler)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
375 ("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
376 else if (nr == 0 || nc == 0 || b_nc == 0)
be176b7e110a [project @ 2007-09-25 19:08:19 by dbateman]
dbateman
parents: 5797
diff changeset
377 retval = RT (nc, b_nc, 0.0);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
378 else
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
379 {
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
380 octave_idx_type nnz_remaining = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
381 CXSPARSE_DNAME () csm;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
382 csm.m = nr;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
383 csm.n = nc;
7520
b166043585a8 use 0 instead of NULL
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
384 csm.x = 0;
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
385 csm.nz = -1;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
386 csm.nzmax = a.nnz ();
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
387 // 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
388 // Prevents the methods below making a copy of the data.
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
389 csm.p = const_cast<octave_idx_type *>(a.cidx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
390 csm.i = const_cast<octave_idx_type *>(a.ridx ());
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
391
5792
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
392 #if defined(CS_VER) && (CS_VER >= 2)
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
393 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm, 0);
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
394 octave_idx_type *p = dm->p;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
395 octave_idx_type *q = dm->q;
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
396 #else
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
397 CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
398 octave_idx_type *p = dm->P;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
399 octave_idx_type *q = dm->Q;
5792
eb90c83b4f91 [project @ 2006-05-04 20:14:49 by dbateman]
dbateman
parents: 5785
diff changeset
400 #endif
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
401 OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
402 for (octave_idx_type i = 0; i < nr; i++)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
403 pinv [p [i]] = i;
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
404 RT btmp;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
405 dmsolve_permute (btmp, b, pinv);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
406 info = 0;
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
407 retval.resize (nc, b_nc);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
408
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
409 // Leading over-determined block
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
410 if (dm->rr [2] < nr && dm->cc [3] < nc)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
411 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
412 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
413 nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
414 nnz_remaining -= m.nnz();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
415 RT mtmp =
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
416 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
417 b_nc), info);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
418 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
419 if (dm->rr [2] > 0 && !info)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
420 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
421 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
422 dm->cc [3], nc, nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
423 nnz_remaining -= m.nnz();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
424 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
425 dm->rr[2], 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
426 btmp.insert (ctmp - m * mtmp, 0, 0);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
427 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
428 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
429
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
430 // Structurally non-singular blocks
5775
ace8d8d26933 [project @ 2006-04-24 19:13:06 by jwe]
jwe
parents: 5697
diff changeset
431 // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
432 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
433 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
434 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
435 dm->cc [2], dm->cc [3], nnz_remaining, false);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
436 nnz_remaining -= m.nnz();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
437 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2],
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
438 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
439 double rcond = 0.0;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
440 MatrixType mtyp (MatrixType::Full);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
441 RT mtmp = m.solve (mtyp, btmp2, info, rcond,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
442 solve_singularity_warning, false);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
443 if (info != 0)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
444 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
445 info = 0;
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
446 mtmp = qrsolve (m, btmp2, info);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
447 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
448
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
449 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
450 if (dm->rr [1] > 0 && !info)
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
451 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
452 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
453 dm->cc [3], nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
454 nnz_remaining -= m.nnz();
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
455 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
456 dm->rr[1], 0, b_nc);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
457 btmp.insert (ctmp - m * mtmp, 0, 0);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
458 }
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
459 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
460
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
461 // Trailing under-determined block
5797
11fcab4c461d [project @ 2006-05-09 06:15:17 by dbateman]
dbateman
parents: 5794
diff changeset
462 if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info)
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
463 {
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
464 ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
465 dm->cc [2], nnz_remaining, true);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
466 RT mtmp =
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
467 qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0,
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
468 b_nc), info);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
469 dmsolve_insert (retval, mtmp, q, 0, 0);
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
470 }
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
471
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
472 CXSPARSE_DNAME (_dfree) (dm);
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
473 }
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
474 return retval;
5684
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
475 #else
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
476 return RT ();
f39adfd0913d [project @ 2006-03-16 18:12:01 by dbateman]
dbateman
parents: 5683
diff changeset
477 #endif
5683
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
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
480 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
481 extern Matrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
482 dmsolve (const SparseMatrix &a, const Matrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
483 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
484
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
485 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
486 dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
487 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
488
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
489 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
490 dmsolve (const SparseComplexMatrix &a, const Matrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
491 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
492
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
493 extern ComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
494 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
495 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
496
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
497 extern SparseMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
498 dmsolve (const SparseMatrix &a, const SparseMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
499 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
500
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
501 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
502 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
503 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
504
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
505 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
506 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
507 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
508
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
509 extern SparseComplexMatrix
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
510 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
10314
07ebe522dac2 untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents: 10158
diff changeset
511 octave_idx_type &info);
5683
a6bfe12c4c72 [project @ 2006-03-16 18:02:27 by dbateman]
dbateman
parents:
diff changeset
512 #endif