5164
|
1 /* |
|
2 |
|
3 Copyright (C) 2004 David Bateman |
|
4 Copyright (C) 1998-2004 Andy Adler |
|
5 |
|
6 Octave is free software; you can redistribute it and/or modify it |
|
7 under the terms of the GNU General Public License as published by the |
|
8 Free Software Foundation; either version 2, or (at your option) any |
|
9 later version. |
|
10 |
|
11 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
12 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
13 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
14 for more details. |
|
15 |
|
16 You should have received a copy of the GNU General Public License |
|
17 along with this program; see the file COPYING. If not, write to the Free |
|
18 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. |
|
19 |
|
20 */ |
|
21 |
|
22 #ifdef HAVE_CONFIG_H |
|
23 #include <config.h> |
|
24 #endif |
|
25 |
|
26 #include <cfloat> |
|
27 |
|
28 #include <iostream> |
|
29 #include <vector> |
|
30 |
|
31 #include "quit.h" |
|
32 #include "lo-ieee.h" |
|
33 #include "lo-mappers.h" |
|
34 #include "f77-fcn.h" |
|
35 #include "dRowVector.h" |
|
36 |
|
37 #include "CSparse.h" |
|
38 #include "boolSparse.h" |
|
39 #include "dSparse.h" |
|
40 #include "oct-spparms.h" |
|
41 #include "SparsedbleLU.h" |
|
42 #include "SparseType.h" |
|
43 |
|
44 // External UMFPACK functions in C |
|
45 extern "C" { |
|
46 #include "umfpack.h" |
|
47 } |
|
48 |
|
49 // Fortran functions we call. |
|
50 extern "C" |
|
51 { |
|
52 F77_RET_T |
|
53 F77_FUNC (dgbtrf, DGBTRF) (const int&, const int&, const int&, |
|
54 const int&, double*, const int&, int*, int&); |
|
55 |
|
56 F77_RET_T |
|
57 F77_FUNC (dgbtrs, DGBTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
58 const int&, const int&, const int&, |
|
59 const double*, const int&, |
|
60 const int*, double*, const int&, int& |
|
61 F77_CHAR_ARG_LEN_DECL); |
|
62 |
|
63 F77_RET_T |
|
64 F77_FUNC (dgbcon, DGBCON) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
65 const int&, const int&, double*, |
|
66 const int&, const int*, const double&, |
|
67 double&, double*, int*, int& |
|
68 F77_CHAR_ARG_LEN_DECL); |
|
69 |
|
70 F77_RET_T |
|
71 F77_FUNC (dpbtrf, DPBTRF) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
72 const int&, double*, const int&, int& |
|
73 F77_CHAR_ARG_LEN_DECL); |
|
74 |
|
75 F77_RET_T |
|
76 F77_FUNC (dpbtrs, DPBTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
77 const int&, const int&, double*, const int&, |
|
78 double*, const int&, int& |
|
79 F77_CHAR_ARG_LEN_DECL); |
|
80 |
|
81 F77_RET_T |
|
82 F77_FUNC (dpbcon, DPBCON) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
83 const int&, double*, const int&, |
|
84 const double&, double&, double*, int*, int& |
|
85 F77_CHAR_ARG_LEN_DECL); |
|
86 F77_RET_T |
|
87 F77_FUNC (dptsv, DPTSV) (const int&, const int&, double*, double*, |
|
88 double*, const int&, int&); |
|
89 |
|
90 F77_RET_T |
|
91 F77_FUNC (dgtsv, DGTSV) (const int&, const int&, double*, double*, |
|
92 double*, double*, const int&, int&); |
|
93 |
|
94 F77_RET_T |
|
95 F77_FUNC (dgttrf, DGTTRF) (const int&, double*, double*, double*, double*, |
|
96 int*, int&); |
|
97 |
|
98 F77_RET_T |
|
99 F77_FUNC (dgttrs, DGTTRS) (F77_CONST_CHAR_ARG_DECL, const int&, |
|
100 const int&, const double*, const double*, |
|
101 const double*, const double*, const int*, |
|
102 double *, const int&, int& |
|
103 F77_CHAR_ARG_LEN_DECL); |
|
104 |
|
105 F77_RET_T |
|
106 F77_FUNC (zptsv, ZPTSV) (const int&, const int&, Complex*, Complex*, |
|
107 Complex*, const int&, int&); |
|
108 |
|
109 F77_RET_T |
|
110 F77_FUNC (zgtsv, ZGTSV) (const int&, const int&, Complex*, Complex*, |
|
111 Complex*, Complex*, const int&, int&); |
|
112 |
|
113 } |
|
114 |
|
115 SparseMatrix::SparseMatrix (const SparseBoolMatrix &a) |
|
116 : MSparse<double> (a.rows (), a.cols (), a.nnz ()) |
|
117 { |
|
118 int nc = cols (); |
|
119 int nz = nnz (); |
|
120 |
|
121 for (int i = 0; i < nc + 1; i++) |
|
122 cidx (i) = a.cidx (i); |
|
123 |
|
124 for (int i = 0; i < nz; i++) |
|
125 { |
|
126 data (i) = a.data (i); |
|
127 ridx (i) = a.ridx (i); |
|
128 } |
|
129 } |
|
130 |
|
131 bool |
|
132 SparseMatrix::operator == (const SparseMatrix& a) const |
|
133 { |
|
134 int nr = rows (); |
|
135 int nc = cols (); |
|
136 int nz = nnz (); |
|
137 int nr_a = a.rows (); |
|
138 int nc_a = a.cols (); |
|
139 int nz_a = a.nnz (); |
|
140 |
|
141 if (nr != nr_a || nc != nc_a || nz != nz_a) |
|
142 return false; |
|
143 |
|
144 for (int i = 0; i < nc + 1; i++) |
|
145 if (cidx(i) != a.cidx(i)) |
|
146 return false; |
|
147 |
|
148 for (int i = 0; i < nz; i++) |
|
149 if (data(i) != a.data(i) || ridx(i) != a.ridx(i)) |
|
150 return false; |
|
151 |
|
152 return true; |
|
153 } |
|
154 |
|
155 bool |
|
156 SparseMatrix::operator != (const SparseMatrix& a) const |
|
157 { |
|
158 return !(*this == a); |
|
159 } |
|
160 |
|
161 bool |
|
162 SparseMatrix::is_symmetric (void) const |
|
163 { |
|
164 if (is_square () && rows () > 0) |
|
165 { |
|
166 for (int i = 0; i < rows (); i++) |
|
167 for (int j = i+1; j < cols (); j++) |
|
168 if (elem (i, j) != elem (j, i)) |
|
169 return false; |
|
170 |
|
171 return true; |
|
172 } |
|
173 |
|
174 return false; |
|
175 } |
|
176 |
|
177 SparseMatrix& |
|
178 SparseMatrix::insert (const SparseMatrix& a, int r, int c) |
|
179 { |
|
180 MSparse<double>::insert (a, r, c); |
|
181 return *this; |
|
182 } |
|
183 |
|
184 SparseMatrix |
|
185 SparseMatrix::max (int dim) const |
|
186 { |
|
187 Array2<int> dummy_idx; |
|
188 return max (dummy_idx, dim); |
|
189 } |
|
190 |
|
191 SparseMatrix |
|
192 SparseMatrix::max (Array2<int>& idx_arg, int dim) const |
|
193 { |
|
194 SparseMatrix result; |
|
195 dim_vector dv = dims (); |
|
196 |
|
197 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
198 return result; |
|
199 |
|
200 int nr = dv(0); |
|
201 int nc = dv(1); |
|
202 |
|
203 if (dim == 0) |
|
204 { |
|
205 idx_arg.resize (1, nc); |
|
206 int nel = 0; |
|
207 for (int j = 0; j < nc; j++) |
|
208 { |
|
209 double tmp_max = octave_NaN; |
|
210 int idx_j = 0; |
|
211 for (int i = cidx(j); i < cidx(j+1); i++) |
|
212 { |
|
213 if (ridx(i) != idx_j) |
|
214 break; |
|
215 else |
|
216 idx_j++; |
|
217 } |
|
218 |
|
219 if (idx_j != nr) |
|
220 tmp_max = 0.; |
|
221 |
|
222 for (int i = cidx(j); i < cidx(j+1); i++) |
|
223 { |
|
224 double tmp = data (i); |
|
225 |
|
226 if (octave_is_NaN_or_NA (tmp)) |
|
227 continue; |
|
228 else if (octave_is_NaN_or_NA (tmp_max) || tmp > tmp_max) |
|
229 { |
|
230 idx_j = ridx (i); |
|
231 tmp_max = tmp; |
|
232 } |
|
233 |
|
234 } |
|
235 |
|
236 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_max) ? 0 : idx_j; |
|
237 if (tmp_max != 0.) |
|
238 nel++; |
|
239 } |
|
240 |
|
241 result = SparseMatrix (1, nc, nel); |
|
242 |
|
243 int ii = 0; |
|
244 result.xcidx (0) = 0; |
|
245 for (int j = 0; j < nc; j++) |
|
246 { |
|
247 double tmp = elem (idx_arg(j), j); |
|
248 if (tmp != 0.) |
|
249 { |
|
250 result.xdata (ii) = tmp; |
|
251 result.xridx (ii++) = 0; |
|
252 } |
|
253 result.xcidx (j+1) = ii; |
|
254 |
|
255 } |
|
256 } |
|
257 else |
|
258 { |
|
259 idx_arg.resize (nr, 1, 0); |
|
260 |
|
261 for (int i = cidx(0); i < cidx(1); i++) |
|
262 idx_arg.elem(ridx(i)) = -1; |
|
263 |
|
264 for (int j = 0; j < nc; j++) |
|
265 for (int i = 0; i < nr; i++) |
|
266 { |
|
267 if (idx_arg.elem(i) != -1) |
|
268 continue; |
|
269 bool found = false; |
|
270 for (int k = cidx(j); k < cidx(j+1); k++) |
|
271 if (ridx(k) == i) |
|
272 { |
|
273 found = true; |
|
274 break; |
|
275 } |
|
276 |
|
277 if (!found) |
|
278 idx_arg.elem(i) = j; |
|
279 |
|
280 } |
|
281 |
|
282 for (int j = 0; j < nc; j++) |
|
283 { |
|
284 for (int i = cidx(j); i < cidx(j+1); i++) |
|
285 { |
|
286 int ir = ridx (i); |
|
287 int ix = idx_arg.elem (ir); |
|
288 double tmp = data (i); |
|
289 |
|
290 if (octave_is_NaN_or_NA (tmp)) |
|
291 continue; |
|
292 else if (ix == -1 || tmp > elem (ir, ix)) |
|
293 idx_arg.elem (ir) = j; |
|
294 } |
|
295 } |
|
296 |
|
297 int nel = 0; |
|
298 for (int j = 0; j < nr; j++) |
|
299 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
|
300 nel++; |
|
301 |
|
302 result = SparseMatrix (nr, 1, nel); |
|
303 |
|
304 int ii = 0; |
|
305 result.xcidx (0) = 0; |
|
306 result.xcidx (1) = nel; |
|
307 for (int j = 0; j < nr; j++) |
|
308 { |
|
309 if (idx_arg(j) == -1) |
|
310 { |
|
311 idx_arg(j) = 0; |
|
312 result.xdata (ii) = octave_NaN; |
|
313 result.xridx (ii++) = j; |
|
314 } |
|
315 else |
|
316 { |
|
317 double tmp = elem (j, idx_arg(j)); |
|
318 if (tmp != 0.) |
|
319 { |
|
320 result.xdata (ii) = tmp; |
|
321 result.xridx (ii++) = j; |
|
322 } |
|
323 } |
|
324 } |
|
325 } |
|
326 |
|
327 return result; |
|
328 } |
|
329 |
|
330 SparseMatrix |
|
331 SparseMatrix::min (int dim) const |
|
332 { |
|
333 Array2<int> dummy_idx; |
|
334 return min (dummy_idx, dim); |
|
335 } |
|
336 |
|
337 SparseMatrix |
|
338 SparseMatrix::min (Array2<int>& idx_arg, int dim) const |
|
339 { |
|
340 SparseMatrix result; |
|
341 dim_vector dv = dims (); |
|
342 |
|
343 if (dv.numel () == 0 || dim > dv.length () || dim < 0) |
|
344 return result; |
|
345 |
|
346 int nr = dv(0); |
|
347 int nc = dv(1); |
|
348 |
|
349 if (dim == 0) |
|
350 { |
|
351 idx_arg.resize (1, nc); |
|
352 int nel = 0; |
|
353 for (int j = 0; j < nc; j++) |
|
354 { |
|
355 double tmp_min = octave_NaN; |
|
356 int idx_j = 0; |
|
357 for (int i = cidx(j); i < cidx(j+1); i++) |
|
358 { |
|
359 if (ridx(i) != idx_j) |
|
360 break; |
|
361 else |
|
362 idx_j++; |
|
363 } |
|
364 |
|
365 if (idx_j != nr) |
|
366 tmp_min = 0.; |
|
367 |
|
368 for (int i = cidx(j); i < cidx(j+1); i++) |
|
369 { |
|
370 double tmp = data (i); |
|
371 |
|
372 if (octave_is_NaN_or_NA (tmp)) |
|
373 continue; |
|
374 else if (octave_is_NaN_or_NA (tmp_min) || tmp < tmp_min) |
|
375 { |
|
376 idx_j = ridx (i); |
|
377 tmp_min = tmp; |
|
378 } |
|
379 |
|
380 } |
|
381 |
|
382 idx_arg.elem (j) = octave_is_NaN_or_NA (tmp_min) ? 0 : idx_j; |
|
383 if (tmp_min != 0.) |
|
384 nel++; |
|
385 } |
|
386 |
|
387 result = SparseMatrix (1, nc, nel); |
|
388 |
|
389 int ii = 0; |
|
390 result.xcidx (0) = 0; |
|
391 for (int j = 0; j < nc; j++) |
|
392 { |
|
393 double tmp = elem (idx_arg(j), j); |
|
394 if (tmp != 0.) |
|
395 { |
|
396 result.xdata (ii) = tmp; |
|
397 result.xridx (ii++) = 0; |
|
398 } |
|
399 result.xcidx (j+1) = ii; |
|
400 |
|
401 } |
|
402 } |
|
403 else |
|
404 { |
|
405 idx_arg.resize (nr, 1, 0); |
|
406 |
|
407 for (int i = cidx(0); i < cidx(1); i++) |
|
408 idx_arg.elem(ridx(i)) = -1; |
|
409 |
|
410 for (int j = 0; j < nc; j++) |
|
411 for (int i = 0; i < nr; i++) |
|
412 { |
|
413 if (idx_arg.elem(i) != -1) |
|
414 continue; |
|
415 bool found = false; |
|
416 for (int k = cidx(j); k < cidx(j+1); k++) |
|
417 if (ridx(k) == i) |
|
418 { |
|
419 found = true; |
|
420 break; |
|
421 } |
|
422 |
|
423 if (!found) |
|
424 idx_arg.elem(i) = j; |
|
425 |
|
426 } |
|
427 |
|
428 for (int j = 0; j < nc; j++) |
|
429 { |
|
430 for (int i = cidx(j); i < cidx(j+1); i++) |
|
431 { |
|
432 int ir = ridx (i); |
|
433 int ix = idx_arg.elem (ir); |
|
434 double tmp = data (i); |
|
435 |
|
436 if (octave_is_NaN_or_NA (tmp)) |
|
437 continue; |
|
438 else if (ix == -1 || tmp < elem (ir, ix)) |
|
439 idx_arg.elem (ir) = j; |
|
440 } |
|
441 } |
|
442 |
|
443 int nel = 0; |
|
444 for (int j = 0; j < nr; j++) |
|
445 if (idx_arg.elem(j) == -1 || elem (j, idx_arg.elem (j)) != 0.) |
|
446 nel++; |
|
447 |
|
448 result = SparseMatrix (nr, 1, nel); |
|
449 |
|
450 int ii = 0; |
|
451 result.xcidx (0) = 0; |
|
452 result.xcidx (1) = nel; |
|
453 for (int j = 0; j < nr; j++) |
|
454 { |
|
455 if (idx_arg(j) == -1) |
|
456 { |
|
457 idx_arg(j) = 0; |
|
458 result.xdata (ii) = octave_NaN; |
|
459 result.xridx (ii++) = j; |
|
460 } |
|
461 else |
|
462 { |
|
463 double tmp = elem (j, idx_arg(j)); |
|
464 if (tmp != 0.) |
|
465 { |
|
466 result.xdata (ii) = tmp; |
|
467 result.xridx (ii++) = j; |
|
468 } |
|
469 } |
|
470 } |
|
471 } |
|
472 |
|
473 return result; |
|
474 } |
|
475 |
|
476 SparseMatrix |
|
477 SparseMatrix::concat (const SparseMatrix& rb, const Array<int>& ra_idx) |
|
478 { |
|
479 // Don't use numel to avoid all possiblity of an overflow |
|
480 if (rb.rows () > 0 && rb.cols () > 0) |
|
481 insert (rb, ra_idx(0), ra_idx(1)); |
|
482 return *this; |
|
483 } |
|
484 |
|
485 SparseComplexMatrix |
|
486 SparseMatrix::concat (const SparseComplexMatrix& rb, const Array<int>& ra_idx) |
|
487 { |
|
488 SparseComplexMatrix retval (*this); |
|
489 if (rb.rows () > 0 && rb.cols () > 0) |
|
490 retval.insert (rb, ra_idx(0), ra_idx(1)); |
|
491 return retval; |
|
492 } |
|
493 |
|
494 SparseMatrix |
|
495 real (const SparseComplexMatrix& a) |
|
496 { |
|
497 int nr = a.rows (); |
|
498 int nc = a.cols (); |
|
499 int nz = a.nnz (); |
|
500 SparseMatrix r (nr, nc, nz); |
|
501 |
|
502 for (int i = 0; i < nc +1; i++) |
|
503 r.cidx(i) = a.cidx(i); |
|
504 |
|
505 for (int i = 0; i < nz; i++) |
|
506 { |
|
507 r.data(i) = real (a.data(i)); |
|
508 r.ridx(i) = a.ridx(i); |
|
509 } |
|
510 |
|
511 return r; |
|
512 } |
|
513 |
|
514 SparseMatrix |
|
515 imag (const SparseComplexMatrix& a) |
|
516 { |
|
517 int nr = a.rows (); |
|
518 int nc = a.cols (); |
|
519 int nz = a.nnz (); |
|
520 SparseMatrix r (nr, nc, nz); |
|
521 |
|
522 for (int i = 0; i < nc +1; i++) |
|
523 r.cidx(i) = a.cidx(i); |
|
524 |
|
525 for (int i = 0; i < nz; i++) |
|
526 { |
|
527 r.data(i) = imag (a.data(i)); |
|
528 r.ridx(i) = a.ridx(i); |
|
529 } |
|
530 |
|
531 return r; |
|
532 } |
|
533 |
|
534 SparseMatrix |
|
535 atan2 (const double& x, const SparseMatrix& y) |
|
536 { |
|
537 int nr = y.rows (); |
|
538 int nc = y.cols (); |
|
539 |
|
540 if (x == 0.) |
|
541 return SparseMatrix (nr, nc); |
|
542 else |
|
543 { |
|
544 // Its going to be basically full, so this is probably the |
|
545 // best way to handle it. |
|
546 Matrix tmp (nr, nc, atan2 (x, 0.)); |
|
547 |
|
548 for (int j = 0; j < nc; j++) |
|
549 for (int i = y.cidx (j); i < y.cidx (j+1); i++) |
|
550 tmp.elem (y.ridx(i), j) = atan2 (x, y.data(i)); |
|
551 |
|
552 return SparseMatrix (tmp); |
|
553 } |
|
554 } |
|
555 |
|
556 SparseMatrix |
|
557 atan2 (const SparseMatrix& x, const double& y) |
|
558 { |
|
559 int nr = x.rows (); |
|
560 int nc = x.cols (); |
|
561 int nz = x.nnz (); |
|
562 |
|
563 SparseMatrix retval (nr, nc, nz); |
|
564 |
|
565 int ii = 0; |
|
566 retval.xcidx(0) = 0; |
|
567 for (int i = 0; i < nc; i++) |
|
568 { |
|
569 for (int j = x.cidx(i); j < x.cidx(i+1); j++) |
|
570 { |
|
571 double tmp = atan2 (x.data(j), y); |
|
572 if (tmp != 0.) |
|
573 { |
|
574 retval.xdata (ii) = tmp; |
|
575 retval.xridx (ii++) = x.ridx (j); |
|
576 } |
|
577 } |
|
578 retval.xcidx (i+1) = ii; |
|
579 } |
|
580 |
|
581 if (ii != nz) |
|
582 { |
|
583 SparseMatrix retval2 (nr, nc, ii); |
|
584 for (int i = 0; i < nc+1; i++) |
|
585 retval2.xcidx (i) = retval.cidx (i); |
|
586 for (int i = 0; i < ii; i++) |
|
587 { |
|
588 retval2.xdata (i) = retval.data (i); |
|
589 retval2.xridx (i) = retval.ridx (i); |
|
590 } |
|
591 return retval2; |
|
592 } |
|
593 else |
|
594 return retval; |
|
595 } |
|
596 |
|
597 SparseMatrix |
|
598 atan2 (const SparseMatrix& x, const SparseMatrix& y) |
|
599 { |
|
600 SparseMatrix r; |
|
601 |
|
602 if ((x.rows() == y.rows()) && (x.cols() == y.cols())) |
|
603 { |
|
604 int x_nr = x.rows (); |
|
605 int x_nc = x.cols (); |
|
606 |
|
607 int y_nr = y.rows (); |
|
608 int y_nc = y.cols (); |
|
609 |
|
610 if (x_nr != y_nr || x_nc != y_nc) |
|
611 gripe_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc); |
|
612 else |
|
613 { |
|
614 r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ())); |
|
615 |
|
616 int jx = 0; |
|
617 r.cidx (0) = 0; |
|
618 for (int i = 0 ; i < x_nc ; i++) |
|
619 { |
|
620 int ja = x.cidx(i); |
|
621 int ja_max = x.cidx(i+1); |
|
622 bool ja_lt_max= ja < ja_max; |
|
623 |
|
624 int jb = y.cidx(i); |
|
625 int jb_max = y.cidx(i+1); |
|
626 bool jb_lt_max = jb < jb_max; |
|
627 |
|
628 while (ja_lt_max || jb_lt_max ) |
|
629 { |
|
630 OCTAVE_QUIT; |
|
631 if ((! jb_lt_max) || |
|
632 (ja_lt_max && (x.ridx(ja) < y.ridx(jb)))) |
|
633 { |
|
634 r.ridx(jx) = x.ridx(ja); |
|
635 r.data(jx) = atan2 (x.data(ja), 0.); |
|
636 jx++; |
|
637 ja++; |
|
638 ja_lt_max= ja < ja_max; |
|
639 } |
|
640 else if (( !ja_lt_max ) || |
|
641 (jb_lt_max && (y.ridx(jb) < x.ridx(ja)) ) ) |
|
642 { |
|
643 jb++; |
|
644 jb_lt_max= jb < jb_max; |
|
645 } |
|
646 else |
|
647 { |
|
648 double tmp = atan2 (x.data(ja), y.data(jb)); |
|
649 if (tmp != 0.) |
|
650 { |
|
651 r.data(jx) = tmp; |
|
652 r.ridx(jx) = x.ridx(ja); |
|
653 jx++; |
|
654 } |
|
655 ja++; |
|
656 ja_lt_max= ja < ja_max; |
|
657 jb++; |
|
658 jb_lt_max= jb < jb_max; |
|
659 } |
|
660 } |
|
661 r.cidx(i+1) = jx; |
|
662 } |
|
663 |
|
664 r.maybe_compress (); |
|
665 } |
|
666 } |
|
667 else |
|
668 (*current_liboctave_error_handler) ("matrix size mismatch"); |
|
669 |
|
670 return r; |
|
671 } |
|
672 |
|
673 SparseMatrix |
|
674 SparseMatrix::inverse (void) const |
|
675 { |
|
676 int info; |
|
677 double rcond; |
|
678 return inverse (info, rcond, 0, 0); |
|
679 } |
|
680 |
|
681 SparseMatrix |
|
682 SparseMatrix::inverse (int& info) const |
|
683 { |
|
684 double rcond; |
|
685 return inverse (info, rcond, 0, 0); |
|
686 } |
|
687 |
|
688 SparseMatrix |
|
689 SparseMatrix::inverse (int& info, double& rcond, int force, int calc_cond) const |
|
690 { |
|
691 info = -1; |
|
692 (*current_liboctave_error_handler) |
|
693 ("SparseMatrix::inverse not implemented yet"); |
|
694 return SparseMatrix (); |
|
695 } |
|
696 |
|
697 DET |
|
698 SparseMatrix::determinant (void) const |
|
699 { |
|
700 int info; |
|
701 double rcond; |
|
702 return determinant (info, rcond, 0); |
|
703 } |
|
704 |
|
705 DET |
|
706 SparseMatrix::determinant (int& info) const |
|
707 { |
|
708 double rcond; |
|
709 return determinant (info, rcond, 0); |
|
710 } |
|
711 |
|
712 DET |
|
713 SparseMatrix::determinant (int& err, double& rcond, int) const |
|
714 { |
|
715 DET retval; |
|
716 |
|
717 int nr = rows (); |
|
718 int nc = cols (); |
|
719 |
|
720 if (nr == 0 || nc == 0 || nr != nc) |
|
721 { |
|
722 double d[2]; |
|
723 d[0] = 1.0; |
|
724 d[1] = 0.0; |
|
725 retval = DET (d); |
|
726 } |
|
727 else |
|
728 { |
|
729 err = 0; |
|
730 |
|
731 // Setup the control parameters |
|
732 Matrix Control (UMFPACK_CONTROL, 1); |
|
733 double *control = Control.fortran_vec (); |
|
734 umfpack_di_defaults (control); |
|
735 |
|
736 double tmp = Voctave_sparse_controls.get_key ("spumoni"); |
|
737 if (!xisnan (tmp)) |
|
738 Control (UMFPACK_PRL) = tmp; |
|
739 |
|
740 tmp = Voctave_sparse_controls.get_key ("piv_tol"); |
|
741 if (!xisnan (tmp)) |
|
742 { |
|
743 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
|
744 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
745 } |
|
746 |
|
747 // Set whether we are allowed to modify Q or not |
|
748 tmp = Voctave_sparse_controls.get_key ("autoamd"); |
|
749 if (!xisnan (tmp)) |
|
750 Control (UMFPACK_FIXQ) = tmp; |
|
751 |
|
752 // Turn-off UMFPACK scaling for LU |
|
753 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; |
|
754 |
|
755 umfpack_di_report_control (control); |
|
756 |
|
757 const int *Ap = cidx (); |
|
758 const int *Ai = ridx (); |
|
759 const double *Ax = data (); |
|
760 |
|
761 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); |
|
762 |
|
763 void *Symbolic; |
|
764 Matrix Info (1, UMFPACK_INFO); |
|
765 double *info = Info.fortran_vec (); |
|
766 int status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL, |
|
767 &Symbolic, control, info); |
|
768 |
|
769 if (status < 0) |
|
770 { |
|
771 (*current_liboctave_error_handler) |
|
772 ("SparseMatrix::determinant symbolic factorization failed"); |
|
773 |
|
774 umfpack_di_report_status (control, status); |
|
775 umfpack_di_report_info (control, info); |
|
776 |
|
777 umfpack_di_free_symbolic (&Symbolic) ; |
|
778 } |
|
779 else |
|
780 { |
|
781 umfpack_di_report_symbolic (Symbolic, control); |
|
782 |
|
783 void *Numeric; |
|
784 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, |
|
785 control, info) ; |
|
786 umfpack_di_free_symbolic (&Symbolic) ; |
|
787 |
|
788 rcond = Info (UMFPACK_RCOND); |
|
789 |
|
790 if (status < 0) |
|
791 { |
|
792 (*current_liboctave_error_handler) |
|
793 ("SparseMatrix::determinant numeric factorization failed"); |
|
794 |
|
795 umfpack_di_report_status (control, status); |
|
796 umfpack_di_report_info (control, info); |
|
797 |
|
798 umfpack_di_free_numeric (&Numeric); |
|
799 } |
|
800 else |
|
801 { |
|
802 umfpack_di_report_numeric (Numeric, control); |
|
803 |
|
804 double d[2]; |
|
805 |
|
806 status = umfpack_di_get_determinant (&d[0], &d[1], Numeric, |
|
807 info); |
|
808 |
|
809 if (status < 0) |
|
810 { |
|
811 (*current_liboctave_error_handler) |
|
812 ("SparseMatrix::determinant error calculating determinant"); |
|
813 |
|
814 umfpack_di_report_status (control, status); |
|
815 umfpack_di_report_info (control, info); |
|
816 |
|
817 umfpack_di_free_numeric (&Numeric); |
|
818 } |
|
819 else |
|
820 retval = DET (d); |
|
821 } |
|
822 } |
|
823 } |
|
824 |
|
825 return retval; |
|
826 } |
|
827 |
|
828 Matrix |
|
829 SparseMatrix::dsolve (SparseType &mattype, const Matrix& b, int& err, |
|
830 double& rcond, solve_singularity_handler) const |
|
831 { |
|
832 Matrix retval; |
|
833 |
|
834 int nr = rows (); |
|
835 int nc = cols (); |
|
836 err = 0; |
|
837 |
|
838 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
839 (*current_liboctave_error_handler) |
|
840 ("matrix dimension mismatch solution of linear equations"); |
|
841 else |
|
842 { |
|
843 // Print spparms("spumoni") info if requested |
|
844 int typ = mattype.type (); |
|
845 mattype.info (); |
|
846 |
|
847 if (typ == SparseType::Diagonal || |
|
848 typ == SparseType::Permuted_Diagonal) |
|
849 { |
|
850 retval.resize (b.rows (), b.cols()); |
|
851 if (typ == SparseType::Diagonal) |
|
852 for (int j = 0; j < b.cols(); j++) |
|
853 for (int i = 0; i < nr; i++) |
|
854 retval(i,j) = b(i,j) / data (i); |
|
855 else |
|
856 for (int j = 0; j < b.cols(); j++) |
|
857 for (int i = 0; i < nr; i++) |
|
858 retval(i,j) = b(ridx(i),j) / data (i); |
|
859 |
|
860 double dmax = 0., dmin = octave_Inf; |
|
861 for (int i = 0; i < nr; i++) |
|
862 { |
|
863 double tmp = fabs(data(i)); |
|
864 if (tmp > dmax) |
|
865 dmax = tmp; |
|
866 if (tmp < dmin) |
|
867 dmin = tmp; |
|
868 } |
|
869 rcond = dmin / dmax; |
|
870 } |
|
871 else |
|
872 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
873 } |
|
874 |
|
875 return retval; |
|
876 } |
|
877 |
|
878 SparseMatrix |
|
879 SparseMatrix::dsolve (SparseType &mattype, const SparseMatrix& b, int& err, |
|
880 double& rcond, solve_singularity_handler) const |
|
881 { |
|
882 SparseMatrix retval; |
|
883 |
|
884 int nr = rows (); |
|
885 int nc = cols (); |
|
886 err = 0; |
|
887 |
|
888 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
889 (*current_liboctave_error_handler) |
|
890 ("matrix dimension mismatch solution of linear equations"); |
|
891 else |
|
892 { |
|
893 // Print spparms("spumoni") info if requested |
|
894 int typ = mattype.type (); |
|
895 mattype.info (); |
|
896 |
|
897 if (typ == SparseType::Diagonal || |
|
898 typ == SparseType::Permuted_Diagonal) |
|
899 { |
|
900 int b_nr = b.rows (); |
|
901 int b_nc = b.cols (); |
|
902 int b_nz = b.nnz (); |
|
903 retval = SparseMatrix (b_nr, b_nc, b_nz); |
|
904 |
|
905 retval.xcidx(0) = 0; |
|
906 int ii = 0; |
|
907 if (typ == SparseType::Diagonal) |
|
908 for (int j = 0; j < b.cols(); j++) |
|
909 { |
|
910 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
911 { |
|
912 retval.xridx (ii) = b.ridx(i); |
|
913 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); |
|
914 } |
|
915 retval.xcidx(j+1) = ii; |
|
916 } |
|
917 else |
|
918 for (int j = 0; j < b.cols(); j++) |
|
919 { |
|
920 for (int i = 0; i < nr; i++) |
|
921 { |
|
922 bool found = false; |
|
923 int k; |
|
924 for (k = b.cidx(j); k < b.cidx(j+1); k++) |
|
925 if (ridx(i) == b.ridx(k)) |
|
926 { |
|
927 found = true; |
|
928 break; |
|
929 } |
|
930 if (found) |
|
931 { |
|
932 retval.xridx (ii) = i; |
|
933 retval.xdata (ii++) = b.data(k) / data (i); |
|
934 } |
|
935 } |
|
936 retval.xcidx(j+1) = ii; |
|
937 } |
|
938 |
|
939 double dmax = 0., dmin = octave_Inf; |
|
940 for (int i = 0; i < nr; i++) |
|
941 { |
|
942 double tmp = fabs(data(i)); |
|
943 if (tmp > dmax) |
|
944 dmax = tmp; |
|
945 if (tmp < dmin) |
|
946 dmin = tmp; |
|
947 } |
|
948 rcond = dmin / dmax; |
|
949 } |
|
950 else |
|
951 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
952 } |
|
953 |
|
954 return retval; |
|
955 } |
|
956 |
|
957 ComplexMatrix |
|
958 SparseMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b, int& err, |
|
959 double& rcond, solve_singularity_handler) const |
|
960 { |
|
961 ComplexMatrix retval; |
|
962 |
|
963 int nr = rows (); |
|
964 int nc = cols (); |
|
965 err = 0; |
|
966 |
|
967 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
968 (*current_liboctave_error_handler) |
|
969 ("matrix dimension mismatch solution of linear equations"); |
|
970 else |
|
971 { |
|
972 // Print spparms("spumoni") info if requested |
|
973 int typ = mattype.type (); |
|
974 mattype.info (); |
|
975 |
|
976 if (typ == SparseType::Diagonal || |
|
977 typ == SparseType::Permuted_Diagonal) |
|
978 { |
|
979 retval.resize (b.rows (), b.cols()); |
|
980 if (typ == SparseType::Diagonal) |
|
981 for (int j = 0; j < b.cols(); j++) |
|
982 for (int i = 0; i < nr; i++) |
|
983 retval(i,j) = b(i,j) / data (i); |
|
984 else |
|
985 for (int j = 0; j < b.cols(); j++) |
|
986 for (int i = 0; i < nr; i++) |
|
987 retval(i,j) = b(ridx(i),j) / data (i); |
|
988 |
|
989 double dmax = 0., dmin = octave_Inf; |
|
990 for (int i = 0; i < nr; i++) |
|
991 { |
|
992 double tmp = fabs(data(i)); |
|
993 if (tmp > dmax) |
|
994 dmax = tmp; |
|
995 if (tmp < dmin) |
|
996 dmin = tmp; |
|
997 } |
|
998 rcond = dmin / dmax; |
|
999 } |
|
1000 else |
|
1001 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1002 } |
|
1003 |
|
1004 return retval; |
|
1005 } |
|
1006 |
|
1007 SparseComplexMatrix |
|
1008 SparseMatrix::dsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
1009 int& err, double& rcond, |
|
1010 solve_singularity_handler) const |
|
1011 { |
|
1012 SparseComplexMatrix retval; |
|
1013 |
|
1014 int nr = rows (); |
|
1015 int nc = cols (); |
|
1016 err = 0; |
|
1017 |
|
1018 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1019 (*current_liboctave_error_handler) |
|
1020 ("matrix dimension mismatch solution of linear equations"); |
|
1021 else |
|
1022 { |
|
1023 // Print spparms("spumoni") info if requested |
|
1024 int typ = mattype.type (); |
|
1025 mattype.info (); |
|
1026 |
|
1027 if (typ == SparseType::Diagonal || |
|
1028 typ == SparseType::Permuted_Diagonal) |
|
1029 { |
|
1030 int b_nr = b.rows (); |
|
1031 int b_nc = b.cols (); |
|
1032 int b_nz = b.nnz (); |
|
1033 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
1034 |
|
1035 retval.xcidx(0) = 0; |
|
1036 int ii = 0; |
|
1037 if (typ == SparseType::Diagonal) |
|
1038 for (int j = 0; j < b.cols(); j++) |
|
1039 { |
|
1040 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1041 { |
|
1042 retval.xridx (ii) = b.ridx(i); |
|
1043 retval.xdata (ii++) = b.data(i) / data (b.ridx (i)); |
|
1044 } |
|
1045 retval.xcidx(j+1) = ii; |
|
1046 } |
|
1047 else |
|
1048 for (int j = 0; j < b.cols(); j++) |
|
1049 { |
|
1050 for (int i = 0; i < nr; i++) |
|
1051 { |
|
1052 bool found = false; |
|
1053 int k; |
|
1054 for (k = b.cidx(j); k < b.cidx(j+1); k++) |
|
1055 if (ridx(i) == b.ridx(k)) |
|
1056 { |
|
1057 found = true; |
|
1058 break; |
|
1059 } |
|
1060 if (found) |
|
1061 { |
|
1062 retval.xridx (ii) = i; |
|
1063 retval.xdata (ii++) = b.data(k) / data (i); |
|
1064 } |
|
1065 } |
|
1066 retval.xcidx(j+1) = ii; |
|
1067 } |
|
1068 |
|
1069 double dmax = 0., dmin = octave_Inf; |
|
1070 for (int i = 0; i < nr; i++) |
|
1071 { |
|
1072 double tmp = fabs(data(i)); |
|
1073 if (tmp > dmax) |
|
1074 dmax = tmp; |
|
1075 if (tmp < dmin) |
|
1076 dmin = tmp; |
|
1077 } |
|
1078 rcond = dmin / dmax; |
|
1079 } |
|
1080 else |
|
1081 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1082 } |
|
1083 |
|
1084 return retval; |
|
1085 } |
|
1086 |
|
1087 Matrix |
|
1088 SparseMatrix::utsolve (SparseType &mattype, const Matrix& b, int& err, |
|
1089 double& rcond, |
|
1090 solve_singularity_handler sing_handler) const |
|
1091 { |
|
1092 Matrix retval; |
|
1093 |
|
1094 int nr = rows (); |
|
1095 int nc = cols (); |
|
1096 err = 0; |
|
1097 |
|
1098 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1099 (*current_liboctave_error_handler) |
|
1100 ("matrix dimension mismatch solution of linear equations"); |
|
1101 else |
|
1102 { |
|
1103 // Print spparms("spumoni") info if requested |
|
1104 int typ = mattype.type (); |
|
1105 mattype.info (); |
|
1106 |
|
1107 if (typ == SparseType::Permuted_Upper || |
|
1108 typ == SparseType::Upper) |
|
1109 { |
|
1110 double anorm = 0.; |
|
1111 double ainvnorm = 0.; |
|
1112 int b_cols = b.cols (); |
|
1113 rcond = 0.; |
|
1114 |
|
1115 // Calculate the 1-norm of matrix for rcond calculation |
|
1116 for (int j = 0; j < nr; j++) |
|
1117 { |
|
1118 double atmp = 0.; |
|
1119 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1120 atmp += fabs(data(i)); |
|
1121 if (atmp > anorm) |
|
1122 anorm = atmp; |
|
1123 } |
|
1124 |
|
1125 if (typ == SparseType::Permuted_Upper) |
|
1126 { |
|
1127 retval.resize (b.rows (), b.cols ()); |
|
1128 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
1129 int *p_perm = mattype.triangular_row_perm (); |
|
1130 int *q_perm = mattype.triangular_col_perm (); |
|
1131 |
|
1132 (*current_liboctave_warning_handler) |
|
1133 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1134 |
|
1135 for (int j = 0; j < b_cols; j++) |
|
1136 { |
|
1137 for (int i = 0; i < nr; i++) |
|
1138 work[i] = b(i,j); |
|
1139 |
|
1140 for (int k = nr-1; k >= 0; k--) |
|
1141 { |
|
1142 int iidx = q_perm[k]; |
|
1143 if (work[iidx] != 0.) |
|
1144 { |
|
1145 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1146 { |
|
1147 err = -2; |
|
1148 goto triangular_error; |
|
1149 } |
|
1150 |
|
1151 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1152 work[iidx] = tmp; |
|
1153 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1154 { |
|
1155 int idx2 = q_perm[ridx(i)]; |
|
1156 work[idx2] = |
|
1157 work[idx2] - tmp * data(i); |
|
1158 } |
|
1159 } |
|
1160 } |
|
1161 |
|
1162 for (int i = 0; i < nr; i++) |
|
1163 retval (i, j) = work[p_perm[i]]; |
|
1164 } |
|
1165 |
|
1166 // Calculation of 1-norm of inv(*this) |
|
1167 for (int i = 0; i < nr; i++) |
|
1168 work[i] = 0.; |
|
1169 |
|
1170 for (int j = 0; j < nr; j++) |
|
1171 { |
|
1172 work[q_perm[j]] = 1.; |
|
1173 |
|
1174 for (int k = j; k >= 0; k--) |
|
1175 { |
|
1176 int iidx = q_perm[k]; |
|
1177 |
|
1178 if (work[iidx] != 0.) |
|
1179 { |
|
1180 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1181 work[iidx] = tmp; |
|
1182 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1183 { |
|
1184 int idx2 = q_perm[ridx(i)]; |
|
1185 work[idx2] = work[idx2] - tmp * data(i); |
|
1186 } |
|
1187 } |
|
1188 } |
|
1189 double atmp = 0; |
|
1190 for (int i = 0; i < j+1; i++) |
|
1191 { |
|
1192 atmp += fabs(work[i]); |
|
1193 work[i] = 0.; |
|
1194 } |
|
1195 if (atmp > ainvnorm) |
|
1196 ainvnorm = atmp; |
|
1197 } |
|
1198 } |
|
1199 else |
|
1200 { |
|
1201 retval = b; |
|
1202 double *x_vec = retval.fortran_vec (); |
|
1203 |
|
1204 for (int j = 0; j < b_cols; j++) |
|
1205 { |
|
1206 int offset = j * nr; |
|
1207 for (int k = nr-1; k >= 0; k--) |
|
1208 { |
|
1209 if (x_vec[k+offset] != 0.) |
|
1210 { |
|
1211 if (ridx(cidx(k+1)-1) != k) |
|
1212 { |
|
1213 err = -2; |
|
1214 goto triangular_error; |
|
1215 } |
|
1216 |
|
1217 double tmp = x_vec[k+offset] / |
|
1218 data(cidx(k+1)-1); |
|
1219 x_vec[k+offset] = tmp; |
|
1220 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1221 { |
|
1222 int iidx = ridx(i); |
|
1223 x_vec[iidx+offset] = |
|
1224 x_vec[iidx+offset] - tmp * data(i); |
|
1225 } |
|
1226 } |
|
1227 } |
|
1228 } |
|
1229 |
|
1230 // Calculation of 1-norm of inv(*this) |
|
1231 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
1232 for (int i = 0; i < nr; i++) |
|
1233 work[i] = 0.; |
|
1234 |
|
1235 for (int j = 0; j < nr; j++) |
|
1236 { |
|
1237 work[j] = 1.; |
|
1238 |
|
1239 for (int k = j; k >= 0; k--) |
|
1240 { |
|
1241 if (work[k] != 0.) |
|
1242 { |
|
1243 double tmp = work[k] / data(cidx(k+1)-1); |
|
1244 work[k] = tmp; |
|
1245 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1246 { |
|
1247 int iidx = ridx(i); |
|
1248 work[iidx] = work[iidx] - tmp * data(i); |
|
1249 } |
|
1250 } |
|
1251 } |
|
1252 double atmp = 0; |
|
1253 for (int i = 0; i < j+1; i++) |
|
1254 { |
|
1255 atmp += fabs(work[i]); |
|
1256 work[i] = 0.; |
|
1257 } |
|
1258 if (atmp > ainvnorm) |
|
1259 ainvnorm = atmp; |
|
1260 } |
|
1261 } |
|
1262 |
|
1263 rcond = 1. / ainvnorm / anorm; |
|
1264 |
|
1265 triangular_error: |
|
1266 if (err != 0) |
|
1267 { |
|
1268 if (sing_handler) |
|
1269 sing_handler (rcond); |
|
1270 else |
|
1271 (*current_liboctave_error_handler) |
|
1272 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1273 rcond); |
|
1274 } |
|
1275 |
|
1276 volatile double rcond_plus_one = rcond + 1.0; |
|
1277 |
|
1278 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1279 { |
|
1280 err = -2; |
|
1281 |
|
1282 if (sing_handler) |
|
1283 sing_handler (rcond); |
|
1284 else |
|
1285 (*current_liboctave_error_handler) |
|
1286 ("matrix singular to machine precision, rcond = %g", |
|
1287 rcond); |
|
1288 } |
|
1289 } |
|
1290 else |
|
1291 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1292 } |
|
1293 |
|
1294 return retval; |
|
1295 } |
|
1296 |
|
1297 SparseMatrix |
|
1298 SparseMatrix::utsolve (SparseType &mattype, const SparseMatrix& b, int& err, |
|
1299 double& rcond, solve_singularity_handler sing_handler) const |
|
1300 { |
|
1301 SparseMatrix retval; |
|
1302 |
|
1303 int nr = rows (); |
|
1304 int nc = cols (); |
|
1305 err = 0; |
|
1306 |
|
1307 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1308 (*current_liboctave_error_handler) |
|
1309 ("matrix dimension mismatch solution of linear equations"); |
|
1310 else |
|
1311 { |
|
1312 // Print spparms("spumoni") info if requested |
|
1313 int typ = mattype.type (); |
|
1314 mattype.info (); |
|
1315 |
|
1316 if (typ == SparseType::Permuted_Upper || |
|
1317 typ == SparseType::Upper) |
|
1318 { |
|
1319 double anorm = 0.; |
|
1320 double ainvnorm = 0.; |
|
1321 rcond = 0.; |
|
1322 |
|
1323 // Calculate the 1-norm of matrix for rcond calculation |
|
1324 for (int j = 0; j < nr; j++) |
|
1325 { |
|
1326 double atmp = 0.; |
|
1327 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1328 atmp += fabs(data(i)); |
|
1329 if (atmp > anorm) |
|
1330 anorm = atmp; |
|
1331 } |
|
1332 |
|
1333 int b_nr = b.rows (); |
|
1334 int b_nc = b.cols (); |
|
1335 int b_nz = b.nnz (); |
|
1336 retval = SparseMatrix (b_nr, b_nc, b_nz); |
|
1337 retval.xcidx(0) = 0; |
|
1338 int ii = 0; |
|
1339 int x_nz = b_nz; |
|
1340 |
|
1341 if (typ == SparseType::Permuted_Upper) |
|
1342 { |
|
1343 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
1344 int *p_perm = mattype.triangular_row_perm (); |
|
1345 int *q_perm = mattype.triangular_col_perm (); |
|
1346 |
|
1347 (*current_liboctave_warning_handler) |
|
1348 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1349 |
|
1350 for (int j = 0; j < b_nc; j++) |
|
1351 { |
|
1352 for (int i = 0; i < nr; i++) |
|
1353 work[i] = 0.; |
|
1354 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1355 work[b.ridx(i)] = b.data(i); |
|
1356 |
|
1357 for (int k = nr-1; k >= 0; k--) |
|
1358 { |
|
1359 int iidx = q_perm[k]; |
|
1360 if (work[iidx] != 0.) |
|
1361 { |
|
1362 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1363 { |
|
1364 err = -2; |
|
1365 goto triangular_error; |
|
1366 } |
|
1367 |
|
1368 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1369 work[iidx] = tmp; |
|
1370 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1371 { |
|
1372 int idx2 = q_perm[ridx(i)]; |
|
1373 work[idx2] = |
|
1374 work[idx2] - tmp * data(i); |
|
1375 } |
|
1376 } |
|
1377 } |
|
1378 |
|
1379 // Count non-zeros in work vector and adjust space in |
|
1380 // retval if needed |
|
1381 int new_nnz = 0; |
|
1382 for (int i = 0; i < nr; i++) |
|
1383 if (work[i] != 0.) |
|
1384 new_nnz++; |
|
1385 |
|
1386 if (ii + new_nnz > x_nz) |
|
1387 { |
|
1388 // Resize the sparse matrix |
|
1389 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1390 retval.change_capacity (sz); |
|
1391 x_nz = sz; |
|
1392 } |
|
1393 |
|
1394 for (int i = 0; i < nr; i++) |
|
1395 if (work[p_perm[i]] != 0.) |
|
1396 { |
|
1397 retval.xridx(ii) = i; |
|
1398 retval.xdata(ii++) = work[p_perm[i]]; |
|
1399 } |
|
1400 retval.xcidx(j+1) = ii; |
|
1401 } |
|
1402 |
|
1403 retval.maybe_compress (); |
|
1404 |
|
1405 // Calculation of 1-norm of inv(*this) |
|
1406 for (int i = 0; i < nr; i++) |
|
1407 work[i] = 0.; |
|
1408 |
|
1409 for (int j = 0; j < nr; j++) |
|
1410 { |
|
1411 work[q_perm[j]] = 1.; |
|
1412 |
|
1413 for (int k = j; k >= 0; k--) |
|
1414 { |
|
1415 int iidx = q_perm[k]; |
|
1416 |
|
1417 if (work[iidx] != 0.) |
|
1418 { |
|
1419 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1420 work[iidx] = tmp; |
|
1421 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1422 { |
|
1423 int idx2 = q_perm[ridx(i)]; |
|
1424 work[idx2] = work[idx2] - tmp * data(i); |
|
1425 } |
|
1426 } |
|
1427 } |
|
1428 double atmp = 0; |
|
1429 for (int i = 0; i < j+1; i++) |
|
1430 { |
|
1431 atmp += fabs(work[i]); |
|
1432 work[i] = 0.; |
|
1433 } |
|
1434 if (atmp > ainvnorm) |
|
1435 ainvnorm = atmp; |
|
1436 } |
|
1437 } |
|
1438 else |
|
1439 { |
|
1440 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
1441 |
|
1442 for (int j = 0; j < b_nc; j++) |
|
1443 { |
|
1444 for (int i = 0; i < nr; i++) |
|
1445 work[i] = 0.; |
|
1446 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1447 work[b.ridx(i)] = b.data(i); |
|
1448 |
|
1449 for (int k = nr-1; k >= 0; k--) |
|
1450 { |
|
1451 if (work[k] != 0.) |
|
1452 { |
|
1453 if (ridx(cidx(k+1)-1) != k) |
|
1454 { |
|
1455 err = -2; |
|
1456 goto triangular_error; |
|
1457 } |
|
1458 |
|
1459 double tmp = work[k] / data(cidx(k+1)-1); |
|
1460 work[k] = tmp; |
|
1461 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1462 { |
|
1463 int iidx = ridx(i); |
|
1464 work[iidx] = work[iidx] - tmp * data(i); |
|
1465 } |
|
1466 } |
|
1467 } |
|
1468 |
|
1469 // Count non-zeros in work vector and adjust space in |
|
1470 // retval if needed |
|
1471 int new_nnz = 0; |
|
1472 for (int i = 0; i < nr; i++) |
|
1473 if (work[i] != 0.) |
|
1474 new_nnz++; |
|
1475 |
|
1476 if (ii + new_nnz > x_nz) |
|
1477 { |
|
1478 // Resize the sparse matrix |
|
1479 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1480 retval.change_capacity (sz); |
|
1481 x_nz = sz; |
|
1482 } |
|
1483 |
|
1484 for (int i = 0; i < nr; i++) |
|
1485 if (work[i] != 0.) |
|
1486 { |
|
1487 retval.xridx(ii) = i; |
|
1488 retval.xdata(ii++) = work[i]; |
|
1489 } |
|
1490 retval.xcidx(j+1) = ii; |
|
1491 } |
|
1492 |
|
1493 retval.maybe_compress (); |
|
1494 |
|
1495 // Calculation of 1-norm of inv(*this) |
|
1496 for (int i = 0; i < nr; i++) |
|
1497 work[i] = 0.; |
|
1498 |
|
1499 for (int j = 0; j < nr; j++) |
|
1500 { |
|
1501 work[j] = 1.; |
|
1502 |
|
1503 for (int k = j; k >= 0; k--) |
|
1504 { |
|
1505 if (work[k] != 0.) |
|
1506 { |
|
1507 double tmp = work[k] / data(cidx(k+1)-1); |
|
1508 work[k] = tmp; |
|
1509 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1510 { |
|
1511 int iidx = ridx(i); |
|
1512 work[iidx] = work[iidx] - tmp * data(i); |
|
1513 } |
|
1514 } |
|
1515 } |
|
1516 double atmp = 0; |
|
1517 for (int i = 0; i < j+1; i++) |
|
1518 { |
|
1519 atmp += fabs(work[i]); |
|
1520 work[i] = 0.; |
|
1521 } |
|
1522 if (atmp > ainvnorm) |
|
1523 ainvnorm = atmp; |
|
1524 } |
|
1525 } |
|
1526 |
|
1527 rcond = 1. / ainvnorm / anorm; |
|
1528 |
|
1529 triangular_error: |
|
1530 if (err != 0) |
|
1531 { |
|
1532 if (sing_handler) |
|
1533 sing_handler (rcond); |
|
1534 else |
|
1535 (*current_liboctave_error_handler) |
|
1536 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1537 rcond); |
|
1538 } |
|
1539 |
|
1540 volatile double rcond_plus_one = rcond + 1.0; |
|
1541 |
|
1542 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1543 { |
|
1544 err = -2; |
|
1545 |
|
1546 if (sing_handler) |
|
1547 sing_handler (rcond); |
|
1548 else |
|
1549 (*current_liboctave_error_handler) |
|
1550 ("matrix singular to machine precision, rcond = %g", |
|
1551 rcond); |
|
1552 } |
|
1553 } |
|
1554 else |
|
1555 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1556 } |
|
1557 return retval; |
|
1558 } |
|
1559 |
|
1560 ComplexMatrix |
|
1561 SparseMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b, int& err, |
|
1562 double& rcond, solve_singularity_handler sing_handler) const |
|
1563 { |
|
1564 ComplexMatrix retval; |
|
1565 |
|
1566 int nr = rows (); |
|
1567 int nc = cols (); |
|
1568 err = 0; |
|
1569 |
|
1570 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1571 (*current_liboctave_error_handler) |
|
1572 ("matrix dimension mismatch solution of linear equations"); |
|
1573 else |
|
1574 { |
|
1575 // Print spparms("spumoni") info if requested |
|
1576 int typ = mattype.type (); |
|
1577 mattype.info (); |
|
1578 |
|
1579 if (typ == SparseType::Permuted_Upper || |
|
1580 typ == SparseType::Upper) |
|
1581 { |
|
1582 double anorm = 0.; |
|
1583 double ainvnorm = 0.; |
|
1584 int b_nc = b.cols (); |
|
1585 rcond = 0.; |
|
1586 |
|
1587 // Calculate the 1-norm of matrix for rcond calculation |
|
1588 for (int j = 0; j < nr; j++) |
|
1589 { |
|
1590 double atmp = 0.; |
|
1591 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1592 atmp += fabs(data(i)); |
|
1593 if (atmp > anorm) |
|
1594 anorm = atmp; |
|
1595 } |
|
1596 |
|
1597 if (typ == SparseType::Permuted_Upper) |
|
1598 { |
|
1599 retval.resize (b.rows (), b.cols ()); |
|
1600 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1601 int *p_perm = mattype.triangular_row_perm (); |
|
1602 int *q_perm = mattype.triangular_col_perm (); |
|
1603 |
|
1604 (*current_liboctave_warning_handler) |
|
1605 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1606 |
|
1607 for (int j = 0; j < b_nc; j++) |
|
1608 { |
|
1609 for (int i = 0; i < nr; i++) |
|
1610 work[i] = b(i,j); |
|
1611 |
|
1612 for (int k = nr-1; k >= 0; k--) |
|
1613 { |
|
1614 int iidx = q_perm[k]; |
|
1615 if (work[iidx] != 0.) |
|
1616 { |
|
1617 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1618 { |
|
1619 err = -2; |
|
1620 goto triangular_error; |
|
1621 } |
|
1622 |
|
1623 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1624 work[iidx] = tmp; |
|
1625 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1626 { |
|
1627 int idx2 = q_perm[ridx(i)]; |
|
1628 work[idx2] = |
|
1629 work[idx2] - tmp * data(i); |
|
1630 } |
|
1631 } |
|
1632 } |
|
1633 |
|
1634 for (int i = 0; i < nr; i++) |
|
1635 retval (i, j) = work[p_perm[i]]; |
|
1636 |
|
1637 } |
|
1638 |
|
1639 // Calculation of 1-norm of inv(*this) |
|
1640 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
1641 for (int i = 0; i < nr; i++) |
|
1642 work2[i] = 0.; |
|
1643 |
|
1644 for (int j = 0; j < nr; j++) |
|
1645 { |
|
1646 work2[q_perm[j]] = 1.; |
|
1647 |
|
1648 for (int k = j; k >= 0; k--) |
|
1649 { |
|
1650 int iidx = q_perm[k]; |
|
1651 |
|
1652 if (work2[iidx] != 0.) |
|
1653 { |
|
1654 double tmp = work2[iidx] / data(cidx(iidx+1)-1); |
|
1655 work2[iidx] = tmp; |
|
1656 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1657 { |
|
1658 int idx2 = q_perm[ridx(i)]; |
|
1659 work2[idx2] = work2[idx2] - tmp * data(i); |
|
1660 } |
|
1661 } |
|
1662 } |
|
1663 double atmp = 0; |
|
1664 for (int i = 0; i < j+1; i++) |
|
1665 { |
|
1666 atmp += fabs(work2[i]); |
|
1667 work2[i] = 0.; |
|
1668 } |
|
1669 if (atmp > ainvnorm) |
|
1670 ainvnorm = atmp; |
|
1671 } |
|
1672 } |
|
1673 else |
|
1674 { |
|
1675 retval = b; |
|
1676 Complex *x_vec = retval.fortran_vec (); |
|
1677 |
|
1678 for (int j = 0; j < b_nc; j++) |
|
1679 { |
|
1680 int offset = j * nr; |
|
1681 for (int k = nr-1; k >= 0; k--) |
|
1682 { |
|
1683 if (x_vec[k+offset] != 0.) |
|
1684 { |
|
1685 if (ridx(cidx(k+1)-1) != k) |
|
1686 { |
|
1687 err = -2; |
|
1688 goto triangular_error; |
|
1689 } |
|
1690 |
|
1691 Complex tmp = x_vec[k+offset] / |
|
1692 data(cidx(k+1)-1); |
|
1693 x_vec[k+offset] = tmp; |
|
1694 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1695 { |
|
1696 int iidx = ridx(i); |
|
1697 x_vec[iidx+offset] = |
|
1698 x_vec[iidx+offset] - tmp * data(i); |
|
1699 } |
|
1700 } |
|
1701 } |
|
1702 } |
|
1703 |
|
1704 // Calculation of 1-norm of inv(*this) |
|
1705 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
1706 for (int i = 0; i < nr; i++) |
|
1707 work[i] = 0.; |
|
1708 |
|
1709 for (int j = 0; j < nr; j++) |
|
1710 { |
|
1711 work[j] = 1.; |
|
1712 |
|
1713 for (int k = j; k >= 0; k--) |
|
1714 { |
|
1715 if (work[k] != 0.) |
|
1716 { |
|
1717 double tmp = work[k] / data(cidx(k+1)-1); |
|
1718 work[k] = tmp; |
|
1719 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1720 { |
|
1721 int iidx = ridx(i); |
|
1722 work[iidx] = work[iidx] - tmp * data(i); |
|
1723 } |
|
1724 } |
|
1725 } |
|
1726 double atmp = 0; |
|
1727 for (int i = 0; i < j+1; i++) |
|
1728 { |
|
1729 atmp += fabs(work[i]); |
|
1730 work[i] = 0.; |
|
1731 } |
|
1732 if (atmp > ainvnorm) |
|
1733 ainvnorm = atmp; |
|
1734 } |
|
1735 } |
|
1736 |
|
1737 rcond = 1. / ainvnorm / anorm; |
|
1738 |
|
1739 triangular_error: |
|
1740 if (err != 0) |
|
1741 { |
|
1742 if (sing_handler) |
|
1743 sing_handler (rcond); |
|
1744 else |
|
1745 (*current_liboctave_error_handler) |
|
1746 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
1747 rcond); |
|
1748 } |
|
1749 |
|
1750 volatile double rcond_plus_one = rcond + 1.0; |
|
1751 |
|
1752 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
1753 { |
|
1754 err = -2; |
|
1755 |
|
1756 if (sing_handler) |
|
1757 sing_handler (rcond); |
|
1758 else |
|
1759 (*current_liboctave_error_handler) |
|
1760 ("matrix singular to machine precision, rcond = %g", |
|
1761 rcond); |
|
1762 } |
|
1763 } |
|
1764 else |
|
1765 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
1766 } |
|
1767 |
|
1768 return retval; |
|
1769 } |
|
1770 |
|
1771 SparseComplexMatrix |
|
1772 SparseMatrix::utsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
1773 int& err, double& rcond, |
|
1774 solve_singularity_handler sing_handler) const |
|
1775 { |
|
1776 SparseComplexMatrix retval; |
|
1777 |
|
1778 int nr = rows (); |
|
1779 int nc = cols (); |
|
1780 err = 0; |
|
1781 |
|
1782 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
1783 (*current_liboctave_error_handler) |
|
1784 ("matrix dimension mismatch solution of linear equations"); |
|
1785 else |
|
1786 { |
|
1787 // Print spparms("spumoni") info if requested |
|
1788 int typ = mattype.type (); |
|
1789 mattype.info (); |
|
1790 |
|
1791 if (typ == SparseType::Permuted_Upper || |
|
1792 typ == SparseType::Upper) |
|
1793 { |
|
1794 double anorm = 0.; |
|
1795 double ainvnorm = 0.; |
|
1796 rcond = 0.; |
|
1797 |
|
1798 // Calculate the 1-norm of matrix for rcond calculation |
|
1799 for (int j = 0; j < nr; j++) |
|
1800 { |
|
1801 double atmp = 0.; |
|
1802 for (int i = cidx(j); i < cidx(j+1); i++) |
|
1803 atmp += fabs(data(i)); |
|
1804 if (atmp > anorm) |
|
1805 anorm = atmp; |
|
1806 } |
|
1807 |
|
1808 int b_nr = b.rows (); |
|
1809 int b_nc = b.cols (); |
|
1810 int b_nz = b.nnz (); |
|
1811 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
1812 retval.xcidx(0) = 0; |
|
1813 int ii = 0; |
|
1814 int x_nz = b_nz; |
|
1815 |
|
1816 if (typ == SparseType::Permuted_Upper) |
|
1817 { |
|
1818 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1819 int *p_perm = mattype.triangular_row_perm (); |
|
1820 int *q_perm = mattype.triangular_col_perm (); |
|
1821 |
|
1822 (*current_liboctave_warning_handler) |
|
1823 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
1824 |
|
1825 for (int j = 0; j < b_nc; j++) |
|
1826 { |
|
1827 for (int i = 0; i < nr; i++) |
|
1828 work[i] = 0.; |
|
1829 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1830 work[b.ridx(i)] = b.data(i); |
|
1831 |
|
1832 for (int k = nr-1; k >= 0; k--) |
|
1833 { |
|
1834 int iidx = q_perm[k]; |
|
1835 if (work[iidx] != 0.) |
|
1836 { |
|
1837 if (ridx(cidx(iidx+1)-1) != iidx) |
|
1838 { |
|
1839 err = -2; |
|
1840 goto triangular_error; |
|
1841 } |
|
1842 |
|
1843 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
1844 work[iidx] = tmp; |
|
1845 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1846 { |
|
1847 int idx2 = q_perm[ridx(i)]; |
|
1848 work[idx2] = |
|
1849 work[idx2] - tmp * data(i); |
|
1850 } |
|
1851 } |
|
1852 } |
|
1853 |
|
1854 // Count non-zeros in work vector and adjust space in |
|
1855 // retval if needed |
|
1856 int new_nnz = 0; |
|
1857 for (int i = 0; i < nr; i++) |
|
1858 if (work[i] != 0.) |
|
1859 new_nnz++; |
|
1860 |
|
1861 if (ii + new_nnz > x_nz) |
|
1862 { |
|
1863 // Resize the sparse matrix |
|
1864 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1865 retval.change_capacity (sz); |
|
1866 x_nz = sz; |
|
1867 } |
|
1868 |
|
1869 for (int i = 0; i < nr; i++) |
|
1870 if (work[p_perm[i]] != 0.) |
|
1871 { |
|
1872 retval.xridx(ii) = i; |
|
1873 retval.xdata(ii++) = work[p_perm[i]]; |
|
1874 } |
|
1875 retval.xcidx(j+1) = ii; |
|
1876 } |
|
1877 |
|
1878 retval.maybe_compress (); |
|
1879 |
|
1880 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
1881 // Calculation of 1-norm of inv(*this) |
|
1882 for (int i = 0; i < nr; i++) |
|
1883 work2[i] = 0.; |
|
1884 |
|
1885 for (int j = 0; j < nr; j++) |
|
1886 { |
|
1887 work2[q_perm[j]] = 1.; |
|
1888 |
|
1889 for (int k = j; k >= 0; k--) |
|
1890 { |
|
1891 int iidx = q_perm[k]; |
|
1892 |
|
1893 if (work2[iidx] != 0.) |
|
1894 { |
|
1895 double tmp = work2[iidx] / data(cidx(iidx+1)-1); |
|
1896 work2[iidx] = tmp; |
|
1897 for (int i = cidx(iidx); i < cidx(iidx+1)-1; i++) |
|
1898 { |
|
1899 int idx2 = q_perm[ridx(i)]; |
|
1900 work2[idx2] = work2[idx2] - tmp * data(i); |
|
1901 } |
|
1902 } |
|
1903 } |
|
1904 double atmp = 0; |
|
1905 for (int i = 0; i < j+1; i++) |
|
1906 { |
|
1907 atmp += fabs(work2[i]); |
|
1908 work2[i] = 0.; |
|
1909 } |
|
1910 if (atmp > ainvnorm) |
|
1911 ainvnorm = atmp; |
|
1912 } |
|
1913 } |
|
1914 else |
|
1915 { |
|
1916 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
1917 |
|
1918 for (int j = 0; j < b_nc; j++) |
|
1919 { |
|
1920 for (int i = 0; i < nr; i++) |
|
1921 work[i] = 0.; |
|
1922 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
1923 work[b.ridx(i)] = b.data(i); |
|
1924 |
|
1925 for (int k = nr-1; k >= 0; k--) |
|
1926 { |
|
1927 if (work[k] != 0.) |
|
1928 { |
|
1929 if (ridx(cidx(k+1)-1) != k) |
|
1930 { |
|
1931 err = -2; |
|
1932 goto triangular_error; |
|
1933 } |
|
1934 |
|
1935 Complex tmp = work[k] / data(cidx(k+1)-1); |
|
1936 work[k] = tmp; |
|
1937 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1938 { |
|
1939 int iidx = ridx(i); |
|
1940 work[iidx] = work[iidx] - tmp * data(i); |
|
1941 } |
|
1942 } |
|
1943 } |
|
1944 |
|
1945 // Count non-zeros in work vector and adjust space in |
|
1946 // retval if needed |
|
1947 int new_nnz = 0; |
|
1948 for (int i = 0; i < nr; i++) |
|
1949 if (work[i] != 0.) |
|
1950 new_nnz++; |
|
1951 |
|
1952 if (ii + new_nnz > x_nz) |
|
1953 { |
|
1954 // Resize the sparse matrix |
|
1955 int sz = new_nnz * (b_nc - j) + x_nz; |
|
1956 retval.change_capacity (sz); |
|
1957 x_nz = sz; |
|
1958 } |
|
1959 |
|
1960 for (int i = 0; i < nr; i++) |
|
1961 if (work[i] != 0.) |
|
1962 { |
|
1963 retval.xridx(ii) = i; |
|
1964 retval.xdata(ii++) = work[i]; |
|
1965 } |
|
1966 retval.xcidx(j+1) = ii; |
|
1967 } |
|
1968 |
|
1969 retval.maybe_compress (); |
|
1970 |
|
1971 // Calculation of 1-norm of inv(*this) |
|
1972 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
1973 for (int i = 0; i < nr; i++) |
|
1974 work2[i] = 0.; |
|
1975 |
|
1976 for (int j = 0; j < nr; j++) |
|
1977 { |
|
1978 work2[j] = 1.; |
|
1979 |
|
1980 for (int k = j; k >= 0; k--) |
|
1981 { |
|
1982 if (work2[k] != 0.) |
|
1983 { |
|
1984 double tmp = work2[k] / data(cidx(k+1)-1); |
|
1985 work2[k] = tmp; |
|
1986 for (int i = cidx(k); i < cidx(k+1)-1; i++) |
|
1987 { |
|
1988 int iidx = ridx(i); |
|
1989 work2[iidx] = work2[iidx] - tmp * data(i); |
|
1990 } |
|
1991 } |
|
1992 } |
|
1993 double atmp = 0; |
|
1994 for (int i = 0; i < j+1; i++) |
|
1995 { |
|
1996 atmp += fabs(work2[i]); |
|
1997 work2[i] = 0.; |
|
1998 } |
|
1999 if (atmp > ainvnorm) |
|
2000 ainvnorm = atmp; |
|
2001 } |
|
2002 } |
|
2003 |
|
2004 rcond = 1. / ainvnorm / anorm; |
|
2005 |
|
2006 triangular_error: |
|
2007 if (err != 0) |
|
2008 { |
|
2009 if (sing_handler) |
|
2010 sing_handler (rcond); |
|
2011 else |
|
2012 (*current_liboctave_error_handler) |
|
2013 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2014 rcond); |
|
2015 } |
|
2016 |
|
2017 volatile double rcond_plus_one = rcond + 1.0; |
|
2018 |
|
2019 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2020 { |
|
2021 err = -2; |
|
2022 |
|
2023 if (sing_handler) |
|
2024 sing_handler (rcond); |
|
2025 else |
|
2026 (*current_liboctave_error_handler) |
|
2027 ("matrix singular to machine precision, rcond = %g", |
|
2028 rcond); |
|
2029 } |
|
2030 } |
|
2031 else |
|
2032 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2033 } |
|
2034 |
|
2035 return retval; |
|
2036 } |
|
2037 |
|
2038 Matrix |
|
2039 SparseMatrix::ltsolve (SparseType &mattype, const Matrix& b, int& err, |
|
2040 double& rcond, |
|
2041 solve_singularity_handler sing_handler) const |
|
2042 { |
|
2043 Matrix retval; |
|
2044 |
|
2045 int nr = rows (); |
|
2046 int nc = cols (); |
|
2047 err = 0; |
|
2048 |
|
2049 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2050 (*current_liboctave_error_handler) |
|
2051 ("matrix dimension mismatch solution of linear equations"); |
|
2052 else |
|
2053 { |
|
2054 // Print spparms("spumoni") info if requested |
|
2055 int typ = mattype.type (); |
|
2056 mattype.info (); |
|
2057 |
|
2058 if (typ == SparseType::Permuted_Lower || |
|
2059 typ == SparseType::Lower) |
|
2060 { |
|
2061 double anorm = 0.; |
|
2062 double ainvnorm = 0.; |
|
2063 int b_cols = b.cols (); |
|
2064 rcond = 0.; |
|
2065 |
|
2066 // Calculate the 1-norm of matrix for rcond calculation |
|
2067 for (int j = 0; j < nr; j++) |
|
2068 { |
|
2069 double atmp = 0.; |
|
2070 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2071 atmp += fabs(data(i)); |
|
2072 if (atmp > anorm) |
|
2073 anorm = atmp; |
|
2074 } |
|
2075 |
|
2076 if (typ == SparseType::Permuted_Lower) |
|
2077 { |
|
2078 retval.resize (b.rows (), b.cols ()); |
|
2079 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
2080 int *p_perm = mattype.triangular_row_perm (); |
|
2081 int *q_perm = mattype.triangular_col_perm (); |
|
2082 |
|
2083 (*current_liboctave_warning_handler) |
|
2084 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2085 |
|
2086 for (int j = 0; j < b_cols; j++) |
|
2087 { |
|
2088 for (int i = 0; i < nr; i++) |
|
2089 work[i] = b(i,j); |
|
2090 |
|
2091 for (int k = 0; k < nr; k++) |
|
2092 { |
|
2093 int iidx = q_perm[k]; |
|
2094 if (work[iidx] != 0.) |
|
2095 { |
|
2096 if (ridx(cidx(iidx)) != iidx) |
|
2097 { |
|
2098 err = -2; |
|
2099 goto triangular_error; |
|
2100 } |
|
2101 |
|
2102 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2103 work[iidx] = tmp; |
|
2104 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2105 { |
|
2106 int idx2 = q_perm[ridx(i)]; |
|
2107 work[idx2] = |
|
2108 work[idx2] - tmp * data(i); |
|
2109 } |
|
2110 } |
|
2111 } |
|
2112 |
|
2113 for (int i = 0; i < nr; i++) |
|
2114 retval (i, j) = work[p_perm[i]]; |
|
2115 |
|
2116 } |
|
2117 |
|
2118 // Calculation of 1-norm of inv(*this) |
|
2119 for (int i = 0; i < nr; i++) |
|
2120 work[i] = 0.; |
|
2121 |
|
2122 for (int j = 0; j < nr; j++) |
|
2123 { |
|
2124 work[q_perm[j]] = 1.; |
|
2125 |
|
2126 for (int k = 0; k < nr; k++) |
|
2127 { |
|
2128 int iidx = q_perm[k]; |
|
2129 |
|
2130 if (work[iidx] != 0.) |
|
2131 { |
|
2132 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2133 work[iidx] = tmp; |
|
2134 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2135 { |
|
2136 int idx2 = q_perm[ridx(i)]; |
|
2137 work[idx2] = work[idx2] - tmp * data(i); |
|
2138 } |
|
2139 } |
|
2140 } |
|
2141 double atmp = 0; |
|
2142 for (int i = 0; i < j+1; i++) |
|
2143 { |
|
2144 atmp += fabs(work[i]); |
|
2145 work[i] = 0.; |
|
2146 } |
|
2147 if (atmp > ainvnorm) |
|
2148 ainvnorm = atmp; |
|
2149 } |
|
2150 } |
|
2151 else |
|
2152 { |
|
2153 retval = b; |
|
2154 double *x_vec = retval.fortran_vec (); |
|
2155 |
|
2156 for (int j = 0; j < b_cols; j++) |
|
2157 { |
|
2158 int offset = j * nr; |
|
2159 for (int k = 0; k < nr; k++) |
|
2160 { |
|
2161 if (x_vec[k+offset] != 0.) |
|
2162 { |
|
2163 if (ridx(cidx(k)) != k) |
|
2164 { |
|
2165 err = -2; |
|
2166 goto triangular_error; |
|
2167 } |
|
2168 |
|
2169 double tmp = x_vec[k+offset] / |
|
2170 data(cidx(k)); |
|
2171 x_vec[k+offset] = tmp; |
|
2172 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2173 { |
|
2174 int iidx = ridx(i); |
|
2175 x_vec[iidx+offset] = |
|
2176 x_vec[iidx+offset] - tmp * data(i); |
|
2177 } |
|
2178 } |
|
2179 } |
|
2180 } |
|
2181 |
|
2182 // Calculation of 1-norm of inv(*this) |
|
2183 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
2184 for (int i = 0; i < nr; i++) |
|
2185 work[i] = 0.; |
|
2186 |
|
2187 for (int j = 0; j < nr; j++) |
|
2188 { |
|
2189 work[j] = 1.; |
|
2190 |
|
2191 for (int k = j; k < nr; k++) |
|
2192 { |
|
2193 |
|
2194 if (work[k] != 0.) |
|
2195 { |
|
2196 double tmp = work[k] / data(cidx(k)); |
|
2197 work[k] = tmp; |
|
2198 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2199 { |
|
2200 int iidx = ridx(i); |
|
2201 work[iidx] = work[iidx] - tmp * data(i); |
|
2202 } |
|
2203 } |
|
2204 } |
|
2205 double atmp = 0; |
|
2206 for (int i = j; i < nr; i++) |
|
2207 { |
|
2208 atmp += fabs(work[i]); |
|
2209 work[i] = 0.; |
|
2210 } |
|
2211 if (atmp > ainvnorm) |
|
2212 ainvnorm = atmp; |
|
2213 } |
|
2214 |
|
2215 } |
|
2216 |
|
2217 rcond = 1. / ainvnorm / anorm; |
|
2218 |
|
2219 triangular_error: |
|
2220 if (err != 0) |
|
2221 { |
|
2222 if (sing_handler) |
|
2223 sing_handler (rcond); |
|
2224 else |
|
2225 (*current_liboctave_error_handler) |
|
2226 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2227 rcond); |
|
2228 } |
|
2229 |
|
2230 volatile double rcond_plus_one = rcond + 1.0; |
|
2231 |
|
2232 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2233 { |
|
2234 err = -2; |
|
2235 |
|
2236 if (sing_handler) |
|
2237 sing_handler (rcond); |
|
2238 else |
|
2239 (*current_liboctave_error_handler) |
|
2240 ("matrix singular to machine precision, rcond = %g", |
|
2241 rcond); |
|
2242 } |
|
2243 } |
|
2244 else |
|
2245 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2246 } |
|
2247 |
|
2248 return retval; |
|
2249 } |
|
2250 |
|
2251 SparseMatrix |
|
2252 SparseMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, int& err, |
|
2253 double& rcond, solve_singularity_handler sing_handler) const |
|
2254 { |
|
2255 SparseMatrix retval; |
|
2256 |
|
2257 int nr = rows (); |
|
2258 int nc = cols (); |
|
2259 err = 0; |
|
2260 |
|
2261 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2262 (*current_liboctave_error_handler) |
|
2263 ("matrix dimension mismatch solution of linear equations"); |
|
2264 else |
|
2265 { |
|
2266 // Print spparms("spumoni") info if requested |
|
2267 int typ = mattype.type (); |
|
2268 mattype.info (); |
|
2269 |
|
2270 if (typ == SparseType::Permuted_Lower || |
|
2271 typ == SparseType::Lower) |
|
2272 { |
|
2273 double anorm = 0.; |
|
2274 double ainvnorm = 0.; |
|
2275 rcond = 0.; |
|
2276 |
|
2277 // Calculate the 1-norm of matrix for rcond calculation |
|
2278 for (int j = 0; j < nr; j++) |
|
2279 { |
|
2280 double atmp = 0.; |
|
2281 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2282 atmp += fabs(data(i)); |
|
2283 if (atmp > anorm) |
|
2284 anorm = atmp; |
|
2285 } |
|
2286 |
|
2287 int b_nr = b.rows (); |
|
2288 int b_nc = b.cols (); |
|
2289 int b_nz = b.nnz (); |
|
2290 retval = SparseMatrix (b_nr, b_nc, b_nz); |
|
2291 retval.xcidx(0) = 0; |
|
2292 int ii = 0; |
|
2293 int x_nz = b_nz; |
|
2294 |
|
2295 if (typ == SparseType::Permuted_Lower) |
|
2296 { |
|
2297 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
2298 int *p_perm = mattype.triangular_row_perm (); |
|
2299 int *q_perm = mattype.triangular_col_perm (); |
|
2300 |
|
2301 (*current_liboctave_warning_handler) |
|
2302 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2303 |
|
2304 for (int j = 0; j < b_nc; j++) |
|
2305 { |
|
2306 for (int i = 0; i < nr; i++) |
|
2307 work[i] = 0.; |
|
2308 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2309 work[b.ridx(i)] = b.data(i); |
|
2310 |
|
2311 for (int k = 0; k < nr; k++) |
|
2312 { |
|
2313 int iidx = q_perm[k]; |
|
2314 if (work[iidx] != 0.) |
|
2315 { |
|
2316 if (ridx(cidx(iidx)) != iidx) |
|
2317 { |
|
2318 err = -2; |
|
2319 goto triangular_error; |
|
2320 } |
|
2321 |
|
2322 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2323 work[iidx] = tmp; |
|
2324 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2325 { |
|
2326 int idx2 = q_perm[ridx(i)]; |
|
2327 work[idx2] = |
|
2328 work[idx2] - tmp * data(i); |
|
2329 } |
|
2330 } |
|
2331 } |
|
2332 |
|
2333 // Count non-zeros in work vector and adjust space in |
|
2334 // retval if needed |
|
2335 int new_nnz = 0; |
|
2336 for (int i = 0; i < nr; i++) |
|
2337 if (work[i] != 0.) |
|
2338 new_nnz++; |
|
2339 |
|
2340 if (ii + new_nnz > x_nz) |
|
2341 { |
|
2342 // Resize the sparse matrix |
|
2343 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2344 retval.change_capacity (sz); |
|
2345 x_nz = sz; |
|
2346 } |
|
2347 |
|
2348 for (int i = 0; i < nr; i++) |
|
2349 if (work[p_perm[i]] != 0.) |
|
2350 { |
|
2351 retval.xridx(ii) = i; |
|
2352 retval.xdata(ii++) = work[p_perm[i]]; |
|
2353 } |
|
2354 retval.xcidx(j+1) = ii; |
|
2355 } |
|
2356 |
|
2357 retval.maybe_compress (); |
|
2358 |
|
2359 // Calculation of 1-norm of inv(*this) |
|
2360 for (int i = 0; i < nr; i++) |
|
2361 work[i] = 0.; |
|
2362 |
|
2363 for (int j = 0; j < nr; j++) |
|
2364 { |
|
2365 work[q_perm[j]] = 1.; |
|
2366 |
|
2367 for (int k = 0; k < nr; k++) |
|
2368 { |
|
2369 int iidx = q_perm[k]; |
|
2370 |
|
2371 if (work[iidx] != 0.) |
|
2372 { |
|
2373 double tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2374 work[iidx] = tmp; |
|
2375 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2376 { |
|
2377 int idx2 = q_perm[ridx(i)]; |
|
2378 work[idx2] = work[idx2] - tmp * data(i); |
|
2379 } |
|
2380 } |
|
2381 } |
|
2382 double atmp = 0; |
|
2383 for (int i = 0; i < j+1; i++) |
|
2384 { |
|
2385 atmp += fabs(work[i]); |
|
2386 work[i] = 0.; |
|
2387 } |
|
2388 if (atmp > ainvnorm) |
|
2389 ainvnorm = atmp; |
|
2390 } |
|
2391 } |
|
2392 else |
|
2393 { |
|
2394 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
2395 |
|
2396 for (int j = 0; j < b_nc; j++) |
|
2397 { |
|
2398 for (int i = 0; i < nr; i++) |
|
2399 work[i] = 0.; |
|
2400 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2401 work[b.ridx(i)] = b.data(i); |
|
2402 |
|
2403 for (int k = 0; k < nr; k++) |
|
2404 { |
|
2405 if (work[k] != 0.) |
|
2406 { |
|
2407 if (ridx(cidx(k)) != k) |
|
2408 { |
|
2409 err = -2; |
|
2410 goto triangular_error; |
|
2411 } |
|
2412 |
|
2413 double tmp = work[k] / data(cidx(k)); |
|
2414 work[k] = tmp; |
|
2415 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2416 { |
|
2417 int iidx = ridx(i); |
|
2418 work[iidx] = work[iidx] - tmp * data(i); |
|
2419 } |
|
2420 } |
|
2421 } |
|
2422 |
|
2423 // Count non-zeros in work vector and adjust space in |
|
2424 // retval if needed |
|
2425 int new_nnz = 0; |
|
2426 for (int i = 0; i < nr; i++) |
|
2427 if (work[i] != 0.) |
|
2428 new_nnz++; |
|
2429 |
|
2430 if (ii + new_nnz > x_nz) |
|
2431 { |
|
2432 // Resize the sparse matrix |
|
2433 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2434 retval.change_capacity (sz); |
|
2435 x_nz = sz; |
|
2436 } |
|
2437 |
|
2438 for (int i = 0; i < nr; i++) |
|
2439 if (work[i] != 0.) |
|
2440 { |
|
2441 retval.xridx(ii) = i; |
|
2442 retval.xdata(ii++) = work[i]; |
|
2443 } |
|
2444 retval.xcidx(j+1) = ii; |
|
2445 } |
|
2446 |
|
2447 retval.maybe_compress (); |
|
2448 |
|
2449 // Calculation of 1-norm of inv(*this) |
|
2450 for (int i = 0; i < nr; i++) |
|
2451 work[i] = 0.; |
|
2452 |
|
2453 for (int j = 0; j < nr; j++) |
|
2454 { |
|
2455 work[j] = 1.; |
|
2456 |
|
2457 for (int k = j; k < nr; k++) |
|
2458 { |
|
2459 |
|
2460 if (work[k] != 0.) |
|
2461 { |
|
2462 double tmp = work[k] / data(cidx(k)); |
|
2463 work[k] = tmp; |
|
2464 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2465 { |
|
2466 int iidx = ridx(i); |
|
2467 work[iidx] = work[iidx] - tmp * data(i); |
|
2468 } |
|
2469 } |
|
2470 } |
|
2471 double atmp = 0; |
|
2472 for (int i = j; i < nr; i++) |
|
2473 { |
|
2474 atmp += fabs(work[i]); |
|
2475 work[i] = 0.; |
|
2476 } |
|
2477 if (atmp > ainvnorm) |
|
2478 ainvnorm = atmp; |
|
2479 } |
|
2480 |
|
2481 } |
|
2482 |
|
2483 rcond = 1. / ainvnorm / anorm; |
|
2484 |
|
2485 triangular_error: |
|
2486 if (err != 0) |
|
2487 { |
|
2488 if (sing_handler) |
|
2489 sing_handler (rcond); |
|
2490 else |
|
2491 (*current_liboctave_error_handler) |
|
2492 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2493 rcond); |
|
2494 } |
|
2495 |
|
2496 volatile double rcond_plus_one = rcond + 1.0; |
|
2497 |
|
2498 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2499 { |
|
2500 err = -2; |
|
2501 |
|
2502 if (sing_handler) |
|
2503 sing_handler (rcond); |
|
2504 else |
|
2505 (*current_liboctave_error_handler) |
|
2506 ("matrix singular to machine precision, rcond = %g", |
|
2507 rcond); |
|
2508 } |
|
2509 } |
|
2510 else |
|
2511 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2512 } |
|
2513 |
|
2514 return retval; |
|
2515 } |
|
2516 |
|
2517 ComplexMatrix |
|
2518 SparseMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b, int& err, |
|
2519 double& rcond, solve_singularity_handler sing_handler) const |
|
2520 { |
|
2521 ComplexMatrix retval; |
|
2522 |
|
2523 int nr = rows (); |
|
2524 int nc = cols (); |
|
2525 err = 0; |
|
2526 |
|
2527 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2528 (*current_liboctave_error_handler) |
|
2529 ("matrix dimension mismatch solution of linear equations"); |
|
2530 else |
|
2531 { |
|
2532 // Print spparms("spumoni") info if requested |
|
2533 int typ = mattype.type (); |
|
2534 mattype.info (); |
|
2535 |
|
2536 if (typ == SparseType::Permuted_Lower || |
|
2537 typ == SparseType::Lower) |
|
2538 { |
|
2539 double anorm = 0.; |
|
2540 double ainvnorm = 0.; |
|
2541 int b_nc = b.cols (); |
|
2542 rcond = 0.; |
|
2543 |
|
2544 // Calculate the 1-norm of matrix for rcond calculation |
|
2545 for (int j = 0; j < nr; j++) |
|
2546 { |
|
2547 double atmp = 0.; |
|
2548 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2549 atmp += fabs(data(i)); |
|
2550 if (atmp > anorm) |
|
2551 anorm = atmp; |
|
2552 } |
|
2553 |
|
2554 if (typ == SparseType::Permuted_Lower) |
|
2555 { |
|
2556 retval.resize (b.rows (), b.cols ()); |
|
2557 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2558 int *p_perm = mattype.triangular_row_perm (); |
|
2559 int *q_perm = mattype.triangular_col_perm (); |
|
2560 |
|
2561 (*current_liboctave_warning_handler) |
|
2562 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2563 |
|
2564 for (int j = 0; j < b_nc; j++) |
|
2565 { |
|
2566 for (int i = 0; i < nr; i++) |
|
2567 work[i] = b(i,j); |
|
2568 |
|
2569 for (int k = 0; k < nr; k++) |
|
2570 { |
|
2571 int iidx = q_perm[k]; |
|
2572 if (work[iidx] != 0.) |
|
2573 { |
|
2574 if (ridx(cidx(iidx)) != iidx) |
|
2575 { |
|
2576 err = -2; |
|
2577 goto triangular_error; |
|
2578 } |
|
2579 |
|
2580 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2581 work[iidx] = tmp; |
|
2582 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2583 { |
|
2584 int idx2 = q_perm[ridx(i)]; |
|
2585 work[idx2] = |
|
2586 work[idx2] - tmp * data(i); |
|
2587 } |
|
2588 } |
|
2589 } |
|
2590 |
|
2591 for (int i = 0; i < nr; i++) |
|
2592 retval (i, j) = work[p_perm[i]]; |
|
2593 |
|
2594 } |
|
2595 |
|
2596 // Calculation of 1-norm of inv(*this) |
|
2597 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
2598 for (int i = 0; i < nr; i++) |
|
2599 work2[i] = 0.; |
|
2600 |
|
2601 for (int j = 0; j < nr; j++) |
|
2602 { |
|
2603 work2[q_perm[j]] = 1.; |
|
2604 |
|
2605 for (int k = 0; k < nr; k++) |
|
2606 { |
|
2607 int iidx = q_perm[k]; |
|
2608 |
|
2609 if (work2[iidx] != 0.) |
|
2610 { |
|
2611 double tmp = work2[iidx] / data(cidx(iidx+1)-1); |
|
2612 work2[iidx] = tmp; |
|
2613 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2614 { |
|
2615 int idx2 = q_perm[ridx(i)]; |
|
2616 work2[idx2] = work2[idx2] - tmp * data(i); |
|
2617 } |
|
2618 } |
|
2619 } |
|
2620 double atmp = 0; |
|
2621 for (int i = 0; i < j+1; i++) |
|
2622 { |
|
2623 atmp += fabs(work2[i]); |
|
2624 work2[i] = 0.; |
|
2625 } |
|
2626 if (atmp > ainvnorm) |
|
2627 ainvnorm = atmp; |
|
2628 } |
|
2629 } |
|
2630 else |
|
2631 { |
|
2632 retval = b; |
|
2633 Complex *x_vec = retval.fortran_vec (); |
|
2634 |
|
2635 for (int j = 0; j < b_nc; j++) |
|
2636 { |
|
2637 int offset = j * nr; |
|
2638 for (int k = 0; k < nr; k++) |
|
2639 { |
|
2640 if (x_vec[k+offset] != 0.) |
|
2641 { |
|
2642 if (ridx(cidx(k)) != k) |
|
2643 { |
|
2644 err = -2; |
|
2645 goto triangular_error; |
|
2646 } |
|
2647 |
|
2648 Complex tmp = x_vec[k+offset] / |
|
2649 data(cidx(k)); |
|
2650 x_vec[k+offset] = tmp; |
|
2651 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2652 { |
|
2653 int iidx = ridx(i); |
|
2654 x_vec[iidx+offset] = |
|
2655 x_vec[iidx+offset] - tmp * data(i); |
|
2656 } |
|
2657 } |
|
2658 } |
|
2659 } |
|
2660 |
|
2661 // Calculation of 1-norm of inv(*this) |
|
2662 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
2663 for (int i = 0; i < nr; i++) |
|
2664 work[i] = 0.; |
|
2665 |
|
2666 for (int j = 0; j < nr; j++) |
|
2667 { |
|
2668 work[j] = 1.; |
|
2669 |
|
2670 for (int k = j; k < nr; k++) |
|
2671 { |
|
2672 |
|
2673 if (work[k] != 0.) |
|
2674 { |
|
2675 double tmp = work[k] / data(cidx(k)); |
|
2676 work[k] = tmp; |
|
2677 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2678 { |
|
2679 int iidx = ridx(i); |
|
2680 work[iidx] = work[iidx] - tmp * data(i); |
|
2681 } |
|
2682 } |
|
2683 } |
|
2684 double atmp = 0; |
|
2685 for (int i = j; i < nr; i++) |
|
2686 { |
|
2687 atmp += fabs(work[i]); |
|
2688 work[i] = 0.; |
|
2689 } |
|
2690 if (atmp > ainvnorm) |
|
2691 ainvnorm = atmp; |
|
2692 } |
|
2693 |
|
2694 } |
|
2695 |
|
2696 rcond = 1. / ainvnorm / anorm; |
|
2697 |
|
2698 triangular_error: |
|
2699 if (err != 0) |
|
2700 { |
|
2701 if (sing_handler) |
|
2702 sing_handler (rcond); |
|
2703 else |
|
2704 (*current_liboctave_error_handler) |
|
2705 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2706 rcond); |
|
2707 } |
|
2708 |
|
2709 volatile double rcond_plus_one = rcond + 1.0; |
|
2710 |
|
2711 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2712 { |
|
2713 err = -2; |
|
2714 |
|
2715 if (sing_handler) |
|
2716 sing_handler (rcond); |
|
2717 else |
|
2718 (*current_liboctave_error_handler) |
|
2719 ("matrix singular to machine precision, rcond = %g", |
|
2720 rcond); |
|
2721 } |
|
2722 } |
|
2723 else |
|
2724 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2725 } |
|
2726 |
|
2727 return retval; |
|
2728 } |
|
2729 |
|
2730 SparseComplexMatrix |
|
2731 SparseMatrix::ltsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
2732 int& err, double& rcond, |
|
2733 solve_singularity_handler sing_handler) const |
|
2734 { |
|
2735 SparseComplexMatrix retval; |
|
2736 |
|
2737 int nr = rows (); |
|
2738 int nc = cols (); |
|
2739 err = 0; |
|
2740 |
|
2741 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
2742 (*current_liboctave_error_handler) |
|
2743 ("matrix dimension mismatch solution of linear equations"); |
|
2744 else |
|
2745 { |
|
2746 // Print spparms("spumoni") info if requested |
|
2747 int typ = mattype.type (); |
|
2748 mattype.info (); |
|
2749 |
|
2750 if (typ == SparseType::Permuted_Lower || |
|
2751 typ == SparseType::Lower) |
|
2752 { |
|
2753 double anorm = 0.; |
|
2754 double ainvnorm = 0.; |
|
2755 rcond = 0.; |
|
2756 |
|
2757 // Calculate the 1-norm of matrix for rcond calculation |
|
2758 for (int j = 0; j < nr; j++) |
|
2759 { |
|
2760 double atmp = 0.; |
|
2761 for (int i = cidx(j); i < cidx(j+1); i++) |
|
2762 atmp += fabs(data(i)); |
|
2763 if (atmp > anorm) |
|
2764 anorm = atmp; |
|
2765 } |
|
2766 |
|
2767 int b_nr = b.rows (); |
|
2768 int b_nc = b.cols (); |
|
2769 int b_nz = b.nnz (); |
|
2770 retval = SparseComplexMatrix (b_nr, b_nc, b_nz); |
|
2771 retval.xcidx(0) = 0; |
|
2772 int ii = 0; |
|
2773 int x_nz = b_nz; |
|
2774 |
|
2775 if (typ == SparseType::Permuted_Lower) |
|
2776 { |
|
2777 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2778 int *p_perm = mattype.triangular_row_perm (); |
|
2779 int *q_perm = mattype.triangular_col_perm (); |
|
2780 |
|
2781 (*current_liboctave_warning_handler) |
|
2782 ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); |
|
2783 |
|
2784 for (int j = 0; j < b_nc; j++) |
|
2785 { |
|
2786 for (int i = 0; i < nr; i++) |
|
2787 work[i] = 0.; |
|
2788 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2789 work[b.ridx(i)] = b.data(i); |
|
2790 |
|
2791 for (int k = 0; k < nr; k++) |
|
2792 { |
|
2793 int iidx = q_perm[k]; |
|
2794 if (work[iidx] != 0.) |
|
2795 { |
|
2796 if (ridx(cidx(iidx)) != iidx) |
|
2797 { |
|
2798 err = -2; |
|
2799 goto triangular_error; |
|
2800 } |
|
2801 |
|
2802 Complex tmp = work[iidx] / data(cidx(iidx+1)-1); |
|
2803 work[iidx] = tmp; |
|
2804 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2805 { |
|
2806 int idx2 = q_perm[ridx(i)]; |
|
2807 work[idx2] = |
|
2808 work[idx2] - tmp * data(i); |
|
2809 } |
|
2810 } |
|
2811 } |
|
2812 |
|
2813 // Count non-zeros in work vector and adjust space in |
|
2814 // retval if needed |
|
2815 int new_nnz = 0; |
|
2816 for (int i = 0; i < nr; i++) |
|
2817 if (work[i] != 0.) |
|
2818 new_nnz++; |
|
2819 |
|
2820 if (ii + new_nnz > x_nz) |
|
2821 { |
|
2822 // Resize the sparse matrix |
|
2823 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2824 retval.change_capacity (sz); |
|
2825 x_nz = sz; |
|
2826 } |
|
2827 |
|
2828 for (int i = 0; i < nr; i++) |
|
2829 if (work[p_perm[i]] != 0.) |
|
2830 { |
|
2831 retval.xridx(ii) = i; |
|
2832 retval.xdata(ii++) = work[p_perm[i]]; |
|
2833 } |
|
2834 retval.xcidx(j+1) = ii; |
|
2835 } |
|
2836 |
|
2837 retval.maybe_compress (); |
|
2838 |
|
2839 // Calculation of 1-norm of inv(*this) |
|
2840 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
2841 for (int i = 0; i < nr; i++) |
|
2842 work2[i] = 0.; |
|
2843 |
|
2844 for (int j = 0; j < nr; j++) |
|
2845 { |
|
2846 work2[q_perm[j]] = 1.; |
|
2847 |
|
2848 for (int k = 0; k < nr; k++) |
|
2849 { |
|
2850 int iidx = q_perm[k]; |
|
2851 |
|
2852 if (work2[iidx] != 0.) |
|
2853 { |
|
2854 double tmp = work2[iidx] / data(cidx(iidx+1)-1); |
|
2855 work2[iidx] = tmp; |
|
2856 for (int i = cidx(iidx)+1; i < cidx(iidx+1); i++) |
|
2857 { |
|
2858 int idx2 = q_perm[ridx(i)]; |
|
2859 work2[idx2] = work2[idx2] - tmp * data(i); |
|
2860 } |
|
2861 } |
|
2862 } |
|
2863 double atmp = 0; |
|
2864 for (int i = 0; i < j+1; i++) |
|
2865 { |
|
2866 atmp += fabs(work2[i]); |
|
2867 work2[i] = 0.; |
|
2868 } |
|
2869 if (atmp > ainvnorm) |
|
2870 ainvnorm = atmp; |
|
2871 } |
|
2872 } |
|
2873 else |
|
2874 { |
|
2875 OCTAVE_LOCAL_BUFFER (Complex, work, nr); |
|
2876 |
|
2877 for (int j = 0; j < b_nc; j++) |
|
2878 { |
|
2879 for (int i = 0; i < nr; i++) |
|
2880 work[i] = 0.; |
|
2881 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
2882 work[b.ridx(i)] = b.data(i); |
|
2883 |
|
2884 for (int k = 0; k < nr; k++) |
|
2885 { |
|
2886 if (work[k] != 0.) |
|
2887 { |
|
2888 if (ridx(cidx(k)) != k) |
|
2889 { |
|
2890 err = -2; |
|
2891 goto triangular_error; |
|
2892 } |
|
2893 |
|
2894 Complex tmp = work[k] / data(cidx(k)); |
|
2895 work[k] = tmp; |
|
2896 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2897 { |
|
2898 int iidx = ridx(i); |
|
2899 work[iidx] = work[iidx] - tmp * data(i); |
|
2900 } |
|
2901 } |
|
2902 } |
|
2903 |
|
2904 // Count non-zeros in work vector and adjust space in |
|
2905 // retval if needed |
|
2906 int new_nnz = 0; |
|
2907 for (int i = 0; i < nr; i++) |
|
2908 if (work[i] != 0.) |
|
2909 new_nnz++; |
|
2910 |
|
2911 if (ii + new_nnz > x_nz) |
|
2912 { |
|
2913 // Resize the sparse matrix |
|
2914 int sz = new_nnz * (b_nc - j) + x_nz; |
|
2915 retval.change_capacity (sz); |
|
2916 x_nz = sz; |
|
2917 } |
|
2918 |
|
2919 for (int i = 0; i < nr; i++) |
|
2920 if (work[i] != 0.) |
|
2921 { |
|
2922 retval.xridx(ii) = i; |
|
2923 retval.xdata(ii++) = work[i]; |
|
2924 } |
|
2925 retval.xcidx(j+1) = ii; |
|
2926 } |
|
2927 |
|
2928 retval.maybe_compress (); |
|
2929 |
|
2930 // Calculation of 1-norm of inv(*this) |
|
2931 OCTAVE_LOCAL_BUFFER (double, work2, nr); |
|
2932 for (int i = 0; i < nr; i++) |
|
2933 work2[i] = 0.; |
|
2934 |
|
2935 for (int j = 0; j < nr; j++) |
|
2936 { |
|
2937 work2[j] = 1.; |
|
2938 |
|
2939 for (int k = j; k < nr; k++) |
|
2940 { |
|
2941 |
|
2942 if (work2[k] != 0.) |
|
2943 { |
|
2944 double tmp = work2[k] / data(cidx(k)); |
|
2945 work2[k] = tmp; |
|
2946 for (int i = cidx(k)+1; i < cidx(k+1); i++) |
|
2947 { |
|
2948 int iidx = ridx(i); |
|
2949 work2[iidx] = work2[iidx] - tmp * data(i); |
|
2950 } |
|
2951 } |
|
2952 } |
|
2953 double atmp = 0; |
|
2954 for (int i = j; i < nr; i++) |
|
2955 { |
|
2956 atmp += fabs(work2[i]); |
|
2957 work2[i] = 0.; |
|
2958 } |
|
2959 if (atmp > ainvnorm) |
|
2960 ainvnorm = atmp; |
|
2961 } |
|
2962 |
|
2963 } |
|
2964 |
|
2965 rcond = 1. / ainvnorm / anorm; |
|
2966 |
|
2967 triangular_error: |
|
2968 if (err != 0) |
|
2969 { |
|
2970 if (sing_handler) |
|
2971 sing_handler (rcond); |
|
2972 else |
|
2973 (*current_liboctave_error_handler) |
|
2974 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
2975 rcond); |
|
2976 } |
|
2977 |
|
2978 volatile double rcond_plus_one = rcond + 1.0; |
|
2979 |
|
2980 if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
2981 { |
|
2982 err = -2; |
|
2983 |
|
2984 if (sing_handler) |
|
2985 sing_handler (rcond); |
|
2986 else |
|
2987 (*current_liboctave_error_handler) |
|
2988 ("matrix singular to machine precision, rcond = %g", |
|
2989 rcond); |
|
2990 } |
|
2991 } |
|
2992 else |
|
2993 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
2994 } |
|
2995 |
|
2996 return retval; |
|
2997 } |
|
2998 |
|
2999 Matrix |
|
3000 SparseMatrix::trisolve (SparseType &mattype, const Matrix& b, int& err, |
|
3001 double& rcond, |
|
3002 solve_singularity_handler sing_handler) const |
|
3003 { |
|
3004 Matrix retval; |
|
3005 |
|
3006 int nr = rows (); |
|
3007 int nc = cols (); |
|
3008 err = 0; |
|
3009 |
|
3010 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3011 (*current_liboctave_error_handler) |
|
3012 ("matrix dimension mismatch solution of linear equations"); |
|
3013 else |
|
3014 { |
|
3015 // Print spparms("spumoni") info if requested |
|
3016 volatile int typ = mattype.type (); |
|
3017 mattype.info (); |
|
3018 |
|
3019 if (typ == SparseType::Tridiagonal_Hermitian) |
|
3020 { |
|
3021 OCTAVE_LOCAL_BUFFER (double, D, nr); |
|
3022 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1); |
|
3023 |
|
3024 if (mattype.is_dense ()) |
|
3025 { |
|
3026 int ii = 0; |
|
3027 |
|
3028 for (int j = 0; j < nc-1; j++) |
|
3029 { |
|
3030 D[j] = data(ii++); |
|
3031 DL[j] = data(ii); |
|
3032 ii += 2; |
|
3033 } |
|
3034 D[nc-1] = data(ii); |
|
3035 } |
|
3036 else |
|
3037 { |
|
3038 D[0] = 0.; |
|
3039 for (int i = 0; i < nr - 1; i++) |
|
3040 { |
|
3041 D[i+1] = 0.; |
|
3042 DL[i] = 0.; |
|
3043 } |
|
3044 |
|
3045 for (int j = 0; j < nc; j++) |
|
3046 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3047 { |
|
3048 if (ridx(i) == j) |
|
3049 D[j] = data(i); |
|
3050 else if (ridx(i) == j + 1) |
|
3051 DL[j] = data(i); |
|
3052 } |
|
3053 } |
|
3054 |
|
3055 int b_nc = b.cols(); |
|
3056 retval = b; |
|
3057 double *result = retval.fortran_vec (); |
|
3058 |
|
3059 F77_XFCN (dptsv, DPTSV, (nr, b_nc, D, DL, result, |
|
3060 b.rows(), err)); |
|
3061 |
|
3062 if (f77_exception_encountered) |
|
3063 (*current_liboctave_error_handler) |
|
3064 ("unrecoverable error in dptsv"); |
|
3065 else if (err != 0) |
|
3066 { |
|
3067 err = 0; |
|
3068 mattype.mark_as_unsymmetric (); |
|
3069 typ = SparseType::Tridiagonal; |
|
3070 } |
|
3071 else |
|
3072 rcond = 1.; |
|
3073 } |
|
3074 |
|
3075 if (typ == SparseType::Tridiagonal) |
|
3076 { |
|
3077 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1); |
|
3078 OCTAVE_LOCAL_BUFFER (double, D, nr); |
|
3079 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1); |
|
3080 |
|
3081 if (mattype.is_dense ()) |
|
3082 { |
|
3083 int ii = 0; |
|
3084 |
|
3085 for (int j = 0; j < nc-1; j++) |
|
3086 { |
|
3087 D[j] = data(ii++); |
|
3088 DL[j] = data(ii++); |
|
3089 DU[j] = data(ii++); |
|
3090 } |
|
3091 D[nc-1] = data(ii); |
|
3092 } |
|
3093 else |
|
3094 { |
|
3095 D[0] = 0.; |
|
3096 for (int i = 0; i < nr - 1; i++) |
|
3097 { |
|
3098 D[i+1] = 0.; |
|
3099 DL[i] = 0.; |
|
3100 DU[i] = 0.; |
|
3101 } |
|
3102 |
|
3103 for (int j = 0; j < nc; j++) |
|
3104 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3105 { |
|
3106 if (ridx(i) == j) |
|
3107 D[j] = data(i); |
|
3108 else if (ridx(i) == j + 1) |
|
3109 DL[j] = data(i); |
|
3110 else if (ridx(i) == j - 1) |
|
3111 DU[j] = data(i); |
|
3112 } |
|
3113 } |
|
3114 |
|
3115 int b_nc = b.cols(); |
|
3116 retval = b; |
|
3117 double *result = retval.fortran_vec (); |
|
3118 |
|
3119 F77_XFCN (dgtsv, DGTSV, (nr, b_nc, DL, D, DU, result, |
|
3120 b.rows(), err)); |
|
3121 |
|
3122 if (f77_exception_encountered) |
|
3123 (*current_liboctave_error_handler) |
|
3124 ("unrecoverable error in dgtsv"); |
|
3125 else if (err != 0) |
|
3126 { |
|
3127 rcond = 0.; |
|
3128 err = -2; |
|
3129 |
|
3130 if (sing_handler) |
|
3131 sing_handler (rcond); |
|
3132 else |
|
3133 (*current_liboctave_error_handler) |
|
3134 ("matrix singular to machine precision"); |
|
3135 |
|
3136 } |
|
3137 else |
|
3138 rcond = 1.; |
|
3139 } |
|
3140 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3141 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3142 } |
|
3143 |
|
3144 return retval; |
|
3145 } |
|
3146 |
|
3147 SparseMatrix |
|
3148 SparseMatrix::trisolve (SparseType &mattype, const SparseMatrix& b, int& err, |
|
3149 double& rcond, solve_singularity_handler sing_handler) const |
|
3150 { |
|
3151 SparseMatrix retval; |
|
3152 |
|
3153 int nr = rows (); |
|
3154 int nc = cols (); |
|
3155 err = 0; |
|
3156 |
|
3157 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3158 (*current_liboctave_error_handler) |
|
3159 ("matrix dimension mismatch solution of linear equations"); |
|
3160 else |
|
3161 { |
|
3162 // Print spparms("spumoni") info if requested |
|
3163 int typ = mattype.type (); |
|
3164 mattype.info (); |
|
3165 |
|
3166 // Note can't treat symmetric case as there is no dpttrf function |
|
3167 if (typ == SparseType::Tridiagonal || |
|
3168 typ == SparseType::Tridiagonal_Hermitian) |
|
3169 { |
|
3170 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2); |
|
3171 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1); |
|
3172 OCTAVE_LOCAL_BUFFER (double, D, nr); |
|
3173 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1); |
|
3174 Array<int> ipvt (nr); |
|
3175 int *pipvt = ipvt.fortran_vec (); |
|
3176 |
|
3177 if (mattype.is_dense ()) |
|
3178 { |
|
3179 int ii = 0; |
|
3180 |
|
3181 for (int j = 0; j < nc-1; j++) |
|
3182 { |
|
3183 D[j] = data(ii++); |
|
3184 DL[j] = data(ii++); |
|
3185 DU[j] = data(ii++); |
|
3186 } |
|
3187 D[nc-1] = data(ii); |
|
3188 } |
|
3189 else |
|
3190 { |
|
3191 D[0] = 0.; |
|
3192 for (int i = 0; i < nr - 1; i++) |
|
3193 { |
|
3194 D[i+1] = 0.; |
|
3195 DL[i] = 0.; |
|
3196 DU[i] = 0.; |
|
3197 } |
|
3198 |
|
3199 for (int j = 0; j < nc; j++) |
|
3200 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3201 { |
|
3202 if (ridx(i) == j) |
|
3203 D[j] = data(i); |
|
3204 else if (ridx(i) == j + 1) |
|
3205 DL[j] = data(i); |
|
3206 else if (ridx(i) == j - 1) |
|
3207 DU[j] = data(i); |
|
3208 } |
|
3209 } |
|
3210 |
|
3211 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); |
|
3212 |
|
3213 if (f77_exception_encountered) |
|
3214 (*current_liboctave_error_handler) |
|
3215 ("unrecoverable error in dgttrf"); |
|
3216 else |
|
3217 { |
|
3218 rcond = 0.0; |
|
3219 if (err != 0) |
|
3220 { |
|
3221 err = -2; |
|
3222 |
|
3223 if (sing_handler) |
|
3224 sing_handler (rcond); |
|
3225 else |
|
3226 (*current_liboctave_error_handler) |
|
3227 ("matrix singular to machine precision"); |
|
3228 |
|
3229 } |
|
3230 else |
|
3231 { |
|
3232 char job = 'N'; |
|
3233 volatile int x_nz = b.nnz (); |
|
3234 int b_nc = b.cols (); |
|
3235 retval = SparseMatrix (nr, b_nc, x_nz); |
|
3236 retval.xcidx(0) = 0; |
|
3237 volatile int ii = 0; |
|
3238 |
|
3239 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
3240 |
|
3241 for (volatile int j = 0; j < b_nc; j++) |
|
3242 { |
|
3243 for (int i = 0; i < nr; i++) |
|
3244 work[i] = 0.; |
|
3245 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
3246 work[b.ridx(i)] = b.data(i); |
|
3247 |
|
3248 F77_XFCN (dgttrs, DGTTRS, |
|
3249 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3250 nr, 1, DL, D, DU, DU2, pipvt, |
|
3251 work, b.rows (), err |
|
3252 F77_CHAR_ARG_LEN (1))); |
|
3253 |
|
3254 if (f77_exception_encountered) |
|
3255 { |
|
3256 (*current_liboctave_error_handler) |
|
3257 ("unrecoverable error in dgttrs"); |
|
3258 break; |
|
3259 } |
|
3260 |
|
3261 // Count non-zeros in work vector and adjust |
|
3262 // space in retval if needed |
|
3263 int new_nnz = 0; |
|
3264 for (int i = 0; i < nr; i++) |
|
3265 if (work[i] != 0.) |
|
3266 new_nnz++; |
|
3267 |
|
3268 if (ii + new_nnz > x_nz) |
|
3269 { |
|
3270 // Resize the sparse matrix |
|
3271 int sz = new_nnz * (b_nc - j) + x_nz; |
|
3272 retval.change_capacity (sz); |
|
3273 x_nz = sz; |
|
3274 } |
|
3275 |
|
3276 for (int i = 0; i < nr; i++) |
|
3277 if (work[i] != 0.) |
|
3278 { |
|
3279 retval.xridx(ii) = i; |
|
3280 retval.xdata(ii++) = work[i]; |
|
3281 } |
|
3282 retval.xcidx(j+1) = ii; |
|
3283 } |
|
3284 |
|
3285 retval.maybe_compress (); |
|
3286 } |
|
3287 } |
|
3288 } |
|
3289 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3290 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3291 } |
|
3292 |
|
3293 return retval; |
|
3294 } |
|
3295 |
|
3296 ComplexMatrix |
|
3297 SparseMatrix::trisolve (SparseType &mattype, const ComplexMatrix& b, int& err, |
|
3298 double& rcond, solve_singularity_handler sing_handler) const |
|
3299 { |
|
3300 ComplexMatrix retval; |
|
3301 |
|
3302 int nr = rows (); |
|
3303 int nc = cols (); |
|
3304 err = 0; |
|
3305 |
|
3306 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3307 (*current_liboctave_error_handler) |
|
3308 ("matrix dimension mismatch solution of linear equations"); |
|
3309 else |
|
3310 { |
|
3311 // Print spparms("spumoni") info if requested |
|
3312 volatile int typ = mattype.type (); |
|
3313 mattype.info (); |
|
3314 |
|
3315 // Note can't treat symmetric case as there is no dpttrf function |
|
3316 if (typ == SparseType::Tridiagonal_Hermitian) |
|
3317 { |
|
3318 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3319 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3320 |
|
3321 if (mattype.is_dense ()) |
|
3322 { |
|
3323 int ii = 0; |
|
3324 |
|
3325 for (int j = 0; j < nc-1; j++) |
|
3326 { |
|
3327 D[j] = data(ii++); |
|
3328 DL[j] = data(ii); |
|
3329 ii += 2; |
|
3330 } |
|
3331 D[nc-1] = data(ii); |
|
3332 } |
|
3333 else |
|
3334 { |
|
3335 D[0] = 0.; |
|
3336 for (int i = 0; i < nr - 1; i++) |
|
3337 { |
|
3338 D[i+1] = 0.; |
|
3339 DL[i] = 0.; |
|
3340 } |
|
3341 |
|
3342 for (int j = 0; j < nc; j++) |
|
3343 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3344 { |
|
3345 if (ridx(i) == j) |
|
3346 D[j] = data(i); |
|
3347 else if (ridx(i) == j + 1) |
|
3348 DL[j] = data(i); |
|
3349 } |
|
3350 } |
|
3351 |
|
3352 int b_nr = b.rows (); |
|
3353 int b_nc = b.cols(); |
|
3354 rcond = 1.; |
|
3355 |
|
3356 retval = b; |
|
3357 Complex *result = retval.fortran_vec (); |
|
3358 |
|
3359 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result, |
|
3360 b_nr, err)); |
|
3361 |
|
3362 if (f77_exception_encountered) |
|
3363 { |
|
3364 (*current_liboctave_error_handler) |
|
3365 ("unrecoverable error in zptsv"); |
|
3366 err = -1; |
|
3367 } |
|
3368 else if (err != 0) |
|
3369 { |
|
3370 err = 0; |
|
3371 mattype.mark_as_unsymmetric (); |
|
3372 typ = SparseType::Tridiagonal; |
|
3373 } |
|
3374 } |
|
3375 |
|
3376 if (typ == SparseType::Tridiagonal) |
|
3377 { |
|
3378 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); |
|
3379 OCTAVE_LOCAL_BUFFER (Complex, D, nr); |
|
3380 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); |
|
3381 |
|
3382 if (mattype.is_dense ()) |
|
3383 { |
|
3384 int ii = 0; |
|
3385 |
|
3386 for (int j = 0; j < nc-1; j++) |
|
3387 { |
|
3388 D[j] = data(ii++); |
|
3389 DL[j] = data(ii++); |
|
3390 DU[j] = data(ii++); |
|
3391 } |
|
3392 D[nc-1] = data(ii); |
|
3393 } |
|
3394 else |
|
3395 { |
|
3396 D[0] = 0.; |
|
3397 for (int i = 0; i < nr - 1; i++) |
|
3398 { |
|
3399 D[i+1] = 0.; |
|
3400 DL[i] = 0.; |
|
3401 DU[i] = 0.; |
|
3402 } |
|
3403 |
|
3404 for (int j = 0; j < nc; j++) |
|
3405 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3406 { |
|
3407 if (ridx(i) == j) |
|
3408 D[j] = data(i); |
|
3409 else if (ridx(i) == j + 1) |
|
3410 DL[j] = data(i); |
|
3411 else if (ridx(i) == j - 1) |
|
3412 DU[j] = data(i); |
|
3413 } |
|
3414 } |
|
3415 |
|
3416 int b_nr = b.rows(); |
|
3417 int b_nc = b.cols(); |
|
3418 rcond = 1.; |
|
3419 |
|
3420 retval = b; |
|
3421 Complex *result = retval.fortran_vec (); |
|
3422 |
|
3423 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result, |
|
3424 b_nr, err)); |
|
3425 |
|
3426 if (f77_exception_encountered) |
|
3427 { |
|
3428 (*current_liboctave_error_handler) |
|
3429 ("unrecoverable error in zgtsv"); |
|
3430 err = -1; |
|
3431 } |
|
3432 else if (err != 0) |
|
3433 { |
|
3434 rcond = 0.; |
|
3435 err = -2; |
|
3436 |
|
3437 if (sing_handler) |
|
3438 sing_handler (rcond); |
|
3439 else |
|
3440 (*current_liboctave_error_handler) |
|
3441 ("matrix singular to machine precision"); |
|
3442 } |
|
3443 } |
|
3444 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3445 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3446 } |
|
3447 |
|
3448 return retval; |
|
3449 } |
|
3450 |
|
3451 SparseComplexMatrix |
|
3452 SparseMatrix::trisolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
3453 int& err, double& rcond, |
|
3454 solve_singularity_handler sing_handler) const |
|
3455 { |
|
3456 SparseComplexMatrix retval; |
|
3457 |
|
3458 int nr = rows (); |
|
3459 int nc = cols (); |
|
3460 err = 0; |
|
3461 |
|
3462 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3463 (*current_liboctave_error_handler) |
|
3464 ("matrix dimension mismatch solution of linear equations"); |
|
3465 else |
|
3466 { |
|
3467 // Print spparms("spumoni") info if requested |
|
3468 int typ = mattype.type (); |
|
3469 mattype.info (); |
|
3470 |
|
3471 // Note can't treat symmetric case as there is no dpttrf function |
|
3472 if (typ == SparseType::Tridiagonal || |
|
3473 typ == SparseType::Tridiagonal_Hermitian) |
|
3474 { |
|
3475 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2); |
|
3476 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1); |
|
3477 OCTAVE_LOCAL_BUFFER (double, D, nr); |
|
3478 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1); |
|
3479 Array<int> ipvt (nr); |
|
3480 int *pipvt = ipvt.fortran_vec (); |
|
3481 |
|
3482 if (mattype.is_dense ()) |
|
3483 { |
|
3484 int ii = 0; |
|
3485 |
|
3486 for (int j = 0; j < nc-1; j++) |
|
3487 { |
|
3488 D[j] = data(ii++); |
|
3489 DL[j] = data(ii++); |
|
3490 DU[j] = data(ii++); |
|
3491 } |
|
3492 D[nc-1] = data(ii); |
|
3493 } |
|
3494 else |
|
3495 { |
|
3496 D[0] = 0.; |
|
3497 for (int i = 0; i < nr - 1; i++) |
|
3498 { |
|
3499 D[i+1] = 0.; |
|
3500 DL[i] = 0.; |
|
3501 DU[i] = 0.; |
|
3502 } |
|
3503 |
|
3504 for (int j = 0; j < nc; j++) |
|
3505 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3506 { |
|
3507 if (ridx(i) == j) |
|
3508 D[j] = data(i); |
|
3509 else if (ridx(i) == j + 1) |
|
3510 DL[j] = data(i); |
|
3511 else if (ridx(i) == j - 1) |
|
3512 DU[j] = data(i); |
|
3513 } |
|
3514 } |
|
3515 |
|
3516 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err)); |
|
3517 |
|
3518 if (f77_exception_encountered) |
|
3519 (*current_liboctave_error_handler) |
|
3520 ("unrecoverable error in dgttrf"); |
|
3521 else |
|
3522 { |
|
3523 rcond = 0.0; |
|
3524 if (err != 0) |
|
3525 { |
|
3526 err = -2; |
|
3527 |
|
3528 if (sing_handler) |
|
3529 sing_handler (rcond); |
|
3530 else |
|
3531 (*current_liboctave_error_handler) |
|
3532 ("matrix singular to machine precision"); |
|
3533 } |
|
3534 else |
|
3535 { |
|
3536 rcond = 1.; |
|
3537 char job = 'N'; |
|
3538 int b_nr = b.rows (); |
|
3539 int b_nc = b.cols (); |
|
3540 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
3541 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
3542 |
|
3543 // Take a first guess that the number of non-zero terms |
|
3544 // will be as many as in b |
|
3545 volatile int x_nz = b.nnz (); |
|
3546 volatile int ii = 0; |
|
3547 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
3548 |
|
3549 retval.xcidx(0) = 0; |
|
3550 for (volatile int j = 0; j < b_nc; j++) |
|
3551 { |
|
3552 |
|
3553 for (int i = 0; i < b_nr; i++) |
|
3554 { |
|
3555 Complex c = b (i,j); |
|
3556 Bx[i] = ::real (c); |
|
3557 Bz[i] = ::imag (c); |
|
3558 } |
|
3559 |
|
3560 |
|
3561 F77_XFCN (dgttrs, DGTTRS, |
|
3562 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3563 nr, 1, DL, D, DU, DU2, pipvt, |
|
3564 Bx, b_nr, err |
|
3565 F77_CHAR_ARG_LEN (1))); |
|
3566 |
|
3567 if (f77_exception_encountered) |
|
3568 { |
|
3569 (*current_liboctave_error_handler) |
|
3570 ("unrecoverable error in dgttrs"); |
|
3571 break; |
|
3572 } |
|
3573 |
|
3574 if (err != 0) |
|
3575 { |
|
3576 (*current_liboctave_error_handler) |
|
3577 ("SparseMatrix::solve solve failed"); |
|
3578 |
|
3579 err = -1; |
|
3580 break; |
|
3581 } |
|
3582 |
|
3583 F77_XFCN (dgttrs, DGTTRS, |
|
3584 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3585 nr, 1, DL, D, DU, DU2, pipvt, |
|
3586 Bz, b_nr, err |
|
3587 F77_CHAR_ARG_LEN (1))); |
|
3588 |
|
3589 if (f77_exception_encountered) |
|
3590 { |
|
3591 (*current_liboctave_error_handler) |
|
3592 ("unrecoverable error in dgttrs"); |
|
3593 break; |
|
3594 } |
|
3595 |
|
3596 if (err != 0) |
|
3597 { |
|
3598 (*current_liboctave_error_handler) |
|
3599 ("SparseMatrix::solve solve failed"); |
|
3600 |
|
3601 err = -1; |
|
3602 break; |
|
3603 } |
|
3604 |
|
3605 // Count non-zeros in work vector and adjust |
|
3606 // space in retval if needed |
|
3607 int new_nnz = 0; |
|
3608 for (int i = 0; i < nr; i++) |
|
3609 if (Bx[i] != 0. || Bz[i] != 0.) |
|
3610 new_nnz++; |
|
3611 |
|
3612 if (ii + new_nnz > x_nz) |
|
3613 { |
|
3614 // Resize the sparse matrix |
|
3615 int sz = new_nnz * (b_nc - j) + x_nz; |
|
3616 retval.change_capacity (sz); |
|
3617 x_nz = sz; |
|
3618 } |
|
3619 |
|
3620 for (int i = 0; i < nr; i++) |
|
3621 if (Bx[i] != 0. || Bz[i] != 0.) |
|
3622 { |
|
3623 retval.xridx(ii) = i; |
|
3624 retval.xdata(ii++) = |
|
3625 Complex (Bx[i], Bz[i]); |
|
3626 } |
|
3627 |
|
3628 retval.xcidx(j+1) = ii; |
|
3629 } |
|
3630 |
|
3631 retval.maybe_compress (); |
|
3632 } |
|
3633 } |
|
3634 } |
|
3635 else if (typ != SparseType::Tridiagonal_Hermitian) |
|
3636 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3637 } |
|
3638 |
|
3639 return retval; |
|
3640 } |
|
3641 |
|
3642 Matrix |
|
3643 SparseMatrix::bsolve (SparseType &mattype, const Matrix& b, int& err, |
|
3644 double& rcond, |
|
3645 solve_singularity_handler sing_handler) const |
|
3646 { |
|
3647 Matrix retval; |
|
3648 |
|
3649 int nr = rows (); |
|
3650 int nc = cols (); |
|
3651 err = 0; |
|
3652 |
|
3653 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3654 (*current_liboctave_error_handler) |
|
3655 ("matrix dimension mismatch solution of linear equations"); |
|
3656 else |
|
3657 { |
|
3658 // Print spparms("spumoni") info if requested |
|
3659 volatile int typ = mattype.type (); |
|
3660 mattype.info (); |
|
3661 |
|
3662 if (typ == SparseType::Banded_Hermitian) |
|
3663 { |
|
3664 int n_lower = mattype.nlower (); |
|
3665 int ldm = n_lower + 1; |
|
3666 Matrix m_band (ldm, nc); |
|
3667 double *tmp_data = m_band.fortran_vec (); |
|
3668 |
|
3669 if (! mattype.is_dense ()) |
|
3670 { |
|
3671 int ii = 0; |
|
3672 |
|
3673 for (int j = 0; j < ldm; j++) |
|
3674 for (int i = 0; i < nc; i++) |
|
3675 tmp_data[ii++] = 0.; |
|
3676 } |
|
3677 |
|
3678 for (int j = 0; j < nc; j++) |
|
3679 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3680 { |
|
3681 int ri = ridx (i); |
|
3682 if (ri >= j) |
|
3683 m_band(ri - j, j) = data(i); |
|
3684 } |
|
3685 |
|
3686 // Calculate the norm of the matrix, for later use. |
|
3687 // double anorm = m_band.abs().sum().row(0).max(); |
|
3688 |
|
3689 char job = 'L'; |
|
3690 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3691 nr, n_lower, tmp_data, ldm, err |
|
3692 F77_CHAR_ARG_LEN (1))); |
|
3693 |
|
3694 if (f77_exception_encountered) |
|
3695 (*current_liboctave_error_handler) |
|
3696 ("unrecoverable error in dpbtrf"); |
|
3697 else |
|
3698 { |
|
3699 rcond = 0.0; |
|
3700 if (err != 0) |
|
3701 { |
|
3702 // Matrix is not positive definite!! Fall through to |
|
3703 // unsymmetric banded solver. |
|
3704 mattype.mark_as_unsymmetric (); |
|
3705 typ = SparseType::Banded; |
|
3706 err = 0; |
|
3707 } |
|
3708 else |
|
3709 { |
|
3710 // Unfortunately, the time to calculate the condition |
|
3711 // number is dominant for narrow banded matrices and |
|
3712 // so we rely on the "err" flag from xPBTRF to flag |
|
3713 // singularity. The commented code below is left here |
|
3714 // for reference |
|
3715 |
|
3716 //Array<double> z (3 * nr); |
|
3717 //double *pz = z.fortran_vec (); |
|
3718 //Array<int> iz (nr); |
|
3719 //int *piz = iz.fortran_vec (); |
|
3720 // |
|
3721 //F77_XFCN (dpbcon, DGBCON, |
|
3722 // (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3723 // nr, n_lower, tmp_data, ldm, |
|
3724 // anorm, rcond, pz, piz, err |
|
3725 // F77_CHAR_ARG_LEN (1))); |
|
3726 // |
|
3727 // |
|
3728 //if (f77_exception_encountered) |
|
3729 // (*current_liboctave_error_handler) |
|
3730 // ("unrecoverable error in dpbcon"); |
|
3731 // |
|
3732 //if (err != 0) |
|
3733 // err = -2; |
|
3734 // |
|
3735 //volatile double rcond_plus_one = rcond + 1.0; |
|
3736 // |
|
3737 //if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
3738 // { |
|
3739 // err = -2; |
|
3740 // |
|
3741 // if (sing_handler) |
|
3742 // sing_handler (rcond); |
|
3743 // else |
|
3744 // (*current_liboctave_error_handler) |
|
3745 // ("matrix singular to machine precision, rcond = %g", |
|
3746 // rcond); |
|
3747 // } |
|
3748 //else |
|
3749 // REST OF CODE, EXCEPT rcond=1 |
|
3750 |
|
3751 rcond = 1.; |
|
3752 retval = b; |
|
3753 double *result = retval.fortran_vec (); |
|
3754 |
|
3755 int b_nc = b.cols (); |
|
3756 |
|
3757 F77_XFCN (dpbtrs, DPBTRS, |
|
3758 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3759 nr, n_lower, b_nc, tmp_data, |
|
3760 ldm, result, b.rows(), err |
|
3761 F77_CHAR_ARG_LEN (1))); |
|
3762 |
|
3763 if (f77_exception_encountered) |
|
3764 (*current_liboctave_error_handler) |
|
3765 ("unrecoverable error in dpbtrs"); |
|
3766 |
|
3767 if (err != 0) |
|
3768 { |
|
3769 (*current_liboctave_error_handler) |
|
3770 ("SparseMatrix::solve solve failed"); |
|
3771 err = -1; |
|
3772 } |
|
3773 } |
|
3774 } |
|
3775 } |
|
3776 |
|
3777 if (typ == SparseType::Banded) |
|
3778 { |
|
3779 // Create the storage for the banded form of the sparse matrix |
|
3780 int n_upper = mattype.nupper (); |
|
3781 int n_lower = mattype.nlower (); |
|
3782 int ldm = n_upper + 2 * n_lower + 1; |
|
3783 |
|
3784 Matrix m_band (ldm, nc); |
|
3785 double *tmp_data = m_band.fortran_vec (); |
|
3786 |
|
3787 if (! mattype.is_dense ()) |
|
3788 { |
|
3789 int ii = 0; |
|
3790 |
|
3791 for (int j = 0; j < ldm; j++) |
|
3792 for (int i = 0; i < nc; i++) |
|
3793 tmp_data[ii++] = 0.; |
|
3794 } |
|
3795 |
|
3796 for (int j = 0; j < nc; j++) |
|
3797 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3798 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
3799 |
|
3800 Array<int> ipvt (nr); |
|
3801 int *pipvt = ipvt.fortran_vec (); |
|
3802 |
|
3803 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
3804 ldm, pipvt, err)); |
|
3805 |
|
3806 if (f77_exception_encountered) |
|
3807 (*current_liboctave_error_handler) |
|
3808 ("unrecoverable error in dgbtrf"); |
|
3809 else |
|
3810 { |
|
3811 // Throw-away extra info LAPACK gives so as to not |
|
3812 // change output. |
|
3813 rcond = 0.0; |
|
3814 if (err != 0) |
|
3815 { |
|
3816 err = -2; |
|
3817 |
|
3818 if (sing_handler) |
|
3819 sing_handler (rcond); |
|
3820 else |
|
3821 (*current_liboctave_error_handler) |
|
3822 ("matrix singular to machine precision"); |
|
3823 |
|
3824 } |
|
3825 else |
|
3826 { |
|
3827 char job = '1'; |
|
3828 |
|
3829 // Unfortunately, the time to calculate the condition |
|
3830 // number is dominant for narrow banded matrices and |
|
3831 // so we rely on the "err" flag from xPBTRF to flag |
|
3832 // singularity. The commented code below is left here |
|
3833 // for reference |
|
3834 |
|
3835 //F77_XFCN (dgbcon, DGBCON, |
|
3836 // (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3837 // nc, n_lower, n_upper, tmp_data, ldm, pipvt, |
|
3838 // anorm, rcond, pz, piz, err |
|
3839 // F77_CHAR_ARG_LEN (1))); |
|
3840 // |
|
3841 //if (f77_exception_encountered) |
|
3842 // (*current_liboctave_error_handler) |
|
3843 // ("unrecoverable error in dgbcon"); |
|
3844 // |
|
3845 // if (err != 0) |
|
3846 // err = -2; |
|
3847 // |
|
3848 //volatile double rcond_plus_one = rcond + 1.0; |
|
3849 // |
|
3850 //if (rcond_plus_one == 1.0 || xisnan (rcond)) |
|
3851 // { |
|
3852 // err = -2; |
|
3853 // |
|
3854 // if (sing_handler) |
|
3855 // sing_handler (rcond); |
|
3856 // else |
|
3857 // (*current_liboctave_error_handler) |
|
3858 // ("matrix singular to machine precision, rcond = %g", |
|
3859 // rcond); |
|
3860 // } |
|
3861 //else |
|
3862 // REST OF CODE, EXCEPT rcond=1 |
|
3863 |
|
3864 rcond = 1.; |
|
3865 retval = b; |
|
3866 double *result = retval.fortran_vec (); |
|
3867 |
|
3868 int b_nc = b.cols (); |
|
3869 |
|
3870 job = 'N'; |
|
3871 F77_XFCN (dgbtrs, DGBTRS, |
|
3872 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3873 nr, n_lower, n_upper, b_nc, tmp_data, |
|
3874 ldm, pipvt, result, b.rows(), err |
|
3875 F77_CHAR_ARG_LEN (1))); |
|
3876 |
|
3877 if (f77_exception_encountered) |
|
3878 (*current_liboctave_error_handler) |
|
3879 ("unrecoverable error in dgbtrs"); |
|
3880 } |
|
3881 } |
|
3882 } |
|
3883 else if (typ != SparseType::Banded_Hermitian) |
|
3884 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
3885 } |
|
3886 |
|
3887 return retval; |
|
3888 } |
|
3889 |
|
3890 SparseMatrix |
|
3891 SparseMatrix::bsolve (SparseType &mattype, const SparseMatrix& b, int& err, |
|
3892 double& rcond, solve_singularity_handler sing_handler) const |
|
3893 { |
|
3894 SparseMatrix retval; |
|
3895 |
|
3896 int nr = rows (); |
|
3897 int nc = cols (); |
|
3898 err = 0; |
|
3899 |
|
3900 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
3901 (*current_liboctave_error_handler) |
|
3902 ("matrix dimension mismatch solution of linear equations"); |
|
3903 else |
|
3904 { |
|
3905 // Print spparms("spumoni") info if requested |
|
3906 volatile int typ = mattype.type (); |
|
3907 mattype.info (); |
|
3908 |
|
3909 if (typ == SparseType::Banded_Hermitian) |
|
3910 { |
|
3911 int n_lower = mattype.nlower (); |
|
3912 int ldm = n_lower + 1; |
|
3913 |
|
3914 Matrix m_band (ldm, nc); |
|
3915 double *tmp_data = m_band.fortran_vec (); |
|
3916 |
|
3917 if (! mattype.is_dense ()) |
|
3918 { |
|
3919 int ii = 0; |
|
3920 |
|
3921 for (int j = 0; j < ldm; j++) |
|
3922 for (int i = 0; i < nc; i++) |
|
3923 tmp_data[ii++] = 0.; |
|
3924 } |
|
3925 |
|
3926 for (int j = 0; j < nc; j++) |
|
3927 for (int i = cidx(j); i < cidx(j+1); i++) |
|
3928 { |
|
3929 int ri = ridx (i); |
|
3930 if (ri >= j) |
|
3931 m_band(ri - j, j) = data(i); |
|
3932 } |
|
3933 |
|
3934 char job = 'L'; |
|
3935 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3936 nr, n_lower, tmp_data, ldm, err |
|
3937 F77_CHAR_ARG_LEN (1))); |
|
3938 |
|
3939 if (f77_exception_encountered) |
|
3940 (*current_liboctave_error_handler) |
|
3941 ("unrecoverable error in dpbtrf"); |
|
3942 else |
|
3943 { |
|
3944 rcond = 0.0; |
|
3945 if (err != 0) |
|
3946 { |
|
3947 mattype.mark_as_unsymmetric (); |
|
3948 typ = SparseType::Banded; |
|
3949 err = 0; |
|
3950 } |
|
3951 else |
|
3952 { |
|
3953 rcond = 1.; |
|
3954 int b_nr = b.rows (); |
|
3955 int b_nc = b.cols (); |
|
3956 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
3957 |
|
3958 // Take a first guess that the number of non-zero terms |
|
3959 // will be as many as in b |
|
3960 volatile int x_nz = b.nnz (); |
|
3961 volatile int ii = 0; |
|
3962 retval = SparseMatrix (b_nr, b_nc, x_nz); |
|
3963 |
|
3964 retval.xcidx(0) = 0; |
|
3965 for (volatile int j = 0; j < b_nc; j++) |
|
3966 { |
|
3967 for (int i = 0; i < b_nr; i++) |
|
3968 Bx[i] = b.elem (i, j); |
|
3969 |
|
3970 F77_XFCN (dpbtrs, DPBTRS, |
|
3971 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
3972 nr, n_lower, 1, tmp_data, |
|
3973 ldm, Bx, b_nr, err |
|
3974 F77_CHAR_ARG_LEN (1))); |
|
3975 |
|
3976 if (f77_exception_encountered) |
|
3977 { |
|
3978 (*current_liboctave_error_handler) |
|
3979 ("unrecoverable error in dpbtrs"); |
|
3980 err = -1; |
|
3981 break; |
|
3982 } |
|
3983 |
|
3984 if (err != 0) |
|
3985 { |
|
3986 (*current_liboctave_error_handler) |
|
3987 ("SparseMatrix::solve solve failed"); |
|
3988 err = -1; |
|
3989 break; |
|
3990 } |
|
3991 |
|
3992 for (int i = 0; i < b_nr; i++) |
|
3993 { |
|
3994 double tmp = Bx[i]; |
|
3995 if (tmp != 0.0) |
|
3996 { |
|
3997 if (ii == x_nz) |
|
3998 { |
|
3999 // Resize the sparse matrix |
|
4000 int sz = x_nz * (b_nc - j) / b_nc; |
|
4001 sz = (sz > 10 ? sz : 10) + x_nz; |
|
4002 retval.change_capacity (sz); |
|
4003 x_nz = sz; |
|
4004 } |
|
4005 retval.xdata(ii) = tmp; |
|
4006 retval.xridx(ii++) = i; |
|
4007 } |
|
4008 } |
|
4009 retval.xcidx(j+1) = ii; |
|
4010 } |
|
4011 |
|
4012 retval.maybe_compress (); |
|
4013 } |
|
4014 } |
|
4015 } |
|
4016 |
|
4017 if (typ == SparseType::Banded) |
|
4018 { |
|
4019 // Create the storage for the banded form of the sparse matrix |
|
4020 int n_upper = mattype.nupper (); |
|
4021 int n_lower = mattype.nlower (); |
|
4022 int ldm = n_upper + 2 * n_lower + 1; |
|
4023 |
|
4024 Matrix m_band (ldm, nc); |
|
4025 double *tmp_data = m_band.fortran_vec (); |
|
4026 |
|
4027 if (! mattype.is_dense ()) |
|
4028 { |
|
4029 int ii = 0; |
|
4030 |
|
4031 for (int j = 0; j < ldm; j++) |
|
4032 for (int i = 0; i < nc; i++) |
|
4033 tmp_data[ii++] = 0.; |
|
4034 } |
|
4035 |
|
4036 for (int j = 0; j < nc; j++) |
|
4037 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4038 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
4039 |
|
4040 Array<int> ipvt (nr); |
|
4041 int *pipvt = ipvt.fortran_vec (); |
|
4042 |
|
4043 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
4044 ldm, pipvt, err)); |
|
4045 |
|
4046 if (f77_exception_encountered) |
|
4047 (*current_liboctave_error_handler) |
|
4048 ("unrecoverable error in dgbtrf"); |
|
4049 else |
|
4050 { |
|
4051 rcond = 0.0; |
|
4052 if (err != 0) |
|
4053 { |
|
4054 err = -2; |
|
4055 |
|
4056 if (sing_handler) |
|
4057 sing_handler (rcond); |
|
4058 else |
|
4059 (*current_liboctave_error_handler) |
|
4060 ("matrix singular to machine precision"); |
|
4061 |
|
4062 } |
|
4063 else |
|
4064 { |
|
4065 char job = 'N'; |
|
4066 volatile int x_nz = b.nnz (); |
|
4067 int b_nc = b.cols (); |
|
4068 retval = SparseMatrix (nr, b_nc, x_nz); |
|
4069 retval.xcidx(0) = 0; |
|
4070 volatile int ii = 0; |
|
4071 |
|
4072 OCTAVE_LOCAL_BUFFER (double, work, nr); |
|
4073 |
|
4074 for (volatile int j = 0; j < b_nc; j++) |
|
4075 { |
|
4076 for (int i = 0; i < nr; i++) |
|
4077 work[i] = 0.; |
|
4078 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
4079 work[b.ridx(i)] = b.data(i); |
|
4080 |
|
4081 F77_XFCN (dgbtrs, DGBTRS, |
|
4082 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4083 nr, n_lower, n_upper, 1, tmp_data, |
|
4084 ldm, pipvt, work, b.rows (), err |
|
4085 F77_CHAR_ARG_LEN (1))); |
|
4086 |
|
4087 if (f77_exception_encountered) |
|
4088 { |
|
4089 (*current_liboctave_error_handler) |
|
4090 ("unrecoverable error in dgbtrs"); |
|
4091 break; |
|
4092 } |
|
4093 |
|
4094 // Count non-zeros in work vector and adjust |
|
4095 // space in retval if needed |
|
4096 int new_nnz = 0; |
|
4097 for (int i = 0; i < nr; i++) |
|
4098 if (work[i] != 0.) |
|
4099 new_nnz++; |
|
4100 |
|
4101 if (ii + new_nnz > x_nz) |
|
4102 { |
|
4103 // Resize the sparse matrix |
|
4104 int sz = new_nnz * (b_nc - j) + x_nz; |
|
4105 retval.change_capacity (sz); |
|
4106 x_nz = sz; |
|
4107 } |
|
4108 |
|
4109 for (int i = 0; i < nr; i++) |
|
4110 if (work[i] != 0.) |
|
4111 { |
|
4112 retval.xridx(ii) = i; |
|
4113 retval.xdata(ii++) = work[i]; |
|
4114 } |
|
4115 retval.xcidx(j+1) = ii; |
|
4116 } |
|
4117 |
|
4118 retval.maybe_compress (); |
|
4119 } |
|
4120 } |
|
4121 } |
|
4122 else if (typ != SparseType::Banded_Hermitian) |
|
4123 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4124 } |
|
4125 |
|
4126 return retval; |
|
4127 } |
|
4128 |
|
4129 ComplexMatrix |
|
4130 SparseMatrix::bsolve (SparseType &mattype, const ComplexMatrix& b, int& err, |
|
4131 double& rcond, solve_singularity_handler sing_handler) const |
|
4132 { |
|
4133 ComplexMatrix retval; |
|
4134 |
|
4135 int nr = rows (); |
|
4136 int nc = cols (); |
|
4137 err = 0; |
|
4138 |
|
4139 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4140 (*current_liboctave_error_handler) |
|
4141 ("matrix dimension mismatch solution of linear equations"); |
|
4142 else |
|
4143 { |
|
4144 // Print spparms("spumoni") info if requested |
|
4145 volatile int typ = mattype.type (); |
|
4146 mattype.info (); |
|
4147 |
|
4148 if (typ == SparseType::Banded_Hermitian) |
|
4149 { |
|
4150 int n_lower = mattype.nlower (); |
|
4151 int ldm = n_lower + 1; |
|
4152 |
|
4153 Matrix m_band (ldm, nc); |
|
4154 double *tmp_data = m_band.fortran_vec (); |
|
4155 |
|
4156 if (! mattype.is_dense ()) |
|
4157 { |
|
4158 int ii = 0; |
|
4159 |
|
4160 for (int j = 0; j < ldm; j++) |
|
4161 for (int i = 0; i < nc; i++) |
|
4162 tmp_data[ii++] = 0.; |
|
4163 } |
|
4164 |
|
4165 for (int j = 0; j < nc; j++) |
|
4166 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4167 { |
|
4168 int ri = ridx (i); |
|
4169 if (ri >= j) |
|
4170 m_band(ri - j, j) = data(i); |
|
4171 } |
|
4172 |
|
4173 char job = 'L'; |
|
4174 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4175 nr, n_lower, tmp_data, ldm, err |
|
4176 F77_CHAR_ARG_LEN (1))); |
|
4177 |
|
4178 if (f77_exception_encountered) |
|
4179 (*current_liboctave_error_handler) |
|
4180 ("unrecoverable error in dpbtrf"); |
|
4181 else |
|
4182 { |
|
4183 rcond = 0.0; |
|
4184 if (err != 0) |
|
4185 { |
|
4186 // Matrix is not positive definite!! Fall through to |
|
4187 // unsymmetric banded solver. |
|
4188 mattype.mark_as_unsymmetric (); |
|
4189 typ = SparseType::Banded; |
|
4190 err = 0; |
|
4191 } |
|
4192 else |
|
4193 { |
|
4194 rcond = 1.; |
|
4195 int b_nr = b.rows (); |
|
4196 int b_nc = b.cols (); |
|
4197 |
|
4198 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
4199 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
4200 |
|
4201 retval.resize (b_nr, b_nc); |
|
4202 |
|
4203 for (volatile int j = 0; j < b_nc; j++) |
|
4204 { |
|
4205 for (int i = 0; i < b_nr; i++) |
|
4206 { |
|
4207 Complex c = b (i,j); |
|
4208 Bx[i] = ::real (c); |
|
4209 Bz[i] = ::imag (c); |
|
4210 } |
|
4211 |
|
4212 F77_XFCN (dpbtrs, DPBTRS, |
|
4213 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4214 nr, n_lower, 1, tmp_data, |
|
4215 ldm, Bx, b_nr, err |
|
4216 F77_CHAR_ARG_LEN (1))); |
|
4217 |
|
4218 if (f77_exception_encountered) |
|
4219 { |
|
4220 (*current_liboctave_error_handler) |
|
4221 ("unrecoverable error in dpbtrs"); |
|
4222 err = -1; |
|
4223 break; |
|
4224 } |
|
4225 |
|
4226 if (err != 0) |
|
4227 { |
|
4228 (*current_liboctave_error_handler) |
|
4229 ("SparseMatrix::solve solve failed"); |
|
4230 err = -1; |
|
4231 break; |
|
4232 } |
|
4233 |
|
4234 F77_XFCN (dpbtrs, DPBTRS, |
|
4235 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4236 nr, n_lower, 1, tmp_data, |
|
4237 ldm, Bz, b.rows(), err |
|
4238 F77_CHAR_ARG_LEN (1))); |
|
4239 |
|
4240 if (f77_exception_encountered) |
|
4241 { |
|
4242 (*current_liboctave_error_handler) |
|
4243 ("unrecoverable error in dpbtrs"); |
|
4244 err = -1; |
|
4245 break; |
|
4246 } |
|
4247 |
|
4248 if (err != 0) |
|
4249 { |
|
4250 (*current_liboctave_error_handler) |
|
4251 ("SparseMatrix::solve solve failed"); |
|
4252 err = -1; |
|
4253 break; |
|
4254 } |
|
4255 |
|
4256 for (int i = 0; i < b_nr; i++) |
|
4257 retval (i, j) = Complex (Bx[i], Bz[i]); |
|
4258 } |
|
4259 } |
|
4260 } |
|
4261 } |
|
4262 |
|
4263 if (typ == SparseType::Banded) |
|
4264 { |
|
4265 // Create the storage for the banded form of the sparse matrix |
|
4266 int n_upper = mattype.nupper (); |
|
4267 int n_lower = mattype.nlower (); |
|
4268 int ldm = n_upper + 2 * n_lower + 1; |
|
4269 |
|
4270 Matrix m_band (ldm, nc); |
|
4271 double *tmp_data = m_band.fortran_vec (); |
|
4272 |
|
4273 if (! mattype.is_dense ()) |
|
4274 { |
|
4275 int ii = 0; |
|
4276 |
|
4277 for (int j = 0; j < ldm; j++) |
|
4278 for (int i = 0; i < nc; i++) |
|
4279 tmp_data[ii++] = 0.; |
|
4280 } |
|
4281 |
|
4282 for (int j = 0; j < nc; j++) |
|
4283 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4284 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
4285 |
|
4286 Array<int> ipvt (nr); |
|
4287 int *pipvt = ipvt.fortran_vec (); |
|
4288 |
|
4289 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
4290 ldm, pipvt, err)); |
|
4291 |
|
4292 if (f77_exception_encountered) |
|
4293 (*current_liboctave_error_handler) |
|
4294 ("unrecoverable error in dgbtrf"); |
|
4295 else |
|
4296 { |
|
4297 rcond = 0.0; |
|
4298 if (err != 0) |
|
4299 { |
|
4300 err = -2; |
|
4301 |
|
4302 if (sing_handler) |
|
4303 sing_handler (rcond); |
|
4304 else |
|
4305 (*current_liboctave_error_handler) |
|
4306 ("matrix singular to machine precision"); |
|
4307 |
|
4308 } |
|
4309 else |
|
4310 { |
|
4311 char job = 'N'; |
|
4312 int b_nc = b.cols (); |
|
4313 retval.resize (nr,b_nc); |
|
4314 |
|
4315 OCTAVE_LOCAL_BUFFER (double, Bz, nr); |
|
4316 OCTAVE_LOCAL_BUFFER (double, Bx, nr); |
|
4317 |
|
4318 for (volatile int j = 0; j < b_nc; j++) |
|
4319 { |
|
4320 for (int i = 0; i < nr; i++) |
|
4321 { |
|
4322 Complex c = b (i, j); |
|
4323 Bx[i] = ::real (c); |
|
4324 Bz[i] = ::imag (c); |
|
4325 } |
|
4326 |
|
4327 F77_XFCN (dgbtrs, DGBTRS, |
|
4328 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4329 nr, n_lower, n_upper, 1, tmp_data, |
|
4330 ldm, pipvt, Bx, b.rows (), err |
|
4331 F77_CHAR_ARG_LEN (1))); |
|
4332 |
|
4333 if (f77_exception_encountered) |
|
4334 { |
|
4335 (*current_liboctave_error_handler) |
|
4336 ("unrecoverable error in dgbtrs"); |
|
4337 break; |
|
4338 } |
|
4339 |
|
4340 F77_XFCN (dgbtrs, DGBTRS, |
|
4341 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4342 nr, n_lower, n_upper, 1, tmp_data, |
|
4343 ldm, pipvt, Bz, b.rows (), err |
|
4344 F77_CHAR_ARG_LEN (1))); |
|
4345 |
|
4346 if (f77_exception_encountered) |
|
4347 { |
|
4348 (*current_liboctave_error_handler) |
|
4349 ("unrecoverable error in dgbtrs"); |
|
4350 break; |
|
4351 } |
|
4352 |
|
4353 for (int i = 0; i < nr; i++) |
|
4354 retval (i, j) = Complex (Bx[i], Bz[i]); |
|
4355 } |
|
4356 } |
|
4357 } |
|
4358 } |
|
4359 else if (typ != SparseType::Banded_Hermitian) |
|
4360 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4361 } |
|
4362 |
|
4363 return retval; |
|
4364 } |
|
4365 |
|
4366 SparseComplexMatrix |
|
4367 SparseMatrix::bsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
4368 int& err, double& rcond, |
|
4369 solve_singularity_handler sing_handler) const |
|
4370 { |
|
4371 SparseComplexMatrix retval; |
|
4372 |
|
4373 int nr = rows (); |
|
4374 int nc = cols (); |
|
4375 err = 0; |
|
4376 |
|
4377 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4378 (*current_liboctave_error_handler) |
|
4379 ("matrix dimension mismatch solution of linear equations"); |
|
4380 else |
|
4381 { |
|
4382 // Print spparms("spumoni") info if requested |
|
4383 volatile int typ = mattype.type (); |
|
4384 mattype.info (); |
|
4385 |
|
4386 if (typ == SparseType::Banded_Hermitian) |
|
4387 { |
|
4388 int n_lower = mattype.nlower (); |
|
4389 int ldm = n_lower + 1; |
|
4390 |
|
4391 Matrix m_band (ldm, nc); |
|
4392 double *tmp_data = m_band.fortran_vec (); |
|
4393 |
|
4394 if (! mattype.is_dense ()) |
|
4395 { |
|
4396 int ii = 0; |
|
4397 |
|
4398 for (int j = 0; j < ldm; j++) |
|
4399 for (int i = 0; i < nc; i++) |
|
4400 tmp_data[ii++] = 0.; |
|
4401 } |
|
4402 |
|
4403 for (int j = 0; j < nc; j++) |
|
4404 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4405 { |
|
4406 int ri = ridx (i); |
|
4407 if (ri >= j) |
|
4408 m_band(ri - j, j) = data(i); |
|
4409 } |
|
4410 |
|
4411 char job = 'L'; |
|
4412 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4413 nr, n_lower, tmp_data, ldm, err |
|
4414 F77_CHAR_ARG_LEN (1))); |
|
4415 |
|
4416 if (f77_exception_encountered) |
|
4417 (*current_liboctave_error_handler) |
|
4418 ("unrecoverable error in dpbtrf"); |
|
4419 else |
|
4420 { |
|
4421 rcond = 0.0; |
|
4422 if (err != 0) |
|
4423 { |
|
4424 // Matrix is not positive definite!! Fall through to |
|
4425 // unsymmetric banded solver. |
|
4426 mattype.mark_as_unsymmetric (); |
|
4427 typ = SparseType::Banded; |
|
4428 |
|
4429 err = 0; |
|
4430 } |
|
4431 else |
|
4432 { |
|
4433 rcond = 1.; |
|
4434 int b_nr = b.rows (); |
|
4435 int b_nc = b.cols (); |
|
4436 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
4437 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
4438 |
|
4439 // Take a first guess that the number of non-zero terms |
|
4440 // will be as many as in b |
|
4441 volatile int x_nz = b.nnz (); |
|
4442 volatile int ii = 0; |
|
4443 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
4444 |
|
4445 retval.xcidx(0) = 0; |
|
4446 for (volatile int j = 0; j < b_nc; j++) |
|
4447 { |
|
4448 |
|
4449 for (int i = 0; i < b_nr; i++) |
|
4450 { |
|
4451 Complex c = b (i,j); |
|
4452 Bx[i] = ::real (c); |
|
4453 Bz[i] = ::imag (c); |
|
4454 } |
|
4455 |
|
4456 F77_XFCN (dpbtrs, DPBTRS, |
|
4457 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4458 nr, n_lower, 1, tmp_data, |
|
4459 ldm, Bx, b_nr, err |
|
4460 F77_CHAR_ARG_LEN (1))); |
|
4461 |
|
4462 if (f77_exception_encountered) |
|
4463 { |
|
4464 (*current_liboctave_error_handler) |
|
4465 ("unrecoverable error in dpbtrs"); |
|
4466 err = -1; |
|
4467 break; |
|
4468 } |
|
4469 |
|
4470 if (err != 0) |
|
4471 { |
|
4472 (*current_liboctave_error_handler) |
|
4473 ("SparseMatrix::solve solve failed"); |
|
4474 err = -1; |
|
4475 break; |
|
4476 } |
|
4477 |
|
4478 F77_XFCN (dpbtrs, DPBTRS, |
|
4479 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4480 nr, n_lower, 1, tmp_data, |
|
4481 ldm, Bz, b_nr, err |
|
4482 F77_CHAR_ARG_LEN (1))); |
|
4483 |
|
4484 if (f77_exception_encountered) |
|
4485 { |
|
4486 (*current_liboctave_error_handler) |
|
4487 ("unrecoverable error in dpbtrs"); |
|
4488 err = -1; |
|
4489 break; |
|
4490 } |
|
4491 |
|
4492 if (err != 0) |
|
4493 { |
|
4494 (*current_liboctave_error_handler) |
|
4495 ("SparseMatrix::solve solve failed"); |
|
4496 |
|
4497 err = -1; |
|
4498 break; |
|
4499 } |
|
4500 |
|
4501 // Count non-zeros in work vector and adjust |
|
4502 // space in retval if needed |
|
4503 int new_nnz = 0; |
|
4504 for (int i = 0; i < nr; i++) |
|
4505 if (Bx[i] != 0. || Bz[i] != 0.) |
|
4506 new_nnz++; |
|
4507 |
|
4508 if (ii + new_nnz > x_nz) |
|
4509 { |
|
4510 // Resize the sparse matrix |
|
4511 int sz = new_nnz * (b_nc - j) + x_nz; |
|
4512 retval.change_capacity (sz); |
|
4513 x_nz = sz; |
|
4514 } |
|
4515 |
|
4516 for (int i = 0; i < nr; i++) |
|
4517 if (Bx[i] != 0. || Bz[i] != 0.) |
|
4518 { |
|
4519 retval.xridx(ii) = i; |
|
4520 retval.xdata(ii++) = |
|
4521 Complex (Bx[i], Bz[i]); |
|
4522 } |
|
4523 |
|
4524 retval.xcidx(j+1) = ii; |
|
4525 } |
|
4526 |
|
4527 retval.maybe_compress (); |
|
4528 } |
|
4529 } |
|
4530 } |
|
4531 |
|
4532 if (typ == SparseType::Banded) |
|
4533 { |
|
4534 // Create the storage for the banded form of the sparse matrix |
|
4535 int n_upper = mattype.nupper (); |
|
4536 int n_lower = mattype.nlower (); |
|
4537 int ldm = n_upper + 2 * n_lower + 1; |
|
4538 |
|
4539 Matrix m_band (ldm, nc); |
|
4540 double *tmp_data = m_band.fortran_vec (); |
|
4541 |
|
4542 if (! mattype.is_dense ()) |
|
4543 { |
|
4544 int ii = 0; |
|
4545 |
|
4546 for (int j = 0; j < ldm; j++) |
|
4547 for (int i = 0; i < nc; i++) |
|
4548 tmp_data[ii++] = 0.; |
|
4549 } |
|
4550 |
|
4551 for (int j = 0; j < nc; j++) |
|
4552 for (int i = cidx(j); i < cidx(j+1); i++) |
|
4553 m_band(ridx(i) - j + n_lower + n_upper, j) = data(i); |
|
4554 |
|
4555 Array<int> ipvt (nr); |
|
4556 int *pipvt = ipvt.fortran_vec (); |
|
4557 |
|
4558 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data, |
|
4559 ldm, pipvt, err)); |
|
4560 |
|
4561 if (f77_exception_encountered) |
|
4562 (*current_liboctave_error_handler) |
|
4563 ("unrecoverable error in dgbtrf"); |
|
4564 else |
|
4565 { |
|
4566 rcond = 0.0; |
|
4567 if (err != 0) |
|
4568 { |
|
4569 err = -2; |
|
4570 |
|
4571 if (sing_handler) |
|
4572 sing_handler (rcond); |
|
4573 else |
|
4574 (*current_liboctave_error_handler) |
|
4575 ("matrix singular to machine precision"); |
|
4576 |
|
4577 } |
|
4578 else |
|
4579 { |
|
4580 char job = 'N'; |
|
4581 volatile int x_nz = b.nnz (); |
|
4582 int b_nc = b.cols (); |
|
4583 retval = SparseComplexMatrix (nr, b_nc, x_nz); |
|
4584 retval.xcidx(0) = 0; |
|
4585 volatile int ii = 0; |
|
4586 |
|
4587 OCTAVE_LOCAL_BUFFER (double, Bx, nr); |
|
4588 OCTAVE_LOCAL_BUFFER (double, Bz, nr); |
|
4589 |
|
4590 for (volatile int j = 0; j < b_nc; j++) |
|
4591 { |
|
4592 for (int i = 0; i < nr; i++) |
|
4593 { |
|
4594 Bx[i] = 0.; |
|
4595 Bz[i] = 0.; |
|
4596 } |
|
4597 for (int i = b.cidx(j); i < b.cidx(j+1); i++) |
|
4598 { |
|
4599 Complex c = b.data(i); |
|
4600 Bx[b.ridx(i)] = ::real (c); |
|
4601 Bz[b.ridx(i)] = ::imag (c); |
|
4602 } |
|
4603 |
|
4604 F77_XFCN (dgbtrs, DGBTRS, |
|
4605 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4606 nr, n_lower, n_upper, 1, tmp_data, |
|
4607 ldm, pipvt, Bx, b.rows (), err |
|
4608 F77_CHAR_ARG_LEN (1))); |
|
4609 |
|
4610 if (f77_exception_encountered) |
|
4611 { |
|
4612 (*current_liboctave_error_handler) |
|
4613 ("unrecoverable error in dgbtrs"); |
|
4614 break; |
|
4615 } |
|
4616 |
|
4617 F77_XFCN (dgbtrs, DGBTRS, |
|
4618 (F77_CONST_CHAR_ARG2 (&job, 1), |
|
4619 nr, n_lower, n_upper, 1, tmp_data, |
|
4620 ldm, pipvt, Bz, b.rows (), err |
|
4621 F77_CHAR_ARG_LEN (1))); |
|
4622 |
|
4623 if (f77_exception_encountered) |
|
4624 { |
|
4625 (*current_liboctave_error_handler) |
|
4626 ("unrecoverable error in dgbtrs"); |
|
4627 break; |
|
4628 } |
|
4629 |
|
4630 // Count non-zeros in work vector and adjust |
|
4631 // space in retval if needed |
|
4632 int new_nnz = 0; |
|
4633 for (int i = 0; i < nr; i++) |
|
4634 if (Bx[i] != 0. || Bz[i] != 0.) |
|
4635 new_nnz++; |
|
4636 |
|
4637 if (ii + new_nnz > x_nz) |
|
4638 { |
|
4639 // Resize the sparse matrix |
|
4640 int sz = new_nnz * (b_nc - j) + x_nz; |
|
4641 retval.change_capacity (sz); |
|
4642 x_nz = sz; |
|
4643 } |
|
4644 |
|
4645 for (int i = 0; i < nr; i++) |
|
4646 if (Bx[i] != 0. || Bz[i] != 0.) |
|
4647 { |
|
4648 retval.xridx(ii) = i; |
|
4649 retval.xdata(ii++) = |
|
4650 Complex (Bx[i], Bz[i]); |
|
4651 } |
|
4652 retval.xcidx(j+1) = ii; |
|
4653 } |
|
4654 |
|
4655 retval.maybe_compress (); |
|
4656 } |
|
4657 } |
|
4658 } |
|
4659 else if (typ != SparseType::Banded_Hermitian) |
|
4660 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4661 } |
|
4662 |
|
4663 return retval; |
|
4664 } |
|
4665 |
|
4666 void * |
|
4667 SparseMatrix::factorize (int& err, double &rcond, Matrix &Control, Matrix &Info, |
|
4668 solve_singularity_handler sing_handler) const |
|
4669 { |
|
4670 // The return values |
|
4671 void *Numeric; |
|
4672 err = 0; |
|
4673 |
|
4674 // Setup the control parameters |
|
4675 Control = Matrix (UMFPACK_CONTROL, 1); |
|
4676 double *control = Control.fortran_vec (); |
|
4677 umfpack_di_defaults (control); |
|
4678 |
|
4679 double tmp = Voctave_sparse_controls.get_key ("spumoni"); |
|
4680 if (!xisnan (tmp)) |
|
4681 Control (UMFPACK_PRL) = tmp; |
|
4682 tmp = Voctave_sparse_controls.get_key ("piv_tol"); |
|
4683 if (!xisnan (tmp)) |
|
4684 { |
|
4685 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; |
|
4686 Control (UMFPACK_PIVOT_TOLERANCE) = tmp; |
|
4687 } |
|
4688 |
|
4689 // Set whether we are allowed to modify Q or not |
|
4690 tmp = Voctave_sparse_controls.get_key ("autoamd"); |
|
4691 if (!xisnan (tmp)) |
|
4692 Control (UMFPACK_FIXQ) = tmp; |
|
4693 |
|
4694 umfpack_di_report_control (control); |
|
4695 |
|
4696 const int *Ap = cidx (); |
|
4697 const int *Ai = ridx (); |
|
4698 const double *Ax = data (); |
|
4699 int nr = rows (); |
|
4700 int nc = cols (); |
|
4701 |
|
4702 umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); |
|
4703 |
|
4704 void *Symbolic; |
|
4705 Info = Matrix (1, UMFPACK_INFO); |
|
4706 double *info = Info.fortran_vec (); |
|
4707 int status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL, |
|
4708 &Symbolic, control, info); |
|
4709 |
|
4710 if (status < 0) |
|
4711 { |
|
4712 (*current_liboctave_error_handler) |
|
4713 ("SparseMatrix::solve symbolic factorization failed"); |
|
4714 err = -1; |
|
4715 |
|
4716 umfpack_di_report_status (control, status); |
|
4717 umfpack_di_report_info (control, info); |
|
4718 |
|
4719 umfpack_di_free_symbolic (&Symbolic) ; |
|
4720 } |
|
4721 else |
|
4722 { |
|
4723 umfpack_di_report_symbolic (Symbolic, control); |
|
4724 |
|
4725 status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, |
|
4726 control, info) ; |
|
4727 umfpack_di_free_symbolic (&Symbolic) ; |
|
4728 |
|
4729 #ifdef HAVE_LSSOLVE |
|
4730 rcond = Info (UMFPACK_RCOND); |
|
4731 volatile double rcond_plus_one = rcond + 1.0; |
|
4732 |
|
4733 if (status == UMFPACK_WARNING_singular_matrix || |
|
4734 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4735 { |
|
4736 umfpack_di_report_numeric (Numeric, control); |
|
4737 |
|
4738 err = -2; |
|
4739 |
|
4740 if (sing_handler) |
|
4741 sing_handler (rcond); |
|
4742 else |
|
4743 (*current_liboctave_error_handler) |
|
4744 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4745 rcond); |
|
4746 |
|
4747 } |
|
4748 else |
|
4749 #endif |
|
4750 if (status < 0) |
|
4751 { |
|
4752 (*current_liboctave_error_handler) |
|
4753 ("SparseMatrix::solve numeric factorization failed"); |
|
4754 |
|
4755 umfpack_di_report_status (control, status); |
|
4756 umfpack_di_report_info (control, info); |
|
4757 |
|
4758 err = -1; |
|
4759 } |
|
4760 else |
|
4761 { |
|
4762 umfpack_di_report_numeric (Numeric, control); |
|
4763 } |
|
4764 } |
|
4765 |
|
4766 if (err != 0) |
|
4767 umfpack_di_free_numeric (&Numeric); |
|
4768 |
|
4769 return Numeric; |
|
4770 } |
|
4771 |
|
4772 Matrix |
|
4773 SparseMatrix::fsolve (SparseType &mattype, const Matrix& b, int& err, |
|
4774 double& rcond, |
|
4775 solve_singularity_handler sing_handler) const |
|
4776 { |
|
4777 Matrix retval; |
|
4778 |
|
4779 int nr = rows (); |
|
4780 int nc = cols (); |
|
4781 err = 0; |
|
4782 |
|
4783 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4784 (*current_liboctave_error_handler) |
|
4785 ("matrix dimension mismatch solution of linear equations"); |
|
4786 else |
|
4787 { |
|
4788 // Print spparms("spumoni") info if requested |
|
4789 int typ = mattype.type (); |
|
4790 mattype.info (); |
|
4791 |
|
4792 if (typ == SparseType::Hermitian) |
|
4793 { |
|
4794 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4795 // through if cholesky factorization fails |
|
4796 |
|
4797 (*current_liboctave_warning_handler) |
|
4798 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4799 |
|
4800 mattype.mark_as_unsymmetric (); |
|
4801 typ = SparseType::Full; |
|
4802 } |
|
4803 |
|
4804 if (typ == SparseType::Full) |
|
4805 { |
|
4806 Matrix Control, Info; |
|
4807 void *Numeric = |
|
4808 factorize (err, rcond, Control, Info, sing_handler); |
|
4809 |
|
4810 if (err == 0) |
|
4811 { |
|
4812 const double *Bx = b.fortran_vec (); |
|
4813 retval.resize (b.rows (), b.cols()); |
|
4814 double *result = retval.fortran_vec (); |
|
4815 int b_nr = b.rows (); |
|
4816 int b_nc = b.cols (); |
|
4817 int status = 0; |
|
4818 double *control = Control.fortran_vec (); |
|
4819 double *info = Info.fortran_vec (); |
|
4820 const int *Ap = cidx (); |
|
4821 const int *Ai = ridx (); |
|
4822 const double *Ax = data (); |
|
4823 |
|
4824 for (int j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) |
|
4825 { |
|
4826 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, |
|
4827 &result[iidx], &Bx[iidx], |
|
4828 Numeric, control, info); |
|
4829 if (status < 0) |
|
4830 { |
|
4831 (*current_liboctave_error_handler) |
|
4832 ("SparseMatrix::solve solve failed"); |
|
4833 |
|
4834 umfpack_di_report_status (control, status); |
|
4835 |
|
4836 err = -1; |
|
4837 |
|
4838 break; |
|
4839 } |
|
4840 } |
|
4841 |
|
4842 #ifndef HAVE_LSSOLVE |
|
4843 rcond = Info (UMFPACK_RCOND); |
|
4844 volatile double rcond_plus_one = rcond + 1.0; |
|
4845 |
|
4846 if (status == UMFPACK_WARNING_singular_matrix || |
|
4847 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4848 { |
|
4849 err = -2; |
|
4850 |
|
4851 if (sing_handler) |
|
4852 sing_handler (rcond); |
|
4853 else |
|
4854 (*current_liboctave_error_handler) |
|
4855 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4856 rcond); |
|
4857 |
|
4858 } |
|
4859 #endif |
|
4860 |
|
4861 umfpack_di_report_info (control, info); |
|
4862 |
|
4863 umfpack_di_free_numeric (&Numeric); |
|
4864 } |
|
4865 } |
|
4866 else if (typ != SparseType::Hermitian) |
|
4867 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
4868 } |
|
4869 |
|
4870 return retval; |
|
4871 } |
|
4872 |
|
4873 SparseMatrix |
|
4874 SparseMatrix::fsolve (SparseType &mattype, const SparseMatrix& b, int& err, double& rcond, |
|
4875 solve_singularity_handler sing_handler) const |
|
4876 { |
|
4877 SparseMatrix retval; |
|
4878 |
|
4879 int nr = rows (); |
|
4880 int nc = cols (); |
|
4881 err = 0; |
|
4882 |
|
4883 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
4884 (*current_liboctave_error_handler) |
|
4885 ("matrix dimension mismatch solution of linear equations"); |
|
4886 else |
|
4887 { |
|
4888 // Print spparms("spumoni") info if requested |
|
4889 int typ = mattype.type (); |
|
4890 mattype.info (); |
|
4891 |
|
4892 if (typ == SparseType::Hermitian) |
|
4893 { |
|
4894 // XXX FIXME XXX Write the cholesky solver and only fall |
|
4895 // through if cholesky factorization fails |
|
4896 |
|
4897 (*current_liboctave_warning_handler) |
|
4898 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
4899 |
|
4900 mattype.mark_as_unsymmetric (); |
|
4901 typ = SparseType::Full; |
|
4902 } |
|
4903 |
|
4904 if (typ == SparseType::Full) |
|
4905 { |
|
4906 Matrix Control, Info; |
|
4907 void *Numeric = factorize (err, rcond, Control, Info, |
|
4908 sing_handler); |
|
4909 |
|
4910 if (err == 0) |
|
4911 { |
|
4912 int b_nr = b.rows (); |
|
4913 int b_nc = b.cols (); |
|
4914 int status = 0; |
|
4915 double *control = Control.fortran_vec (); |
|
4916 double *info = Info.fortran_vec (); |
|
4917 const int *Ap = cidx (); |
|
4918 const int *Ai = ridx (); |
|
4919 const double *Ax = data (); |
|
4920 |
|
4921 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
4922 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr); |
|
4923 |
|
4924 // Take a first guess that the number of non-zero terms |
|
4925 // will be as many as in b |
|
4926 int x_nz = b.nnz (); |
|
4927 int ii = 0; |
|
4928 retval = SparseMatrix (b_nr, b_nc, x_nz); |
|
4929 |
|
4930 retval.xcidx(0) = 0; |
|
4931 for (int j = 0; j < b_nc; j++) |
|
4932 { |
|
4933 |
|
4934 for (int i = 0; i < b_nr; i++) |
|
4935 Bx[i] = b.elem (i, j); |
|
4936 |
|
4937 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, Xx, |
|
4938 Bx, Numeric, control, |
|
4939 info); |
|
4940 if (status < 0) |
|
4941 { |
|
4942 (*current_liboctave_error_handler) |
|
4943 ("SparseMatrix::solve solve failed"); |
|
4944 |
|
4945 umfpack_di_report_status (control, status); |
|
4946 |
|
4947 err = -1; |
|
4948 |
|
4949 break; |
|
4950 } |
|
4951 |
|
4952 for (int i = 0; i < b_nr; i++) |
|
4953 { |
|
4954 double tmp = Xx[i]; |
|
4955 if (tmp != 0.0) |
|
4956 { |
|
4957 if (ii == x_nz) |
|
4958 { |
|
4959 // Resize the sparse matrix |
|
4960 int sz = x_nz * (b_nc - j) / b_nc; |
|
4961 sz = (sz > 10 ? sz : 10) + x_nz; |
|
4962 retval.change_capacity (sz); |
|
4963 x_nz = sz; |
|
4964 } |
|
4965 retval.xdata(ii) = tmp; |
|
4966 retval.xridx(ii++) = i; |
|
4967 } |
|
4968 } |
|
4969 retval.xcidx(j+1) = ii; |
|
4970 } |
|
4971 |
|
4972 retval.maybe_compress (); |
|
4973 |
|
4974 #ifndef HAVE_LSSOLVE |
|
4975 rcond = Info (UMFPACK_RCOND); |
|
4976 volatile double rcond_plus_one = rcond + 1.0; |
|
4977 |
|
4978 if (status == UMFPACK_WARNING_singular_matrix || |
|
4979 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
4980 { |
|
4981 err = -2; |
|
4982 |
|
4983 if (sing_handler) |
|
4984 sing_handler (rcond); |
|
4985 else |
|
4986 (*current_liboctave_error_handler) |
|
4987 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
4988 rcond); |
|
4989 |
|
4990 } |
|
4991 #endif |
|
4992 |
|
4993 umfpack_di_report_info (control, info); |
|
4994 |
|
4995 umfpack_di_free_numeric (&Numeric); |
|
4996 } |
|
4997 } |
|
4998 else if (typ != SparseType::Hermitian) |
|
4999 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
5000 } |
|
5001 |
|
5002 return retval; |
|
5003 } |
|
5004 |
|
5005 ComplexMatrix |
|
5006 SparseMatrix::fsolve (SparseType &mattype, const ComplexMatrix& b, int& err, double& rcond, |
|
5007 solve_singularity_handler sing_handler) const |
|
5008 { |
|
5009 ComplexMatrix retval; |
|
5010 |
|
5011 int nr = rows (); |
|
5012 int nc = cols (); |
|
5013 err = 0; |
|
5014 |
|
5015 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
5016 (*current_liboctave_error_handler) |
|
5017 ("matrix dimension mismatch solution of linear equations"); |
|
5018 else |
|
5019 { |
|
5020 // Print spparms("spumoni") info if requested |
|
5021 int typ = mattype.type (); |
|
5022 mattype.info (); |
|
5023 |
|
5024 if (typ == SparseType::Hermitian) |
|
5025 { |
|
5026 // XXX FIXME XXX Write the cholesky solver and only fall |
|
5027 // through if cholesky factorization fails |
|
5028 |
|
5029 (*current_liboctave_warning_handler) |
|
5030 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
5031 |
|
5032 mattype.mark_as_unsymmetric (); |
|
5033 typ = SparseType::Full; |
|
5034 } |
|
5035 |
|
5036 if (typ == SparseType::Full) |
|
5037 { |
|
5038 Matrix Control, Info; |
|
5039 void *Numeric = factorize (err, rcond, Control, Info, |
|
5040 sing_handler); |
|
5041 |
|
5042 if (err == 0) |
|
5043 { |
|
5044 int b_nr = b.rows (); |
|
5045 int b_nc = b.cols (); |
|
5046 int status = 0; |
|
5047 double *control = Control.fortran_vec (); |
|
5048 double *info = Info.fortran_vec (); |
|
5049 const int *Ap = cidx (); |
|
5050 const int *Ai = ridx (); |
|
5051 const double *Ax = data (); |
|
5052 |
|
5053 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
5054 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
5055 |
|
5056 retval.resize (b_nr, b_nc); |
|
5057 |
|
5058 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr); |
|
5059 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr); |
|
5060 |
|
5061 for (int j = 0; j < b_nc; j++) |
|
5062 { |
|
5063 for (int i = 0; i < b_nr; i++) |
|
5064 { |
|
5065 Complex c = b (i,j); |
|
5066 Bx[i] = ::real (c); |
|
5067 Bz[i] = ::imag (c); |
|
5068 } |
|
5069 |
|
5070 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, |
|
5071 Xx, Bx, Numeric, control, |
|
5072 info); |
|
5073 int status2 = umfpack_di_solve (UMFPACK_A, Ap, Ai, |
|
5074 Ax, Xz, Bz, Numeric, |
|
5075 control, info) ; |
|
5076 |
|
5077 if (status < 0 || status2 < 0) |
|
5078 { |
|
5079 (*current_liboctave_error_handler) |
|
5080 ("SparseMatrix::solve solve failed"); |
|
5081 |
|
5082 umfpack_di_report_status (control, status); |
|
5083 |
|
5084 err = -1; |
|
5085 |
|
5086 break; |
|
5087 } |
|
5088 |
|
5089 for (int i = 0; i < b_nr; i++) |
|
5090 retval (i, j) = Complex (Xx[i], Xz[i]); |
|
5091 } |
|
5092 |
|
5093 #ifndef HAVE_LSSOLVE |
|
5094 rcond = Info (UMFPACK_RCOND); |
|
5095 volatile double rcond_plus_one = rcond + 1.0; |
|
5096 |
|
5097 if (status == UMFPACK_WARNING_singular_matrix || |
|
5098 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
5099 { |
|
5100 err = -2; |
|
5101 |
|
5102 if (sing_handler) |
|
5103 sing_handler (rcond); |
|
5104 else |
|
5105 (*current_liboctave_error_handler) |
|
5106 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
5107 rcond); |
|
5108 |
|
5109 } |
|
5110 #endif |
|
5111 |
|
5112 umfpack_di_report_info (control, info); |
|
5113 |
|
5114 umfpack_di_free_numeric (&Numeric); |
|
5115 } |
|
5116 } |
|
5117 else if (typ != SparseType::Hermitian) |
|
5118 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
5119 } |
|
5120 |
|
5121 return retval; |
|
5122 } |
|
5123 |
|
5124 SparseComplexMatrix |
|
5125 SparseMatrix::fsolve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5126 int& err, double& rcond, |
|
5127 solve_singularity_handler sing_handler) const |
|
5128 { |
|
5129 SparseComplexMatrix retval; |
|
5130 |
|
5131 int nr = rows (); |
|
5132 int nc = cols (); |
|
5133 err = 0; |
|
5134 |
|
5135 if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) |
|
5136 (*current_liboctave_error_handler) |
|
5137 ("matrix dimension mismatch solution of linear equations"); |
|
5138 else |
|
5139 { |
|
5140 // Print spparms("spumoni") info if requested |
|
5141 int typ = mattype.type (); |
|
5142 mattype.info (); |
|
5143 |
|
5144 if (typ == SparseType::Hermitian) |
|
5145 { |
|
5146 // XXX FIXME XXX Write the cholesky solver and only fall |
|
5147 // through if cholesky factorization fails |
|
5148 |
|
5149 (*current_liboctave_warning_handler) |
|
5150 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); |
|
5151 |
|
5152 mattype.mark_as_unsymmetric (); |
|
5153 typ = SparseType::Full; |
|
5154 } |
|
5155 |
|
5156 if (typ == SparseType::Full) |
|
5157 { |
|
5158 Matrix Control, Info; |
|
5159 void *Numeric = factorize (err, rcond, Control, Info, |
|
5160 sing_handler); |
|
5161 |
|
5162 if (err == 0) |
|
5163 { |
|
5164 int b_nr = b.rows (); |
|
5165 int b_nc = b.cols (); |
|
5166 int status = 0; |
|
5167 double *control = Control.fortran_vec (); |
|
5168 double *info = Info.fortran_vec (); |
|
5169 const int *Ap = cidx (); |
|
5170 const int *Ai = ridx (); |
|
5171 const double *Ax = data (); |
|
5172 |
|
5173 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); |
|
5174 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); |
|
5175 |
|
5176 // Take a first guess that the number of non-zero terms |
|
5177 // will be as many as in b |
|
5178 int x_nz = b.nnz (); |
|
5179 int ii = 0; |
|
5180 retval = SparseComplexMatrix (b_nr, b_nc, x_nz); |
|
5181 |
|
5182 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr); |
|
5183 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr); |
|
5184 |
|
5185 retval.xcidx(0) = 0; |
|
5186 for (int j = 0; j < b_nc; j++) |
|
5187 { |
|
5188 for (int i = 0; i < b_nr; i++) |
|
5189 { |
|
5190 Complex c = b (i,j); |
|
5191 Bx[i] = ::real (c); |
|
5192 Bz[i] = ::imag (c); |
|
5193 } |
|
5194 |
|
5195 status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, Xx, |
|
5196 Bx, Numeric, control, |
|
5197 info); |
|
5198 int status2 = umfpack_di_solve (UMFPACK_A, Ap, Ai, |
|
5199 Ax, Xz, Bz, Numeric, |
|
5200 control, info) ; |
|
5201 |
|
5202 if (status < 0 || status2 < 0) |
|
5203 { |
|
5204 (*current_liboctave_error_handler) |
|
5205 ("SparseMatrix::solve solve failed"); |
|
5206 |
|
5207 umfpack_di_report_status (control, status); |
|
5208 |
|
5209 err = -1; |
|
5210 |
|
5211 break; |
|
5212 } |
|
5213 |
|
5214 for (int i = 0; i < b_nr; i++) |
|
5215 { |
|
5216 Complex tmp = Complex (Xx[i], Xz[i]); |
|
5217 if (tmp != 0.0) |
|
5218 { |
|
5219 if (ii == x_nz) |
|
5220 { |
|
5221 // Resize the sparse matrix |
|
5222 int sz = x_nz * (b_nc - j) / b_nc; |
|
5223 sz = (sz > 10 ? sz : 10) + x_nz; |
|
5224 retval.change_capacity (sz); |
|
5225 x_nz = sz; |
|
5226 } |
|
5227 retval.xdata(ii) = tmp; |
|
5228 retval.xridx(ii++) = i; |
|
5229 } |
|
5230 } |
|
5231 retval.xcidx(j+1) = ii; |
|
5232 } |
|
5233 |
|
5234 retval.maybe_compress (); |
|
5235 |
|
5236 #ifndef HAVE_LSSOLVE |
|
5237 rcond = Info (UMFPACK_RCOND); |
|
5238 volatile double rcond_plus_one = rcond + 1.0; |
|
5239 |
|
5240 if (status == UMFPACK_WARNING_singular_matrix || |
|
5241 rcond_plus_one == 1.0 || xisnan (rcond)) |
|
5242 { |
|
5243 err = -2; |
|
5244 |
|
5245 if (sing_handler) |
|
5246 sing_handler (rcond); |
|
5247 else |
|
5248 (*current_liboctave_error_handler) |
|
5249 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", |
|
5250 rcond); |
|
5251 |
|
5252 } |
|
5253 #endif |
|
5254 |
|
5255 umfpack_di_report_info (control, info); |
|
5256 |
|
5257 umfpack_di_free_numeric (&Numeric); |
|
5258 } |
|
5259 } |
|
5260 else if (typ != SparseType::Hermitian) |
|
5261 (*current_liboctave_error_handler) ("incorrect matrix type"); |
|
5262 } |
|
5263 |
|
5264 return retval; |
|
5265 } |
|
5266 |
|
5267 Matrix |
|
5268 SparseMatrix::solve (SparseType &mattype, const Matrix& b) const |
|
5269 { |
|
5270 int info; |
|
5271 double rcond; |
|
5272 return solve (mattype, b, info, rcond, 0); |
|
5273 } |
|
5274 |
|
5275 Matrix |
|
5276 SparseMatrix::solve (SparseType &mattype, const Matrix& b, int& info) const |
|
5277 { |
|
5278 double rcond; |
|
5279 return solve (mattype, b, info, rcond, 0); |
|
5280 } |
|
5281 |
|
5282 Matrix |
|
5283 SparseMatrix::solve (SparseType &mattype, const Matrix& b, int& info, |
|
5284 double& rcond) const |
|
5285 { |
|
5286 return solve (mattype, b, info, rcond, 0); |
|
5287 } |
|
5288 |
|
5289 Matrix |
|
5290 SparseMatrix::solve (SparseType &mattype, const Matrix& b, int& err, |
|
5291 double& rcond, |
|
5292 solve_singularity_handler sing_handler) const |
|
5293 { |
|
5294 int typ = mattype.type (); |
|
5295 |
|
5296 if (typ == SparseType::Unknown) |
|
5297 typ = mattype.type (*this); |
|
5298 |
|
5299 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5300 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5301 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5302 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5303 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5304 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5305 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5306 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5307 else if (typ == SparseType::Tridiagonal || |
|
5308 typ == SparseType::Tridiagonal_Hermitian) |
|
5309 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5310 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5311 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5312 else |
|
5313 { |
|
5314 (*current_liboctave_error_handler) |
|
5315 ("matrix dimension mismatch solution of linear equations"); |
|
5316 return Matrix (); |
|
5317 } |
|
5318 } |
|
5319 |
|
5320 SparseMatrix |
|
5321 SparseMatrix::solve (SparseType &mattype, const SparseMatrix& b) const |
|
5322 { |
|
5323 int info; |
|
5324 double rcond; |
|
5325 return solve (mattype, b, info, rcond, 0); |
|
5326 } |
|
5327 |
|
5328 SparseMatrix |
|
5329 SparseMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5330 int& info) const |
|
5331 { |
|
5332 double rcond; |
|
5333 return solve (mattype, b, info, rcond, 0); |
|
5334 } |
|
5335 |
|
5336 SparseMatrix |
|
5337 SparseMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5338 int& info, double& rcond) const |
|
5339 { |
|
5340 return solve (mattype, b, info, rcond, 0); |
|
5341 } |
|
5342 |
|
5343 SparseMatrix |
|
5344 SparseMatrix::solve (SparseType &mattype, const SparseMatrix& b, |
|
5345 int& err, double& rcond, |
|
5346 solve_singularity_handler sing_handler) const |
|
5347 { |
|
5348 int typ = mattype.type (); |
|
5349 |
|
5350 if (typ == SparseType::Unknown) |
|
5351 typ = mattype.type (*this); |
|
5352 |
|
5353 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5354 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5355 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5356 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5357 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5358 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5359 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5360 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5361 else if (typ == SparseType::Tridiagonal || |
|
5362 typ == SparseType::Tridiagonal_Hermitian) |
|
5363 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5364 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5365 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5366 else |
|
5367 { |
|
5368 (*current_liboctave_error_handler) |
|
5369 ("matrix dimension mismatch solution of linear equations"); |
|
5370 return SparseMatrix (); |
|
5371 } |
|
5372 } |
|
5373 |
|
5374 ComplexMatrix |
|
5375 SparseMatrix::solve (SparseType &mattype, const ComplexMatrix& b) const |
|
5376 { |
|
5377 int info; |
|
5378 double rcond; |
|
5379 return solve (mattype, b, info, rcond, 0); |
|
5380 } |
|
5381 |
|
5382 ComplexMatrix |
|
5383 SparseMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5384 int& info) const |
|
5385 { |
|
5386 double rcond; |
|
5387 return solve (mattype, b, info, rcond, 0); |
|
5388 } |
|
5389 |
|
5390 ComplexMatrix |
|
5391 SparseMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5392 int& info, double& rcond) const |
|
5393 { |
|
5394 return solve (mattype, b, info, rcond, 0); |
|
5395 } |
|
5396 |
|
5397 ComplexMatrix |
|
5398 SparseMatrix::solve (SparseType &mattype, const ComplexMatrix& b, |
|
5399 int& err, double& rcond, |
|
5400 solve_singularity_handler sing_handler) const |
|
5401 { |
|
5402 int typ = mattype.type (); |
|
5403 |
|
5404 if (typ == SparseType::Unknown) |
|
5405 typ = mattype.type (*this); |
|
5406 |
|
5407 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5408 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5409 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5410 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5411 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5412 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5413 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5414 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5415 else if (typ == SparseType::Tridiagonal || |
|
5416 typ == SparseType::Tridiagonal_Hermitian) |
|
5417 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5418 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5419 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5420 else |
|
5421 { |
|
5422 (*current_liboctave_error_handler) |
|
5423 ("matrix dimension mismatch solution of linear equations"); |
|
5424 return ComplexMatrix (); |
|
5425 } |
|
5426 } |
|
5427 |
|
5428 SparseComplexMatrix |
|
5429 SparseMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b) const |
|
5430 { |
|
5431 int info; |
|
5432 double rcond; |
|
5433 return solve (mattype, b, info, rcond, 0); |
|
5434 } |
|
5435 |
|
5436 SparseComplexMatrix |
|
5437 SparseMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5438 int& info) const |
|
5439 { |
|
5440 double rcond; |
|
5441 return solve (mattype, b, info, rcond, 0); |
|
5442 } |
|
5443 |
|
5444 SparseComplexMatrix |
|
5445 SparseMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5446 int& info, double& rcond) const |
|
5447 { |
|
5448 return solve (mattype, b, info, rcond, 0); |
|
5449 } |
|
5450 |
|
5451 SparseComplexMatrix |
|
5452 SparseMatrix::solve (SparseType &mattype, const SparseComplexMatrix& b, |
|
5453 int& err, double& rcond, |
|
5454 solve_singularity_handler sing_handler) const |
|
5455 { |
|
5456 int typ = mattype.type (); |
|
5457 |
|
5458 if (typ == SparseType::Unknown) |
|
5459 typ = mattype.type (*this); |
|
5460 |
|
5461 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) |
|
5462 return dsolve (mattype, b, err, rcond, sing_handler); |
|
5463 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper) |
|
5464 return utsolve (mattype, b, err, rcond, sing_handler); |
|
5465 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower) |
|
5466 return ltsolve (mattype, b, err, rcond, sing_handler); |
|
5467 else if (typ == SparseType::Banded || typ == SparseType::Banded_Hermitian) |
|
5468 return bsolve (mattype, b, err, rcond, sing_handler); |
|
5469 else if (typ == SparseType::Tridiagonal || |
|
5470 typ == SparseType::Tridiagonal_Hermitian) |
|
5471 return trisolve (mattype, b, err, rcond, sing_handler); |
|
5472 else if (typ == SparseType::Full || typ == SparseType::Hermitian) |
|
5473 return fsolve (mattype, b, err, rcond, sing_handler); |
|
5474 else |
|
5475 { |
|
5476 (*current_liboctave_error_handler) |
|
5477 ("matrix dimension mismatch solution of linear equations"); |
|
5478 return SparseComplexMatrix (); |
|
5479 } |
|
5480 } |
|
5481 |
|
5482 ColumnVector |
|
5483 SparseMatrix::solve (SparseType &mattype, const ColumnVector& b) const |
|
5484 { |
|
5485 int info; double rcond; |
|
5486 return solve (mattype, b, info, rcond); |
|
5487 } |
|
5488 |
|
5489 ColumnVector |
|
5490 SparseMatrix::solve (SparseType &mattype, const ColumnVector& b, int& info) const |
|
5491 { |
|
5492 double rcond; |
|
5493 return solve (mattype, b, info, rcond); |
|
5494 } |
|
5495 |
|
5496 ColumnVector |
|
5497 SparseMatrix::solve (SparseType &mattype, const ColumnVector& b, int& info, double& rcond) const |
|
5498 { |
|
5499 return solve (mattype, b, info, rcond, 0); |
|
5500 } |
|
5501 |
|
5502 ColumnVector |
|
5503 SparseMatrix::solve (SparseType &mattype, const ColumnVector& b, int& info, double& rcond, |
|
5504 solve_singularity_handler sing_handler) const |
|
5505 { |
|
5506 Matrix tmp (b); |
|
5507 return solve (mattype, tmp, info, rcond, sing_handler).column (0); |
|
5508 } |
|
5509 |
|
5510 ComplexColumnVector |
|
5511 SparseMatrix::solve (SparseType &mattype, const ComplexColumnVector& b) const |
|
5512 { |
|
5513 int info; |
|
5514 double rcond; |
|
5515 return solve (mattype, b, info, rcond, 0); |
|
5516 } |
|
5517 |
|
5518 ComplexColumnVector |
|
5519 SparseMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, int& info) const |
|
5520 { |
|
5521 double rcond; |
|
5522 return solve (mattype, b, info, rcond, 0); |
|
5523 } |
|
5524 |
|
5525 ComplexColumnVector |
|
5526 SparseMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, int& info, |
|
5527 double& rcond) const |
|
5528 { |
|
5529 return solve (mattype, b, info, rcond, 0); |
|
5530 } |
|
5531 |
|
5532 ComplexColumnVector |
|
5533 SparseMatrix::solve (SparseType &mattype, const ComplexColumnVector& b, int& info, double& rcond, |
|
5534 solve_singularity_handler sing_handler) const |
|
5535 { |
|
5536 ComplexMatrix tmp (b); |
|
5537 return solve (mattype, tmp, info, rcond, sing_handler).column (0); |
|
5538 } |
|
5539 |
|
5540 Matrix |
|
5541 SparseMatrix::solve (const Matrix& b) const |
|
5542 { |
|
5543 int info; |
|
5544 double rcond; |
|
5545 return solve (b, info, rcond, 0); |
|
5546 } |
|
5547 |
|
5548 Matrix |
|
5549 SparseMatrix::solve (const Matrix& b, int& info) const |
|
5550 { |
|
5551 double rcond; |
|
5552 return solve (b, info, rcond, 0); |
|
5553 } |
|
5554 |
|
5555 Matrix |
|
5556 SparseMatrix::solve (const Matrix& b, int& info, |
|
5557 double& rcond) const |
|
5558 { |
|
5559 return solve (b, info, rcond, 0); |
|
5560 } |
|
5561 |
|
5562 Matrix |
|
5563 SparseMatrix::solve (const Matrix& b, int& err, |
|
5564 double& rcond, |
|
5565 solve_singularity_handler sing_handler) const |
|
5566 { |
|
5567 SparseType mattype (*this); |
|
5568 return solve (mattype, b, err, rcond, sing_handler); |
|
5569 } |
|
5570 |
|
5571 SparseMatrix |
|
5572 SparseMatrix::solve (const SparseMatrix& b) const |
|
5573 { |
|
5574 int info; |
|
5575 double rcond; |
|
5576 return solve (b, info, rcond, 0); |
|
5577 } |
|
5578 |
|
5579 SparseMatrix |
|
5580 SparseMatrix::solve (const SparseMatrix& b, |
|
5581 int& info) const |
|
5582 { |
|
5583 double rcond; |
|
5584 return solve (b, info, rcond, 0); |
|
5585 } |
|
5586 |
|
5587 SparseMatrix |
|
5588 SparseMatrix::solve (const SparseMatrix& b, |
|
5589 int& info, double& rcond) const |
|
5590 { |
|
5591 return solve (b, info, rcond, 0); |
|
5592 } |
|
5593 |
|
5594 SparseMatrix |
|
5595 SparseMatrix::solve (const SparseMatrix& b, |
|
5596 int& err, double& rcond, |
|
5597 solve_singularity_handler sing_handler) const |
|
5598 { |
|
5599 SparseType mattype (*this); |
|
5600 return solve (mattype, b, err, rcond, sing_handler); |
|
5601 } |
|
5602 |
|
5603 ComplexMatrix |
|
5604 SparseMatrix::solve (const ComplexMatrix& b, |
|
5605 int& info) const |
|
5606 { |
|
5607 double rcond; |
|
5608 return solve (b, info, rcond, 0); |
|
5609 } |
|
5610 |
|
5611 ComplexMatrix |
|
5612 SparseMatrix::solve (const ComplexMatrix& b, |
|
5613 int& info, double& rcond) const |
|
5614 { |
|
5615 return solve (b, info, rcond, 0); |
|
5616 } |
|
5617 |
|
5618 ComplexMatrix |
|
5619 SparseMatrix::solve (const ComplexMatrix& b, |
|
5620 int& err, double& rcond, |
|
5621 solve_singularity_handler sing_handler) const |
|
5622 { |
|
5623 SparseType mattype (*this); |
|
5624 return solve (mattype, b, err, rcond, sing_handler); |
|
5625 } |
|
5626 |
|
5627 SparseComplexMatrix |
|
5628 SparseMatrix::solve (const SparseComplexMatrix& b) const |
|
5629 { |
|
5630 int info; |
|
5631 double rcond; |
|
5632 return solve (b, info, rcond, 0); |
|
5633 } |
|
5634 |
|
5635 SparseComplexMatrix |
|
5636 SparseMatrix::solve (const SparseComplexMatrix& b, |
|
5637 int& info) const |
|
5638 { |
|
5639 double rcond; |
|
5640 return solve (b, info, rcond, 0); |
|
5641 } |
|
5642 |
|
5643 SparseComplexMatrix |
|
5644 SparseMatrix::solve (const SparseComplexMatrix& b, |
|
5645 int& info, double& rcond) const |
|
5646 { |
|
5647 return solve (b, info, rcond, 0); |
|
5648 } |
|
5649 |
|
5650 SparseComplexMatrix |
|
5651 SparseMatrix::solve (const SparseComplexMatrix& b, |
|
5652 int& err, double& rcond, |
|
5653 solve_singularity_handler sing_handler) const |
|
5654 { |
|
5655 SparseType mattype (*this); |
|
5656 return solve (mattype, b, err, rcond, sing_handler); |
|
5657 } |
|
5658 |
|
5659 ColumnVector |
|
5660 SparseMatrix::solve (const ColumnVector& b) const |
|
5661 { |
|
5662 int info; double rcond; |
|
5663 return solve (b, info, rcond); |
|
5664 } |
|
5665 |
|
5666 ColumnVector |
|
5667 SparseMatrix::solve (const ColumnVector& b, int& info) const |
|
5668 { |
|
5669 double rcond; |
|
5670 return solve (b, info, rcond); |
|
5671 } |
|
5672 |
|
5673 ColumnVector |
|
5674 SparseMatrix::solve (const ColumnVector& b, int& info, double& rcond) const |
|
5675 { |
|
5676 return solve (b, info, rcond, 0); |
|
5677 } |
|
5678 |
|
5679 ColumnVector |
|
5680 SparseMatrix::solve (const ColumnVector& b, int& info, double& rcond, |
|
5681 solve_singularity_handler sing_handler) const |
|
5682 { |
|
5683 Matrix tmp (b); |
|
5684 return solve (tmp, info, rcond, sing_handler).column (0); |
|
5685 } |
|
5686 |
|
5687 ComplexColumnVector |
|
5688 SparseMatrix::solve (const ComplexColumnVector& b) const |
|
5689 { |
|
5690 int info; |
|
5691 double rcond; |
|
5692 return solve (b, info, rcond, 0); |
|
5693 } |
|
5694 |
|
5695 ComplexColumnVector |
|
5696 SparseMatrix::solve (const ComplexColumnVector& b, int& info) const |
|
5697 { |
|
5698 double rcond; |
|
5699 return solve (b, info, rcond, 0); |
|
5700 } |
|
5701 |
|
5702 ComplexColumnVector |
|
5703 SparseMatrix::solve (const ComplexColumnVector& b, int& info, |
|
5704 double& rcond) const |
|
5705 { |
|
5706 return solve (b, info, rcond, 0); |
|
5707 } |
|
5708 |
|
5709 ComplexColumnVector |
|
5710 SparseMatrix::solve (const ComplexColumnVector& b, int& info, double& rcond, |
|
5711 solve_singularity_handler sing_handler) const |
|
5712 { |
|
5713 ComplexMatrix tmp (b); |
|
5714 return solve (tmp, info, rcond, sing_handler).column (0); |
|
5715 } |
|
5716 |
|
5717 Matrix |
|
5718 SparseMatrix::lssolve (const Matrix& b) const |
|
5719 { |
|
5720 int info; |
|
5721 int rank; |
|
5722 return lssolve (b, info, rank); |
|
5723 } |
|
5724 |
|
5725 Matrix |
|
5726 SparseMatrix::lssolve (const Matrix& b, int& info) const |
|
5727 { |
|
5728 int rank; |
|
5729 return lssolve (b, info, rank); |
|
5730 } |
|
5731 |
|
5732 Matrix |
|
5733 SparseMatrix::lssolve (const Matrix& b, int& info, int& rank) const |
|
5734 { |
|
5735 info = -1; |
|
5736 (*current_liboctave_error_handler) |
|
5737 ("SparseMatrix::lssolve not implemented yet"); |
|
5738 return Matrix (); |
|
5739 } |
|
5740 |
|
5741 SparseMatrix |
|
5742 SparseMatrix::lssolve (const SparseMatrix& b) const |
|
5743 { |
|
5744 int info; |
|
5745 int rank; |
|
5746 return lssolve (b, info, rank); |
|
5747 } |
|
5748 |
|
5749 SparseMatrix |
|
5750 SparseMatrix::lssolve (const SparseMatrix& b, int& info) const |
|
5751 { |
|
5752 int rank; |
|
5753 return lssolve (b, info, rank); |
|
5754 } |
|
5755 |
|
5756 SparseMatrix |
|
5757 SparseMatrix::lssolve (const SparseMatrix& b, int& info, int& rank) const |
|
5758 { |
|
5759 info = -1; |
|
5760 (*current_liboctave_error_handler) |
|
5761 ("SparseMatrix::lssolve not implemented yet"); |
|
5762 return SparseMatrix (); |
|
5763 } |
|
5764 |
|
5765 ComplexMatrix |
|
5766 SparseMatrix::lssolve (const ComplexMatrix& b) const |
|
5767 { |
|
5768 int info; |
|
5769 int rank; |
|
5770 return lssolve (b, info, rank); |
|
5771 } |
|
5772 |
|
5773 ComplexMatrix |
|
5774 SparseMatrix::lssolve (const ComplexMatrix& b, int& info) const |
|
5775 { |
|
5776 int rank; |
|
5777 return lssolve (b, info, rank); |
|
5778 } |
|
5779 |
|
5780 ComplexMatrix |
|
5781 SparseMatrix::lssolve (const ComplexMatrix& b, int& info, int& rank) const |
|
5782 { |
|
5783 info = -1; |
|
5784 (*current_liboctave_error_handler) |
|
5785 ("SparseMatrix::lssolve not implemented yet"); |
|
5786 return ComplexMatrix (); |
|
5787 } |
|
5788 |
|
5789 SparseComplexMatrix |
|
5790 SparseMatrix::lssolve (const SparseComplexMatrix& b) const |
|
5791 { |
|
5792 int info; |
|
5793 int rank; |
|
5794 return lssolve (b, info, rank); |
|
5795 } |
|
5796 |
|
5797 SparseComplexMatrix |
|
5798 SparseMatrix::lssolve (const SparseComplexMatrix& b, int& info) const |
|
5799 { |
|
5800 int rank; |
|
5801 return lssolve (b, info, rank); |
|
5802 } |
|
5803 |
|
5804 SparseComplexMatrix |
|
5805 SparseMatrix::lssolve (const SparseComplexMatrix& b, int& info, |
|
5806 int& rank) const |
|
5807 { |
|
5808 info = -1; |
|
5809 (*current_liboctave_error_handler) |
|
5810 ("SparseMatrix::lssolve not implemented yet"); |
|
5811 return SparseComplexMatrix (); |
|
5812 } |
|
5813 |
|
5814 ColumnVector |
|
5815 SparseMatrix::lssolve (const ColumnVector& b) const |
|
5816 { |
|
5817 int info; |
|
5818 int rank; |
|
5819 return lssolve (b, info, rank); |
|
5820 } |
|
5821 |
|
5822 ColumnVector |
|
5823 SparseMatrix::lssolve (const ColumnVector& b, int& info) const |
|
5824 { |
|
5825 int rank; |
|
5826 return lssolve (b, info, rank); |
|
5827 } |
|
5828 |
|
5829 ColumnVector |
|
5830 SparseMatrix::lssolve (const ColumnVector& b, int& info, int& rank) const |
|
5831 { |
|
5832 Matrix tmp (b); |
|
5833 return lssolve (tmp, info, rank).column (0); |
|
5834 } |
|
5835 |
|
5836 ComplexColumnVector |
|
5837 SparseMatrix::lssolve (const ComplexColumnVector& b) const |
|
5838 { |
|
5839 int info; |
|
5840 int rank; |
|
5841 return lssolve (b, info, rank); |
|
5842 } |
|
5843 |
|
5844 ComplexColumnVector |
|
5845 SparseMatrix::lssolve (const ComplexColumnVector& b, int& info) const |
|
5846 { |
|
5847 int rank; |
|
5848 return lssolve (b, info, rank); |
|
5849 } |
|
5850 |
|
5851 ComplexColumnVector |
|
5852 SparseMatrix::lssolve (const ComplexColumnVector& b, int& info, |
|
5853 int& rank) const |
|
5854 { |
|
5855 ComplexMatrix tmp (b); |
|
5856 return lssolve (tmp, info, rank).column (0); |
|
5857 } |
|
5858 |
|
5859 // other operations. |
|
5860 |
|
5861 SparseMatrix |
|
5862 SparseMatrix::map (d_d_Mapper f) const |
|
5863 { |
|
5864 int nr = rows (); |
|
5865 int nc = cols (); |
|
5866 int nz = nnz (); |
|
5867 bool f_zero = (f(0.0) == 0.0); |
|
5868 |
|
5869 // Count number of non-zero elements |
|
5870 int nel = (f_zero ? 0 : nr*nc - nz); |
|
5871 for (int i = 0; i < nz; i++) |
|
5872 if (f (data(i)) != 0.0) |
|
5873 nel++; |
|
5874 |
|
5875 SparseMatrix retval (nr, nc, nel); |
|
5876 |
|
5877 if (f_zero) |
|
5878 { |
|
5879 int ii = 0; |
|
5880 for (int j = 0; j < nc; j++) |
|
5881 { |
|
5882 for (int i = 0; i < nr; i++) |
|
5883 { |
|
5884 double tmp = f (elem (i, j)); |
|
5885 if (tmp != 0.0) |
|
5886 { |
|
5887 retval.data(ii) = tmp; |
|
5888 retval.ridx(ii++) = i; |
|
5889 } |
|
5890 } |
|
5891 retval.cidx(j+1) = ii; |
|
5892 } |
|
5893 } |
|
5894 else |
|
5895 { |
|
5896 int ii = 0; |
|
5897 for (int j = 0; j < nc; j++) |
|
5898 { |
|
5899 for (int i = cidx(j); i < cidx(j+1); i++) |
|
5900 { |
|
5901 retval.data(ii) = f (elem(i)); |
|
5902 retval.ridx(ii++) = ridx(i); |
|
5903 } |
|
5904 retval.cidx(j+1) = ii; |
|
5905 } |
|
5906 } |
|
5907 |
|
5908 return retval; |
|
5909 } |
|
5910 |
|
5911 SparseBoolMatrix |
|
5912 SparseMatrix::map (b_d_Mapper f) const |
|
5913 { |
|
5914 int nr = rows (); |
|
5915 int nc = cols (); |
|
5916 int nz = nnz (); |
|
5917 bool f_zero = f(0.0); |
|
5918 |
|
5919 // Count number of non-zero elements |
|
5920 int nel = (f_zero ? 0 : nr*nc - nz); |
|
5921 for (int i = 0; i < nz; i++) |
|
5922 if (f (data(i)) != 0.0) |
|
5923 nel++; |
|
5924 |
|
5925 SparseBoolMatrix retval (nr, nc, nel); |
|
5926 |
|
5927 if (f_zero) |
|
5928 { |
|
5929 int ii = 0; |
|
5930 for (int j = 0; j < nc; j++) |
|
5931 { |
|
5932 for (int i = 0; i < nr; i++) |
|
5933 { |
|
5934 bool tmp = f (elem (i, j)); |
|
5935 if (tmp) |
|
5936 { |
|
5937 retval.data(ii) = tmp; |
|
5938 retval.ridx(ii++) = i; |
|
5939 } |
|
5940 } |
|
5941 retval.cidx(j+1) = ii; |
|
5942 } |
|
5943 } |
|
5944 else |
|
5945 { |
|
5946 int ii = 0; |
|
5947 for (int j = 0; j < nc; j++) |
|
5948 { |
|
5949 for (int i = cidx(j); i < cidx(j+1); i++) |
|
5950 { |
|
5951 retval.data(ii) = f (elem(i)); |
|
5952 retval.ridx(ii++) = ridx(i); |
|
5953 } |
|
5954 retval.cidx(j+1) = ii; |
|
5955 } |
|
5956 } |
|
5957 |
|
5958 return retval; |
|
5959 } |
|
5960 |
|
5961 SparseMatrix& |
|
5962 SparseMatrix::apply (d_d_Mapper f) |
|
5963 { |
|
5964 *this = map (f); |
|
5965 return *this; |
|
5966 } |
|
5967 |
|
5968 bool |
|
5969 SparseMatrix::any_element_is_negative (bool neg_zero) const |
|
5970 { |
|
5971 int nel = nnz (); |
|
5972 |
|
5973 if (neg_zero) |
|
5974 { |
|
5975 for (int i = 0; i < nel; i++) |
|
5976 if (lo_ieee_signbit (data (i))) |
|
5977 return true; |
|
5978 } |
|
5979 else |
|
5980 { |
|
5981 for (int i = 0; i < nel; i++) |
|
5982 if (data (i) < 0) |
|
5983 return true; |
|
5984 } |
|
5985 |
|
5986 return false; |
|
5987 } |
|
5988 |
|
5989 bool |
|
5990 SparseMatrix::any_element_is_inf_or_nan (void) const |
|
5991 { |
|
5992 int nel = nnz (); |
|
5993 |
|
5994 for (int i = 0; i < nel; i++) |
|
5995 { |
|
5996 double val = data (i); |
|
5997 if (xisinf (val) || xisnan (val)) |
|
5998 return true; |
|
5999 } |
|
6000 |
|
6001 return false; |
|
6002 } |
|
6003 |
|
6004 bool |
|
6005 SparseMatrix::all_elements_are_int_or_inf_or_nan (void) const |
|
6006 { |
|
6007 int nel = nnz (); |
|
6008 |
|
6009 for (int i = 0; i < nel; i++) |
|
6010 { |
|
6011 double val = data (i); |
|
6012 if (xisnan (val) || D_NINT (val) == val) |
|
6013 continue; |
|
6014 else |
|
6015 return false; |
|
6016 } |
|
6017 |
|
6018 return true; |
|
6019 } |
|
6020 |
|
6021 // Return nonzero if any element of M is not an integer. Also extract |
|
6022 // the largest and smallest values and return them in MAX_VAL and MIN_VAL. |
|
6023 |
|
6024 bool |
|
6025 SparseMatrix::all_integers (double& max_val, double& min_val) const |
|
6026 { |
|
6027 int nel = nnz (); |
|
6028 |
|
6029 if (nel == 0) |
|
6030 return false; |
|
6031 |
|
6032 max_val = data (0); |
|
6033 min_val = data (0); |
|
6034 |
|
6035 for (int i = 0; i < nel; i++) |
|
6036 { |
|
6037 double val = data (i); |
|
6038 |
|
6039 if (val > max_val) |
|
6040 max_val = val; |
|
6041 |
|
6042 if (val < min_val) |
|
6043 min_val = val; |
|
6044 |
|
6045 if (D_NINT (val) != val) |
|
6046 return false; |
|
6047 } |
|
6048 |
|
6049 return true; |
|
6050 } |
|
6051 |
|
6052 bool |
|
6053 SparseMatrix::too_large_for_float (void) const |
|
6054 { |
|
6055 int nel = nnz (); |
|
6056 |
|
6057 for (int i = 0; i < nel; i++) |
|
6058 { |
|
6059 double val = data (i); |
|
6060 |
|
6061 if (val > FLT_MAX || val < FLT_MIN) |
|
6062 return true; |
|
6063 } |
|
6064 |
|
6065 return false; |
|
6066 } |
|
6067 |
|
6068 SparseBoolMatrix |
|
6069 SparseMatrix::operator ! (void) const |
|
6070 { |
|
6071 int nr = rows (); |
|
6072 int nc = cols (); |
|
6073 int nz1 = nnz (); |
|
6074 int nz2 = nr*nc - nz1; |
|
6075 |
|
6076 SparseBoolMatrix r (nr, nc, nz2); |
|
6077 |
|
6078 int ii = 0; |
|
6079 int jj = 0; |
|
6080 r.cidx (0) = 0; |
|
6081 for (int i = 0; i < nc; i++) |
|
6082 { |
|
6083 for (int j = 0; j < nr; j++) |
|
6084 { |
|
6085 if (jj < cidx(i+1) && ridx(jj) == j) |
|
6086 jj++; |
|
6087 else |
|
6088 { |
|
6089 r.data(ii) = true; |
|
6090 r.ridx(ii++) = j; |
|
6091 } |
|
6092 } |
|
6093 r.cidx (i+1) = ii; |
|
6094 } |
|
6095 |
|
6096 return r; |
|
6097 } |
|
6098 |
|
6099 // XXX FIXME XXX Do these really belong here? Maybe they should be |
|
6100 // in a base class? |
|
6101 |
|
6102 SparseBoolMatrix |
|
6103 SparseMatrix::all (int dim) const |
|
6104 { |
|
6105 SPARSE_ALL_OP (dim); |
|
6106 } |
|
6107 |
|
6108 SparseBoolMatrix |
|
6109 SparseMatrix::any (int dim) const |
|
6110 { |
|
6111 SPARSE_ANY_OP (dim); |
|
6112 } |
|
6113 |
|
6114 SparseMatrix |
|
6115 SparseMatrix::cumprod (int dim) const |
|
6116 { |
|
6117 SPARSE_CUMPROD (SparseMatrix, double, cumprod); |
|
6118 } |
|
6119 |
|
6120 SparseMatrix |
|
6121 SparseMatrix::cumsum (int dim) const |
|
6122 { |
|
6123 SPARSE_CUMSUM (SparseMatrix, double, cumsum); |
|
6124 } |
|
6125 |
|
6126 SparseMatrix |
|
6127 SparseMatrix::prod (int dim) const |
|
6128 { |
|
6129 SPARSE_REDUCTION_OP (SparseMatrix, double, *=, 1.0, 1.0); |
|
6130 } |
|
6131 |
|
6132 SparseMatrix |
|
6133 SparseMatrix::sum (int dim) const |
|
6134 { |
|
6135 SPARSE_REDUCTION_OP (SparseMatrix, double, +=, 0.0, 0.0); |
|
6136 } |
|
6137 |
|
6138 SparseMatrix |
|
6139 SparseMatrix::sumsq (int dim) const |
|
6140 { |
|
6141 #define ROW_EXPR \ |
|
6142 double d = elem (i, j); \ |
|
6143 tmp[i] += d * d |
|
6144 |
|
6145 #define COL_EXPR \ |
|
6146 double d = elem (i, j); \ |
|
6147 tmp[j] += d * d |
|
6148 |
|
6149 SPARSE_BASE_REDUCTION_OP (SparseMatrix, double, ROW_EXPR, COL_EXPR, |
|
6150 0.0, 0.0); |
|
6151 |
|
6152 #undef ROW_EXPR |
|
6153 #undef COL_EXPR |
|
6154 } |
|
6155 |
|
6156 SparseMatrix |
|
6157 SparseMatrix::abs (void) const |
|
6158 { |
|
6159 int nz = nnz (); |
|
6160 |
|
6161 SparseMatrix retval (*this); |
|
6162 |
|
6163 for (int i = 0; i < nz; i++) |
|
6164 retval.data(i) = fabs(retval.data(i)); |
|
6165 |
|
6166 return retval; |
|
6167 } |
|
6168 |
|
6169 SparseMatrix |
|
6170 SparseMatrix::diag (int k) const |
|
6171 { |
|
6172 int nnr = rows (); |
|
6173 int nnc = cols (); |
|
6174 |
|
6175 if (k > 0) |
|
6176 nnc -= k; |
|
6177 else if (k < 0) |
|
6178 nnr += k; |
|
6179 |
|
6180 SparseMatrix d; |
|
6181 |
|
6182 if (nnr > 0 && nnc > 0) |
|
6183 { |
|
6184 int ndiag = (nnr < nnc) ? nnr : nnc; |
|
6185 |
|
6186 // Count the number of non-zero elements |
|
6187 int nel = 0; |
|
6188 if (k > 0) |
|
6189 { |
|
6190 for (int i = 0; i < ndiag; i++) |
|
6191 if (elem (i, i+k) != 0.) |
|
6192 nel++; |
|
6193 } |
|
6194 else if ( k < 0) |
|
6195 { |
|
6196 for (int i = 0; i < ndiag; i++) |
|
6197 if (elem (i-k, i) != 0.) |
|
6198 nel++; |
|
6199 } |
|
6200 else |
|
6201 { |
|
6202 for (int i = 0; i < ndiag; i++) |
|
6203 if (elem (i, i) != 0.) |
|
6204 nel++; |
|
6205 } |
|
6206 |
|
6207 d = SparseMatrix (ndiag, 1, nel); |
|
6208 d.xcidx (0) = 0; |
|
6209 d.xcidx (1) = nel; |
|
6210 |
|
6211 int ii = 0; |
|
6212 if (k > 0) |
|
6213 { |
|
6214 for (int i = 0; i < ndiag; i++) |
|
6215 { |
|
6216 double tmp = elem (i, i+k); |
|
6217 if (tmp != 0.) |
|
6218 { |
|
6219 d.xdata (ii) = tmp; |
|
6220 d.xridx (ii++) = i; |
|
6221 } |
|
6222 } |
|
6223 } |
|
6224 else if ( k < 0) |
|
6225 { |
|
6226 for (int i = 0; i < ndiag; i++) |
|
6227 { |
|
6228 double tmp = elem (i-k, i); |
|
6229 if (tmp != 0.) |
|
6230 { |
|
6231 d.xdata (ii) = tmp; |
|
6232 d.xridx (ii++) = i; |
|
6233 } |
|
6234 } |
|
6235 } |
|
6236 else |
|
6237 { |
|
6238 for (int i = 0; i < ndiag; i++) |
|
6239 { |
|
6240 double tmp = elem (i, i); |
|
6241 if (tmp != 0.) |
|
6242 { |
|
6243 d.xdata (ii) = tmp; |
|
6244 d.xridx (ii++) = i; |
|
6245 } |
|
6246 } |
|
6247 } |
|
6248 } |
|
6249 else |
|
6250 (*current_liboctave_error_handler) |
|
6251 ("diag: requested diagonal out of range"); |
|
6252 |
|
6253 return d; |
|
6254 } |
|
6255 |
|
6256 Matrix |
|
6257 SparseMatrix::matrix_value (void) const |
|
6258 { |
|
6259 int nr = rows (); |
|
6260 int nc = cols (); |
|
6261 |
|
6262 Matrix retval (nr, nc, 0.0); |
|
6263 for (int j = 0; j < nc; j++) |
|
6264 for (int i = cidx(j); i < cidx(j+1); i++) |
|
6265 retval.elem (ridx(i), j) = data (i); |
|
6266 |
|
6267 return retval; |
|
6268 } |
|
6269 |
|
6270 std::ostream& |
|
6271 operator << (std::ostream& os, const SparseMatrix& a) |
|
6272 { |
|
6273 int nc = a.cols (); |
|
6274 |
|
6275 // add one to the printed indices to go from |
|
6276 // zero-based to one-based arrays |
|
6277 for (int j = 0; j < nc; j++) { |
|
6278 OCTAVE_QUIT; |
|
6279 for (int i = a.cidx(j); i < a.cidx(j+1); i++) { |
|
6280 os << a.ridx(i) + 1 << " " << j + 1 << " "; |
|
6281 octave_write_double (os, a.data(i)); |
|
6282 os << "\n"; |
|
6283 } |
|
6284 } |
|
6285 |
|
6286 return os; |
|
6287 } |
|
6288 |
|
6289 std::istream& |
|
6290 operator >> (std::istream& is, SparseMatrix& a) |
|
6291 { |
|
6292 int nr = a.rows (); |
|
6293 int nc = a.cols (); |
|
6294 int nz = a.nnz (); |
|
6295 |
|
6296 if (nr < 1 || nc < 1) |
|
6297 is.clear (std::ios::badbit); |
|
6298 else |
|
6299 { |
|
6300 int itmp, jtmp, jold = 0; |
|
6301 double tmp; |
|
6302 int ii = 0; |
|
6303 |
|
6304 a.cidx (0) = 0; |
|
6305 for (int i = 0; i < nz; i++) |
|
6306 { |
|
6307 is >> itmp; |
|
6308 itmp--; |
|
6309 is >> jtmp; |
|
6310 jtmp--; |
|
6311 tmp = octave_read_double (is); |
|
6312 |
|
6313 if (is) |
|
6314 { |
|
6315 if (jold != jtmp) |
|
6316 { |
|
6317 for (int j = jold; j < jtmp; j++) |
|
6318 a.cidx(j+1) = ii; |
|
6319 |
|
6320 jold = jtmp; |
|
6321 } |
|
6322 a.data (ii) = tmp; |
|
6323 a.ridx (ii++) = itmp; |
|
6324 } |
|
6325 else |
|
6326 goto done; |
|
6327 } |
|
6328 |
|
6329 for (int j = jold; j < nc; j++) |
|
6330 a.cidx(j+1) = ii; |
|
6331 } |
|
6332 |
|
6333 done: |
|
6334 |
|
6335 return is; |
|
6336 } |
|
6337 |
|
6338 SparseMatrix |
|
6339 SparseMatrix::squeeze (void) const |
|
6340 { |
|
6341 return MSparse<double>::squeeze (); |
|
6342 } |
|
6343 |
|
6344 SparseMatrix |
|
6345 SparseMatrix::index (idx_vector& i, int resize_ok) const |
|
6346 { |
|
6347 return MSparse<double>::index (i, resize_ok); |
|
6348 } |
|
6349 |
|
6350 SparseMatrix |
|
6351 SparseMatrix::index (idx_vector& i, idx_vector& j, int resize_ok) const |
|
6352 { |
|
6353 return MSparse<double>::index (i, j, resize_ok); |
|
6354 } |
|
6355 |
|
6356 SparseMatrix |
|
6357 SparseMatrix::index (Array<idx_vector>& ra_idx, int resize_ok) const |
|
6358 { |
|
6359 return MSparse<double>::index (ra_idx, resize_ok); |
|
6360 } |
|
6361 |
|
6362 SparseMatrix |
|
6363 SparseMatrix::reshape (const dim_vector& new_dims) const |
|
6364 { |
|
6365 return MSparse<double>::reshape (new_dims); |
|
6366 } |
|
6367 |
|
6368 SparseMatrix |
|
6369 SparseMatrix::permute (const Array<int>& vec, bool inv) const |
|
6370 { |
|
6371 return MSparse<double>::permute (vec, inv); |
|
6372 } |
|
6373 |
|
6374 SparseMatrix |
|
6375 SparseMatrix::ipermute (const Array<int>& vec) const |
|
6376 { |
|
6377 return MSparse<double>::ipermute (vec); |
|
6378 } |
|
6379 |
|
6380 // matrix by matrix -> matrix operations |
|
6381 |
|
6382 SparseMatrix |
|
6383 operator * (const SparseMatrix& m, const SparseMatrix& a) |
|
6384 { |
|
6385 #ifdef HAVE_SPARSE_BLAS |
|
6386 // XXX FIXME XXX Isn't there a sparse BLAS ?? |
|
6387 #else |
|
6388 // Use Andy's sparse matrix multiply function |
|
6389 SPARSE_SPARSE_MUL (SparseMatrix, double); |
|
6390 #endif |
|
6391 } |
|
6392 |
|
6393 // XXX FIXME XXX -- it would be nice to share code among the min/max |
|
6394 // functions below. |
|
6395 |
|
6396 #define EMPTY_RETURN_CHECK(T) \ |
|
6397 if (nr == 0 || nc == 0) \ |
|
6398 return T (nr, nc); |
|
6399 |
|
6400 SparseMatrix |
|
6401 min (double d, const SparseMatrix& m) |
|
6402 { |
|
6403 SparseMatrix result; |
|
6404 |
|
6405 int nr = m.rows (); |
|
6406 int nc = m.columns (); |
|
6407 |
|
6408 EMPTY_RETURN_CHECK (SparseMatrix); |
|
6409 |
|
6410 // Count the number of non-zero elements |
|
6411 if (d < 0.) |
|
6412 { |
|
6413 result = SparseMatrix (nr, nc, d); |
|
6414 for (int j = 0; j < nc; j++) |
|
6415 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6416 { |
|
6417 double tmp = xmin (d, m.data (i)); |
|
6418 if (tmp != 0.) |
|
6419 { |
|
6420 int idx = m.ridx(i) + j * nr; |
|
6421 result.xdata(idx) = tmp; |
|
6422 result.xridx(idx) = m.ridx(i); |
|
6423 } |
|
6424 } |
|
6425 } |
|
6426 else |
|
6427 { |
|
6428 int nel = 0; |
|
6429 for (int j = 0; j < nc; j++) |
|
6430 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6431 if (xmin (d, m.data (i)) != 0.) |
|
6432 nel++; |
|
6433 |
|
6434 result = SparseMatrix (nr, nc, nel); |
|
6435 |
|
6436 int ii = 0; |
|
6437 result.xcidx(0) = 0; |
|
6438 for (int j = 0; j < nc; j++) |
|
6439 { |
|
6440 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6441 { |
|
6442 double tmp = xmin (d, m.data (i)); |
|
6443 |
|
6444 if (tmp != 0.) |
|
6445 { |
|
6446 result.xdata(ii) = tmp; |
|
6447 result.xridx(ii++) = m.ridx(i); |
|
6448 } |
|
6449 } |
|
6450 result.xcidx(j+1) = ii; |
|
6451 } |
|
6452 } |
|
6453 |
|
6454 return result; |
|
6455 } |
|
6456 |
|
6457 SparseMatrix |
|
6458 min (const SparseMatrix& m, double d) |
|
6459 { |
|
6460 return min (d, m); |
|
6461 } |
|
6462 |
|
6463 SparseMatrix |
|
6464 min (const SparseMatrix& a, const SparseMatrix& b) |
|
6465 { |
|
6466 SparseMatrix r; |
|
6467 |
|
6468 if ((a.rows() == b.rows()) && (a.cols() == b.cols())) |
|
6469 { |
|
6470 int a_nr = a.rows (); |
|
6471 int a_nc = a.cols (); |
|
6472 |
|
6473 int b_nr = b.rows (); |
|
6474 int b_nc = b.cols (); |
|
6475 |
|
6476 if (a_nr != b_nr || a_nc != b_nc) |
|
6477 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc); |
|
6478 else |
|
6479 { |
|
6480 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
6481 |
|
6482 int jx = 0; |
|
6483 r.cidx (0) = 0; |
|
6484 for (int i = 0 ; i < a_nc ; i++) |
|
6485 { |
|
6486 int ja = a.cidx(i); |
|
6487 int ja_max = a.cidx(i+1); |
|
6488 bool ja_lt_max= ja < ja_max; |
|
6489 |
|
6490 int jb = b.cidx(i); |
|
6491 int jb_max = b.cidx(i+1); |
|
6492 bool jb_lt_max = jb < jb_max; |
|
6493 |
|
6494 while (ja_lt_max || jb_lt_max ) |
|
6495 { |
|
6496 OCTAVE_QUIT; |
|
6497 if ((! jb_lt_max) || |
|
6498 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
6499 { |
|
6500 double tmp = xmin (a.data(ja), 0.); |
|
6501 if (tmp != 0.) |
|
6502 { |
|
6503 r.ridx(jx) = a.ridx(ja); |
|
6504 r.data(jx) = tmp; |
|
6505 jx++; |
|
6506 } |
|
6507 ja++; |
|
6508 ja_lt_max= ja < ja_max; |
|
6509 } |
|
6510 else if (( !ja_lt_max ) || |
|
6511 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
6512 { |
|
6513 double tmp = xmin (0., b.data(jb)); |
|
6514 if (tmp != 0.) |
|
6515 { |
|
6516 r.ridx(jx) = b.ridx(jb); |
|
6517 r.data(jx) = tmp; |
|
6518 jx++; |
|
6519 } |
|
6520 jb++; |
|
6521 jb_lt_max= jb < jb_max; |
|
6522 } |
|
6523 else |
|
6524 { |
|
6525 double tmp = xmin (a.data(ja), b.data(jb)); |
|
6526 if (tmp != 0.) |
|
6527 { |
|
6528 r.data(jx) = tmp; |
|
6529 r.ridx(jx) = a.ridx(ja); |
|
6530 jx++; |
|
6531 } |
|
6532 ja++; |
|
6533 ja_lt_max= ja < ja_max; |
|
6534 jb++; |
|
6535 jb_lt_max= jb < jb_max; |
|
6536 } |
|
6537 } |
|
6538 r.cidx(i+1) = jx; |
|
6539 } |
|
6540 |
|
6541 r.maybe_compress (); |
|
6542 } |
|
6543 } |
|
6544 else |
|
6545 (*current_liboctave_error_handler) ("matrix size mismatch"); |
|
6546 |
|
6547 return r; |
|
6548 } |
|
6549 |
|
6550 SparseMatrix |
|
6551 max (double d, const SparseMatrix& m) |
|
6552 { |
|
6553 SparseMatrix result; |
|
6554 |
|
6555 int nr = m.rows (); |
|
6556 int nc = m.columns (); |
|
6557 |
|
6558 EMPTY_RETURN_CHECK (SparseMatrix); |
|
6559 |
|
6560 // Count the number of non-zero elements |
|
6561 if (d > 0.) |
|
6562 { |
|
6563 result = SparseMatrix (nr, nc, d); |
|
6564 for (int j = 0; j < nc; j++) |
|
6565 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6566 { |
|
6567 double tmp = xmax (d, m.data (i)); |
|
6568 |
|
6569 if (tmp != 0.) |
|
6570 { |
|
6571 int idx = m.ridx(i) + j * nr; |
|
6572 result.xdata(idx) = tmp; |
|
6573 result.xridx(idx) = m.ridx(i); |
|
6574 } |
|
6575 } |
|
6576 } |
|
6577 else |
|
6578 { |
|
6579 int nel = 0; |
|
6580 for (int j = 0; j < nc; j++) |
|
6581 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6582 if (xmax (d, m.data (i)) != 0.) |
|
6583 nel++; |
|
6584 |
|
6585 result = SparseMatrix (nr, nc, nel); |
|
6586 |
|
6587 int ii = 0; |
|
6588 result.xcidx(0) = 0; |
|
6589 for (int j = 0; j < nc; j++) |
|
6590 { |
|
6591 for (int i = m.cidx(j); i < m.cidx(j+1); i++) |
|
6592 { |
|
6593 double tmp = xmax (d, m.data (i)); |
|
6594 if (tmp != 0.) |
|
6595 { |
|
6596 result.xdata(ii) = tmp; |
|
6597 result.xridx(ii++) = m.ridx(i); |
|
6598 } |
|
6599 } |
|
6600 result.xcidx(j+1) = ii; |
|
6601 } |
|
6602 } |
|
6603 |
|
6604 return result; |
|
6605 } |
|
6606 |
|
6607 SparseMatrix |
|
6608 max (const SparseMatrix& m, double d) |
|
6609 { |
|
6610 return max (d, m); |
|
6611 } |
|
6612 |
|
6613 SparseMatrix |
|
6614 max (const SparseMatrix& a, const SparseMatrix& b) |
|
6615 { |
|
6616 SparseMatrix r; |
|
6617 |
|
6618 if ((a.rows() == b.rows()) && (a.cols() == b.cols())) |
|
6619 { |
|
6620 int a_nr = a.rows (); |
|
6621 int a_nc = a.cols (); |
|
6622 |
|
6623 int b_nr = b.rows (); |
|
6624 int b_nc = b.cols (); |
|
6625 |
|
6626 if (a_nr != b_nr || a_nc != b_nc) |
|
6627 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc); |
|
6628 else |
|
6629 { |
|
6630 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); |
|
6631 |
|
6632 int jx = 0; |
|
6633 r.cidx (0) = 0; |
|
6634 for (int i = 0 ; i < a_nc ; i++) |
|
6635 { |
|
6636 int ja = a.cidx(i); |
|
6637 int ja_max = a.cidx(i+1); |
|
6638 bool ja_lt_max= ja < ja_max; |
|
6639 |
|
6640 int jb = b.cidx(i); |
|
6641 int jb_max = b.cidx(i+1); |
|
6642 bool jb_lt_max = jb < jb_max; |
|
6643 |
|
6644 while (ja_lt_max || jb_lt_max ) |
|
6645 { |
|
6646 OCTAVE_QUIT; |
|
6647 if ((! jb_lt_max) || |
|
6648 (ja_lt_max && (a.ridx(ja) < b.ridx(jb)))) |
|
6649 { |
|
6650 double tmp = xmax (a.data(ja), 0.); |
|
6651 if (tmp != 0.) |
|
6652 { |
|
6653 r.ridx(jx) = a.ridx(ja); |
|
6654 r.data(jx) = tmp; |
|
6655 jx++; |
|
6656 } |
|
6657 ja++; |
|
6658 ja_lt_max= ja < ja_max; |
|
6659 } |
|
6660 else if (( !ja_lt_max ) || |
|
6661 (jb_lt_max && (b.ridx(jb) < a.ridx(ja)) ) ) |
|
6662 { |
|
6663 double tmp = xmax (0., b.data(jb)); |
|
6664 if (tmp != 0.) |
|
6665 { |
|
6666 r.ridx(jx) = b.ridx(jb); |
|
6667 r.data(jx) = tmp; |
|
6668 jx++; |
|
6669 } |
|
6670 jb++; |
|
6671 jb_lt_max= jb < jb_max; |
|
6672 } |
|
6673 else |
|
6674 { |
|
6675 double tmp = xmax (a.data(ja), b.data(jb)); |
|
6676 if (tmp != 0.) |
|
6677 { |
|
6678 r.data(jx) = tmp; |
|
6679 r.ridx(jx) = a.ridx(ja); |
|
6680 jx++; |
|
6681 } |
|
6682 ja++; |
|
6683 ja_lt_max= ja < ja_max; |
|
6684 jb++; |
|
6685 jb_lt_max= jb < jb_max; |
|
6686 } |
|
6687 } |
|
6688 r.cidx(i+1) = jx; |
|
6689 } |
|
6690 |
|
6691 r.maybe_compress (); |
|
6692 } |
|
6693 } |
|
6694 else |
|
6695 (*current_liboctave_error_handler) ("matrix size mismatch"); |
|
6696 |
|
6697 return r; |
|
6698 } |
|
6699 |
|
6700 SPARSE_SMS_CMP_OPS (SparseMatrix, 0.0, , double, 0.0, ) |
|
6701 SPARSE_SMS_BOOL_OPS (SparseMatrix, double, 0.0) |
|
6702 |
|
6703 SPARSE_SSM_CMP_OPS (double, 0.0, , SparseMatrix, 0.0, ) |
|
6704 SPARSE_SSM_BOOL_OPS (double, SparseMatrix, 0.0) |
|
6705 |
|
6706 SPARSE_SMSM_CMP_OPS (SparseMatrix, 0.0, , SparseMatrix, 0.0, ) |
|
6707 SPARSE_SMSM_BOOL_OPS (SparseMatrix, SparseMatrix, 0.0) |
|
6708 |
|
6709 /* |
|
6710 ;;; Local Variables: *** |
|
6711 ;;; mode: C++ *** |
|
6712 ;;; End: *** |
|
6713 */ |