Mercurial > hg > octave-lyh

comparison liboctave/sparse-dmsolve.cc @ 10314:07ebe522dac2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
untabify liboctave C++ sources
author John W. Eaton <jwe@octave.org>
date Thu, 11 Feb 2010 12:23:32 -0500
parents 4c0cdbe0acca
children 12884915a8e4
comparison
equal deleted inserted replaced
10313:f3b65e1ae355 10314:07ebe522dac2
35 #include "oct-locbuf.h" 35 #include "oct-locbuf.h"
36 36
37 template <class T> 37 template <class T>
38 static MSparse<T> 38 static MSparse<T>
39 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv, 39 dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
40 const octave_idx_type *Q, octave_idx_type rst, 40 const octave_idx_type *Q, octave_idx_type rst,
41 octave_idx_type rend, octave_idx_type cst, 41 octave_idx_type rend, octave_idx_type cst,
42 octave_idx_type cend, octave_idx_type maxnz = -1, 42 octave_idx_type cend, octave_idx_type maxnz = -1,
43 bool lazy = false) 43 bool lazy = false)
44 { 44 {
45 octave_idx_type nz = (rend - rst) * (cend - cst); 45 octave_idx_type nz = (rend - rst) * (cend - cst);
46 maxnz = (maxnz < 0 ? A.nnz () : maxnz); 46 maxnz = (maxnz < 0 ? A.nnz () : maxnz);
47 MSparse<T> B (rend - rst, cend - cst, (nz < maxnz ? nz : maxnz)); 47 MSparse<T> B (rend - rst, cend - cst, (nz < maxnz ? nz : maxnz));
48 // Some sparse functions can support lazy indexing (where elements 48 // Some sparse functions can support lazy indexing (where elements
51 // win here in terms of speed. 51 // win here in terms of speed.
52 if (lazy) 52 if (lazy)
53 { 53 {
54 nz = 0; 54 nz = 0;
55 for (octave_idx_type j = cst ; j < cend ; j++) 55 for (octave_idx_type j = cst ; j < cend ; j++)
56 { 56 {
57 octave_idx_type qq = (Q ? Q [j] : j); 57 octave_idx_type qq = (Q ? Q [j] : j);
58 B.xcidx (j - cst) = nz; 58 B.xcidx (j - cst) = nz;
59 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++) 59 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
60 { 60 {
61 octave_quit (); 61 octave_quit ();
62 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p)); 62 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
63 if (r >= rst && r < rend) 63 if (r >= rst && r < rend)
64 { 64 {
65 B.xdata (nz) = A.data (p); 65 B.xdata (nz) = A.data (p);
66 B.xridx (nz++) = r - rst ; 66 B.xridx (nz++) = r - rst ;
67 } 67 }
68 } 68 }
69 } 69 }
70 B.xcidx (cend - cst) = nz ; 70 B.xcidx (cend - cst) = nz ;
71 } 71 }
72 else 72 else
73 { 73 {
74 OCTAVE_LOCAL_BUFFER (T, X, rend - rst); 74 OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
75 octave_sort<octave_idx_type> sort; 75 octave_sort<octave_idx_type> sort;
76 octave_idx_type *ri = B.xridx(); 76 octave_idx_type *ri = B.xridx();
77 nz = 0; 77 nz = 0;
78 for (octave_idx_type j = cst ; j < cend ; j++) 78 for (octave_idx_type j = cst ; j < cend ; j++)
79 { 79 {
80 octave_idx_type qq = (Q ? Q [j] : j); 80 octave_idx_type qq = (Q ? Q [j] : j);
81 B.xcidx (j - cst) = nz; 81 B.xcidx (j - cst) = nz;
82 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++) 82 for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
83 { 83 {
84 octave_quit (); 84 octave_quit ();
85 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p)); 85 octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
86 if (r >= rst && r < rend) 86 if (r >= rst && r < rend)
87 { 87 {
88 X [r-rst] = A.data (p); 88 X [r-rst] = A.data (p);
89 B.xridx (nz++) = r - rst ; 89 B.xridx (nz++) = r - rst ;
90 } 90 }
91 } 91 }
92 sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst)); 92 sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
93 for (octave_idx_type p = B.cidx (j - cst); p < nz; p++) 93 for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
94 B.xdata (p) = X [B.xridx (p)]; 94 B.xdata (p) = X [B.xridx (p)];
95 } 95 }
96 B.xcidx (cend - cst) = nz ; 96 B.xcidx (cend - cst) = nz ;
97 } 97 }
98 98
99 return B; 99 return B;
100 } 100 }
101 101
102 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 102 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
103 static MSparse<double> 103 static MSparse<double>
104 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv, 104 dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
105 const octave_idx_type *Q, octave_idx_type rst, 105 const octave_idx_type *Q, octave_idx_type rst,
106 octave_idx_type rend, octave_idx_type cst, 106 octave_idx_type rend, octave_idx_type cst,
107 octave_idx_type cend, octave_idx_type maxnz, 107 octave_idx_type cend, octave_idx_type maxnz,
108 bool lazy); 108 bool lazy);
109 109
110 static MSparse<Complex> 110 static MSparse<Complex>
111 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv, 111 dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
112 const octave_idx_type *Q, octave_idx_type rst, 112 const octave_idx_type *Q, octave_idx_type rst,
113 octave_idx_type rend, octave_idx_type cst, 113 octave_idx_type rend, octave_idx_type cst,
114 octave_idx_type cend, octave_idx_type maxnz, 114 octave_idx_type cend, octave_idx_type maxnz,
115 bool lazy); 115 bool lazy);
116 #endif 116 #endif
117 117
118 template <class T> 118 template <class T>
119 static MArray2<T> 119 static MArray2<T>
120 dmsolve_extract (const MArray2<T> &m, const octave_idx_type *, 120 dmsolve_extract (const MArray2<T> &m, const octave_idx_type *,
121 const octave_idx_type *, octave_idx_type r1, 121 const octave_idx_type *, octave_idx_type r1,
122 octave_idx_type r2, octave_idx_type c1, 122 octave_idx_type r2, octave_idx_type c1,
123 octave_idx_type c2) 123 octave_idx_type c2)
124 { 124 {
125 r2 -= 1; 125 r2 -= 1;
126 c2 -= 1; 126 c2 -= 1;
127 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; } 127 if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
128 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; } 128 if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
140 } 140 }
141 141
142 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 142 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
143 static MArray2<double> 143 static MArray2<double>
144 dmsolve_extract (const MArray2<double> &m, const octave_idx_type *, 144 dmsolve_extract (const MArray2<double> &m, const octave_idx_type *,
145 const octave_idx_type *, octave_idx_type r1, 145 const octave_idx_type *, octave_idx_type r1,
146 octave_idx_type r2, octave_idx_type c1, 146 octave_idx_type r2, octave_idx_type c1,
147 octave_idx_type c2) 147 octave_idx_type c2)
148 148
149 static MArray2<Complex> 149 static MArray2<Complex>
150 dmsolve_extract (const MArray2<Complex> &m, const octave_idx_type *, 150 dmsolve_extract (const MArray2<Complex> &m, const octave_idx_type *,
151 const octave_idx_type *, octave_idx_type r1, 151 const octave_idx_type *, octave_idx_type r1,
152 octave_idx_type r2, octave_idx_type c1, 152 octave_idx_type r2, octave_idx_type c1,
153 octave_idx_type c2) 153 octave_idx_type c2)
154 #endif 154 #endif
155 155
156 template <class T> 156 template <class T>
157 static void 157 static void
158 dmsolve_insert (MArray2<T> &a, const MArray2<T> &b, const octave_idx_type *Q, 158 dmsolve_insert (MArray2<T> &a, const MArray2<T> &b, const octave_idx_type *Q,
159 octave_idx_type r, octave_idx_type c) 159 octave_idx_type r, octave_idx_type c)
160 { 160 {
161 T *ax = a.fortran_vec(); 161 T *ax = a.fortran_vec();
162 const T *bx = b.fortran_vec(); 162 const T *bx = b.fortran_vec();
163 octave_idx_type anr = a.rows(); 163 octave_idx_type anr = a.rows();
164 octave_idx_type nr = b.rows(); 164 octave_idx_type nr = b.rows();
166 for (octave_idx_type j = 0; j < nc; j++) 166 for (octave_idx_type j = 0; j < nc; j++)
167 { 167 {
168 octave_idx_type aoff = (c + j) * anr; 168 octave_idx_type aoff = (c + j) * anr;
169 octave_idx_type boff = j * nr; 169 octave_idx_type boff = j * nr;
170 for (octave_idx_type i = 0; i < nr; i++) 170 for (octave_idx_type i = 0; i < nr; i++)
171 { 171 {
172 octave_quit (); 172 octave_quit ();
173 ax [Q [r + i] + aoff] = bx [i + boff]; 173 ax [Q [r + i] + aoff] = bx [i + boff];
174 } 174 }
175 } 175 }
176 } 176 }
177 177
178 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 178 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
179 static void 179 static void
180 dmsolve_insert (MArray2<double> &a, const MArray2<double> &b, 180 dmsolve_insert (MArray2<double> &a, const MArray2<double> &b,
181 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c); 181 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
182 182
183 static void 183 static void
184 dmsolve_insert (MArray2<Complex> &a, const MArray2<Complex> &b, 184 dmsolve_insert (MArray2<Complex> &a, const MArray2<Complex> &b,
185 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c); 185 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
186 #endif 186 #endif
187 187
188 template <class T> 188 template <class T>
189 static void 189 static void
190 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q, 190 dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
191 octave_idx_type r, octave_idx_type c) 191 octave_idx_type r, octave_idx_type c)
192 { 192 {
193 octave_idx_type b_rows = b.rows (); 193 octave_idx_type b_rows = b.rows ();
194 octave_idx_type b_cols = b.cols (); 194 octave_idx_type b_cols = b.cols ();
195 octave_idx_type nr = a.rows (); 195 octave_idx_type nr = a.rows ();
196 octave_idx_type nc = a.cols (); 196 octave_idx_type nc = a.cols ();
206 nel += a.xcidx(nc) - a.xcidx(c + b_cols); 206 nel += a.xcidx(nc) - a.xcidx(c + b_cols);
207 207
208 for (octave_idx_type i = c; i < c + b_cols; i++) 208 for (octave_idx_type i = c; i < c + b_cols; i++)
209 for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++) 209 for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++)
210 if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows) 210 if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows)
211 nel++; 211 nel++;
212 212
213 OCTAVE_LOCAL_BUFFER (T, X, nr); 213 OCTAVE_LOCAL_BUFFER (T, X, nr);
214 octave_sort<octave_idx_type> sort; 214 octave_sort<octave_idx_type> sort;
215 MSparse<T> tmp (a); 215 MSparse<T> tmp (a);
216 a = MSparse<T> (nr, nc, nel); 216 a = MSparse<T> (nr, nc, nel);
229 for (octave_idx_type i = c; i < c + b_cols; i++) 229 for (octave_idx_type i = c; i < c + b_cols; i++)
230 { 230 {
231 octave_quit (); 231 octave_quit ();
232 232
233 for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++) 233 for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
234 if (Qinv [tmp.xridx(j)] < r || Qinv [tmp.xridx(j)] >= r + b_rows) 234 if (Qinv [tmp.xridx(j)] < r || Qinv [tmp.xridx(j)] >= r + b_rows)
235 { 235 {
236 X [tmp.xridx(j)] = tmp.xdata(j); 236 X [tmp.xridx(j)] = tmp.xdata(j);
237 a.xridx(ii++) = tmp.xridx(j); 237 a.xridx(ii++) = tmp.xridx(j);
238 } 238 }
239 239
240 octave_quit (); 240 octave_quit ();
241 241
242 for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++) 242 for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
243 { 243 {
244 X [Q [r + b.ridx(j)]] = b.data(j); 244 X [Q [r + b.ridx(j)]] = b.data(j);
245 a.xridx(ii++) = Q [r + b.ridx(j)]; 245 a.xridx(ii++) = Q [r + b.ridx(j)];
246 } 246 }
247 247
248 sort.sort (ri + a.xcidx (i), ii - a.xcidx (i)); 248 sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
249 for (octave_idx_type p = a.xcidx (i); p < ii; p++) 249 for (octave_idx_type p = a.xcidx (i); p < ii; p++)
250 a.xdata (p) = X [a.xridx (p)]; 250 a.xdata (p) = X [a.xridx (p)];
251 a.xcidx(i+1) = ii; 251 a.xcidx(i+1) = ii;
252 } 252 }
253 253
254 for (octave_idx_type i = c + b_cols; i < nc; i++) 254 for (octave_idx_type i = c + b_cols; i < nc; i++)
255 { 255 {
256 for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++) 256 for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
257 { 257 {
258 a.xdata(ii) = tmp.xdata(j); 258 a.xdata(ii) = tmp.xdata(j);
259 a.xridx(ii++) = tmp.xridx(j); 259 a.xridx(ii++) = tmp.xridx(j);
260 } 260 }
261 a.xcidx(i+1) = ii; 261 a.xcidx(i+1) = ii;
262 } 262 }
263 } 263 }
264 264
265 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 265 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
266 static void 266 static void
267 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b, 267 dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
268 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c); 268 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
269 269
270 static void 270 static void
271 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b, 271 dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
272 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c); 272 const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
273 #endif 273 #endif
274 274
275 template <class T, class RT> 275 template <class T, class RT>
276 static void 276 static void
277 dmsolve_permute (MArray2<RT> &a, const MArray2<T>& b, const octave_idx_type *p) 277 dmsolve_permute (MArray2<RT> &a, const MArray2<T>& b, const octave_idx_type *p)
283 RT *Btx = a.fortran_vec(); 283 RT *Btx = a.fortran_vec();
284 for (octave_idx_type j = 0; j < b_nc; j++) 284 for (octave_idx_type j = 0; j < b_nc; j++)
285 { 285 {
286 octave_idx_type off = j * b_nr; 286 octave_idx_type off = j * b_nr;
287 for (octave_idx_type i = 0; i < b_nr; i++) 287 for (octave_idx_type i = 0; i < b_nr; i++)
288 { 288 {
289 octave_quit (); 289 octave_quit ();
290 Btx [p [i] + off] = Bx [ i + off]; 290 Btx [p [i] + off] = Bx [ i + off];
291 } 291 }
292 } 292 }
293 } 293 }
294 294
295 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 295 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
296 static void 296 static void
297 dmsolve_permute (MArray2<double> &a, const MArray2<double>& b, 297 dmsolve_permute (MArray2<double> &a, const MArray2<double>& b,
298 const octave_idx_type *p); 298 const octave_idx_type *p);
299 299
300 static void 300 static void
301 dmsolve_permute (MArray2<Complex> &a, const MArray2<double>& b, 301 dmsolve_permute (MArray2<Complex> &a, const MArray2<double>& b,
302 const octave_idx_type *p); 302 const octave_idx_type *p);
303 303
304 static void 304 static void
305 dmsolve_permute (MArray2<Complex> &a, const MArray2<Complex>& b, 305 dmsolve_permute (MArray2<Complex> &a, const MArray2<Complex>& b,
306 const octave_idx_type *p); 306 const octave_idx_type *p);
307 #endif 307 #endif
308 308
309 template <class T, class RT> 309 template <class T, class RT>
310 static void 310 static void
311 dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p) 311 dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p)
320 OCTAVE_LOCAL_BUFFER (RT, X, b_nr); 320 OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
321 a.xcidx(0) = 0; 321 a.xcidx(0) = 0;
322 for (octave_idx_type j = 0; j < b_nc; j++) 322 for (octave_idx_type j = 0; j < b_nc; j++)
323 { 323 {
324 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) 324 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
325 { 325 {
326 octave_quit (); 326 octave_quit ();
327 octave_idx_type r = p [b.ridx (i)]; 327 octave_idx_type r = p [b.ridx (i)];
328 X [r] = b.data (i); 328 X [r] = b.data (i);
329 a.xridx(nz++) = p [b.ridx (i)]; 329 a.xridx(nz++) = p [b.ridx (i)];
330 } 330 }
331 sort.sort (ri + a.xcidx (j), nz - a.xcidx (j)); 331 sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
332 for (octave_idx_type i = a.cidx (j); i < nz; i++) 332 for (octave_idx_type i = a.cidx (j); i < nz; i++)
333 { 333 {
334 octave_quit (); 334 octave_quit ();
335 a.xdata (i) = X [a.xridx (i)]; 335 a.xdata (i) = X [a.xridx (i)];
336 } 336 }
337 a.xcidx(j+1) = nz; 337 a.xcidx(j+1) = nz;
338 } 338 }
339 } 339 }
340 340
341 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 341 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
342 static void 342 static void
343 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b, 343 dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
344 const octave_idx_type *p); 344 const octave_idx_type *p);
345 345
346 static void 346 static void
347 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b, 347 dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
348 const octave_idx_type *p); 348 const octave_idx_type *p);
349 349
350 static void 350 static void
351 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b, 351 dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
352 const octave_idx_type *p); 352 const octave_idx_type *p);
353 #endif 353 #endif
354 354
355 static void 355 static void
356 solve_singularity_warning (double) 356 solve_singularity_warning (double)
357 { 357 {
398 octave_idx_type *p = dm->P; 398 octave_idx_type *p = dm->P;
399 octave_idx_type *q = dm->Q; 399 octave_idx_type *q = dm->Q;
400 #endif 400 #endif
401 OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr); 401 OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
402 for (octave_idx_type i = 0; i < nr; i++) 402 for (octave_idx_type i = 0; i < nr; i++)
403 pinv [p [i]] = i; 403 pinv [p [i]] = i;
404 RT btmp; 404 RT btmp;
405 dmsolve_permute (btmp, b, pinv); 405 dmsolve_permute (btmp, b, pinv);
406 info = 0; 406 info = 0;
407 retval.resize (nc, b_nc); 407 retval.resize (nc, b_nc);
408 408
409 // Leading over-determined block 409 // Leading over-determined block
410 if (dm->rr [2] < nr && dm->cc [3] < nc) 410 if (dm->rr [2] < nr && dm->cc [3] < nc)
411 { 411 {
412 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc, 412 ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
413 nnz_remaining, true); 413 nnz_remaining, true);
414 nnz_remaining -= m.nnz(); 414 nnz_remaining -= m.nnz();
415 RT mtmp = 415 RT mtmp =
416 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0, 416 qrsolve (m, dmsolve_extract (btmp, 0, 0, dm->rr[2], b_nr, 0,
417 b_nc), info); 417 b_nc), info);
418 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0); 418 dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
419 if (dm->rr [2] > 0 && !info) 419 if (dm->rr [2] > 0 && !info)
420 { 420 {
421 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2], 421 m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
422 dm->cc [3], nc, nnz_remaining, true); 422 dm->cc [3], nc, nnz_remaining, true);
423 nnz_remaining -= m.nnz(); 423 nnz_remaining -= m.nnz();
424 RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 424 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
425 dm->rr[2], 0, b_nc); 425 dm->rr[2], 0, b_nc);
426 btmp.insert (ctmp - m * mtmp, 0, 0); 426 btmp.insert (ctmp - m * mtmp, 0, 0);
427 } 427 }
428 } 428 }
429 429
430 // Structurally non-singular blocks 430 // Structurally non-singular blocks
431 // FIXME Should use fine Dulmange-Mendelsohn decomposition here. 431 // FIXME Should use fine Dulmange-Mendelsohn decomposition here.
432 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info) 432 if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] && !info)
433 { 433 {
434 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2], 434 ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
435 dm->cc [2], dm->cc [3], nnz_remaining, false); 435 dm->cc [2], dm->cc [3], nnz_remaining, false);
436 nnz_remaining -= m.nnz(); 436 nnz_remaining -= m.nnz();
437 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2], 437 RT btmp2 = dmsolve_extract (btmp, 0, 0, dm->rr [1], dm->rr [2],
438 0, b_nc); 438 0, b_nc);
439 double rcond = 0.0; 439 double rcond = 0.0;
440 MatrixType mtyp (MatrixType::Full); 440 MatrixType mtyp (MatrixType::Full);
441 RT mtmp = m.solve (mtyp, btmp2, info, rcond, 441 RT mtmp = m.solve (mtyp, btmp2, info, rcond,
442 solve_singularity_warning, false); 442 solve_singularity_warning, false);
443 if (info != 0) 443 if (info != 0)
444 { 444 {
445 info = 0; 445 info = 0;
446 mtmp = qrsolve (m, btmp2, info); 446 mtmp = qrsolve (m, btmp2, info);
447 } 447 }
448 448
449 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0); 449 dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
450 if (dm->rr [1] > 0 && !info) 450 if (dm->rr [1] > 0 && !info)
451 { 451 {
452 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2], 452 m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
453 dm->cc [3], nnz_remaining, true); 453 dm->cc [3], nnz_remaining, true);
454 nnz_remaining -= m.nnz(); 454 nnz_remaining -= m.nnz();
455 RT ctmp = dmsolve_extract (btmp, 0, 0, 0, 455 RT ctmp = dmsolve_extract (btmp, 0, 0, 0,
456 dm->rr[1], 0, b_nc); 456 dm->rr[1], 0, b_nc);
457 btmp.insert (ctmp - m * mtmp, 0, 0); 457 btmp.insert (ctmp - m * mtmp, 0, 0);
458 } 458 }
459 } 459 }
460 460
461 // Trailing under-determined block 461 // Trailing under-determined block
462 if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info) 462 if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info)
463 { 463 {
464 ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0, 464 ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0,
465 dm->cc [2], nnz_remaining, true); 465 dm->cc [2], nnz_remaining, true);
466 RT mtmp = 466 RT mtmp =
467 qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0, 467 qrsolve (m, dmsolve_extract(btmp, 0, 0, 0, dm->rr [1] , 0,
468 b_nc), info); 468 b_nc), info);
469 dmsolve_insert (retval, mtmp, q, 0, 0); 469 dmsolve_insert (retval, mtmp, q, 0, 0);
470 } 470 }
471 471
472 CXSPARSE_DNAME (_dfree) (dm); 472 CXSPARSE_DNAME (_dfree) (dm);
473 } 473 }
474 return retval; 474 return retval;
475 #else 475 #else
478 } 478 }
479 479
480 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) 480 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
481 extern Matrix 481 extern Matrix
482 dmsolve (const SparseMatrix &a, const Matrix &b, 482 dmsolve (const SparseMatrix &a, const Matrix &b,
483 octave_idx_type &info); 483 octave_idx_type &info);
484 484
485 extern ComplexMatrix 485 extern ComplexMatrix
486 dmsolve (const SparseMatrix &a, const ComplexMatrix &b, 486 dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
487 octave_idx_type &info); 487 octave_idx_type &info);
488 488
489 extern ComplexMatrix 489 extern ComplexMatrix
490 dmsolve (const SparseComplexMatrix &a, const Matrix &b, 490 dmsolve (const SparseComplexMatrix &a, const Matrix &b,
491 octave_idx_type &info); 491 octave_idx_type &info);
492 492
493 extern ComplexMatrix 493 extern ComplexMatrix
494 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b, 494 dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
495 octave_idx_type &info); 495 octave_idx_type &info);
496 496
497 extern SparseMatrix 497 extern SparseMatrix
498 dmsolve (const SparseMatrix &a, const SparseMatrix &b, 498 dmsolve (const SparseMatrix &a, const SparseMatrix &b,
499 octave_idx_type &info); 499 octave_idx_type &info);
500 500
501 extern SparseComplexMatrix 501 extern SparseComplexMatrix
502 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b, 502 dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
503 octave_idx_type &info); 503 octave_idx_type &info);
504 504
505 extern SparseComplexMatrix 505 extern SparseComplexMatrix
506 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b, 506 dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
507 octave_idx_type &info); 507 octave_idx_type &info);
508 508
509 extern SparseComplexMatrix 509 extern SparseComplexMatrix
510 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b, 510 dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
511 octave_idx_type &info); 511 octave_idx_type &info);
512 #endif 512 #endif