GNU Octave 11.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
dSparse.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2026 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include <cmath>
31#include <istream>
32#include <limits>
33#include <ostream>
34
35#include "dRowVector.h"
36#include "lapack-proto.h"
37#include "lo-ieee.h"
38#include "mappers.h"
39#include "oct-locbuf.h"
40#include "quit.h"
41
42#include "dDiagMatrix.h"
43#include "CSparse.h"
44#include "boolSparse.h"
45#include "dSparse.h"
46#include "oct-spparms.h"
47#include "sparse-lu.h"
48#include "MatrixType.h"
49#include "oct-sparse.h"
50#include "sparse-util.h"
51#include "sparse-chol.h"
52#include "sparse-qr.h"
53
54#include "Sparse-op-defs.h"
55
56#include "Sparse-diag-op-defs.h"
57
58#include "Sparse-perm-op-defs.h"
59
60#include "sparse-dmsolve.h"
61
63 : MSparse<double> (a.rows (), a.cols (), a.nnz ())
64{
65 octave_idx_type nc = cols ();
66 octave_idx_type nz = a.nnz ();
67
68 for (octave_idx_type i = 0; i < nc + 1; i++)
69 cidx (i) = a.cidx (i);
70
71 for (octave_idx_type i = 0; i < nz; i++)
72 {
73 data (i) = a.data (i);
74 ridx (i) = a.ridx (i);
75 }
76}
77
79 : MSparse<double> (a.rows (), a.cols (), a.length ())
80{
81 octave_idx_type j = 0;
82 octave_idx_type l = a.length ();
83 for (octave_idx_type i = 0; i < l; i++)
84 {
85 cidx (i) = j;
86 if (a(i, i) != 0.0)
87 {
88 data (j) = a(i, i);
89 ridx (j) = i;
90 j++;
91 }
92 }
93 for (octave_idx_type i = l; i <= a.cols (); i++)
94 cidx (i) = j;
95}
96
97bool
99{
100 octave_idx_type nr = rows ();
101 octave_idx_type nc = cols ();
102 octave_idx_type nz = nnz ();
103 octave_idx_type nr_a = a.rows ();
104 octave_idx_type nc_a = a.cols ();
105 octave_idx_type nz_a = a.nnz ();
106
107 if (nr != nr_a || nc != nc_a || nz != nz_a)
108 return false;
109
110 for (octave_idx_type i = 0; i < nc + 1; i++)
111 if (cidx (i) != a.cidx (i))
112 return false;
113
114 for (octave_idx_type i = 0; i < nz; i++)
115 if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
116 return false;
117
118 return true;
119}
120
121bool
123{
124 return !(*this == a);
125}
126
127bool
130 octave_idx_type nr = rows ();
131 octave_idx_type nc = cols ();
132
133 if (nr == nc && nr > 0)
134 {
135 for (octave_idx_type j = 0; j < nc; j++)
136 {
137 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
138 {
139 octave_idx_type ri = ridx (i);
140
141 if (ri != j)
142 {
143 bool found = false;
144
145 for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
146 {
147 if (ridx (k) == j)
148 {
149 if (data (i) == data (k))
150 found = true;
151 break;
152 }
153 }
154
155 if (! found)
156 return false;
157 }
158 }
159 }
160
161 return true;
162 }
163
164 return false;
165}
166
170{
171 MSparse<double>::insert (a, r, c);
172 return *this;
173}
174
177{
178 MSparse<double>::insert (a, indx);
179 return *this;
180}
181
183SparseMatrix::max (int dim, bool nanflag, bool realabs) const
184{
185 Array<octave_idx_type> dummy_idx;
186 return max (dummy_idx, dim, nanflag, realabs);
187}
188
191 bool nanflag, bool realabs) const
192{
193 SparseMatrix result;
194 const dim_vector& dv = dims ();
195 octave_idx_type nr = dv(0);
196 octave_idx_type nc = dv(1);
197
198 if (dim >= dv.ndims ())
199 {
200 idx_arg.resize (dim_vector (nr, nc), 0);
201 return *this;
202 }
203
204 if (dim < 0)
205 dim = dv.first_non_singleton ();
206
207 if (dim == 0)
208 {
209 idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
210
211 if (nr == 0 || nc == 0 || dim >= dv.ndims ())
212 return SparseMatrix (nr == 0 ? 0 : 1, nc);
213
214 octave_idx_type nel = 0;
215 if (realabs)
216 {
217 for (octave_idx_type j = 0; j < nc; j++)
218 {
219 double tmp_max = octave::numeric_limits<double>::NaN ();
220 octave_idx_type idx_j = 0;
221 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
222 {
223 if (ridx (i) != idx_j)
224 break;
225 else
226 idx_j++;
227 }
228
229 if (idx_j != nr)
230 tmp_max = 0.;
231
232 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
233 {
234 double tmp = data (i);
235
236 if (octave::math::isnan (tmp) && ! nanflag)
237 {
238 idx_j = ridx (i);
239 tmp_max = tmp;
240 break;
241 }
242 else if (octave::math::isnan (tmp))
243 continue;
244 else if (octave::math::isnan (tmp_max) || tmp > tmp_max)
245 {
246 idx_j = ridx (i);
247 tmp_max = tmp;
248 }
249
250 }
251
252 if (octave::math::isnan (tmp_max) && ! nanflag)
253 idx_arg.elem (j) = idx_j;
254 else
255 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
256 if (tmp_max != 0.)
257 nel++;
258 }
259 }
260 else
261 {
262 for (octave_idx_type j = 0; j < nc; j++)
263 {
264 double tmp_max = octave::numeric_limits<double>::NaN ();
265 double abs_max = octave::numeric_limits<double>::NaN ();
266 octave_idx_type idx_j = 0;
267 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
268 {
269 if (ridx (i) != idx_j)
270 break;
271 else
272 idx_j++;
273 }
274
275 if (idx_j != nr)
276 {
277 tmp_max = 0.;
278 abs_max = 0.;
279 }
280 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
281 {
282 double tmp = data (i);
283
284 if (octave::math::isnan (tmp) && ! nanflag)
285 {
286 idx_j = ridx (i);
287 tmp_max = tmp;
288 break;
289 }
290 else if (octave::math::isnan (tmp))
291 continue;
292
293 double abs_tmp = std::abs (tmp);
294 if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
295 {
296 idx_j = ridx (i);
297 tmp_max = tmp;
298 abs_max = abs_tmp;
299 }
300 else if (abs_tmp == abs_max && tmp > tmp_max)
301 {
302 idx_j = ridx (i);
303 tmp_max = tmp;
304 abs_max = abs_tmp;
305 }
306
307 }
308
309 if (octave::math::isnan (tmp_max) && ! nanflag)
310 idx_arg.elem (j) = idx_j;
311 else
312 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
313 if (tmp_max != 0.)
314 nel++;
315 }
316 }
317
318 result = SparseMatrix (1, nc, nel);
319
320 octave_idx_type ii = 0;
321 result.xcidx (0) = 0;
322 for (octave_idx_type j = 0; j < nc; j++)
323 {
324 double tmp = elem (idx_arg(j), j);
325 if (tmp != 0.)
326 {
327 result.xdata (ii) = tmp;
328 result.xridx (ii++) = 0;
329 }
330 result.xcidx (j+1) = ii;
331
332 }
333 }
334 else
335 {
336 idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
337
338 if (nr == 0 || nc == 0 || dim >= dv.ndims ())
339 return SparseMatrix (nr, nc == 0 ? 0 : 1);
340
342
343 for (octave_idx_type i = 0; i < nr; i++)
344 found[i] = 0;
345
346 for (octave_idx_type j = 0; j < nc; j++)
347 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
348 if (found[ridx (i)] == -j)
349 found[ridx (i)] = -j - 1;
350
351 for (octave_idx_type i = 0; i < nr; i++)
352 if (found[i] > -nc && found[i] < 0)
353 idx_arg.elem (i) = -found[i];
354
355 if (realabs)
356 {
357 for (octave_idx_type j = 0; j < nc; j++)
358 {
359 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
360 {
361 octave_idx_type ir = ridx (i);
362 octave_idx_type ix = idx_arg.elem (ir);
363 double tmp = data (i);
364
365 if (octave::math::isnan (tmp) && ! nanflag)
366 {
367 idx_arg.elem (ir) = j;
368 break;
369 }
370 else if (octave::math::isnan (tmp))
371 continue;
372 else if (ix == -1 || octave::math::isnan (elem (ir, ix)) ||
373 tmp > elem (ir, ix))
374 idx_arg.elem (ir) = j;
375 }
376 }
377 }
378 else
379 {
380 for (octave_idx_type j = 0; j < nc; j++)
381 {
382 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
383 {
384 octave_idx_type ir = ridx (i);
385 octave_idx_type ix = idx_arg.elem (ir);
386 double tmp = data (i);
387
388 if (octave::math::isnan (tmp) && ! nanflag)
389 {
390 idx_arg.elem (ir) = j;
391 break;
392 }
393 else if (octave::math::isnan (tmp))
394 continue;
395 else if (ix == -1 || octave::math::isnan (elem (ir, ix)) ||
396 std::abs (tmp) > std::abs (elem (ir, ix)))
397 idx_arg.elem (ir) = j;
398 else if (std::abs (tmp) == std::abs (elem (ir, ix)) &&
399 tmp > elem (ir, ix))
400 idx_arg.elem (ir) = j;
401 }
402 }
403 }
404
405 octave_idx_type nel = 0;
406 for (octave_idx_type j = 0; j < nr; j++)
407 if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
408 nel++;
409
410 result = SparseMatrix (nr, 1, nel);
411
412 octave_idx_type ii = 0;
413 result.xcidx (0) = 0;
414 result.xcidx (1) = nel;
415 for (octave_idx_type j = 0; j < nr; j++)
416 {
417 if (idx_arg(j) == -1)
418 {
419 idx_arg(j) = 0;
420 result.xdata (ii) = octave::numeric_limits<double>::NaN ();
421 result.xridx (ii++) = j;
422 }
423 else
424 {
425 double tmp = elem (j, idx_arg(j));
426 if (tmp != 0.)
427 {
428 result.xdata (ii) = tmp;
429 result.xridx (ii++) = j;
430 }
431 }
432 }
433 }
434
435 return result;
436}
437
439SparseMatrix::min (int dim, bool nanflag, bool realabs) const
440{
441 Array<octave_idx_type> dummy_idx;
442 return min (dummy_idx, dim, nanflag, realabs);
443}
444
447 bool nanflag, bool realabs) const
448{
449 SparseMatrix result;
450 const dim_vector& dv = dims ();
451 octave_idx_type nr = dv(0);
452 octave_idx_type nc = dv(1);
453
454 if (dim >= dv.ndims ())
455 {
456 idx_arg.resize (dim_vector (nr, nc), 0);
457 return *this;
458 }
459
460 if (dim < 0)
461 dim = dv.first_non_singleton ();
462
463 if (dim == 0)
464 {
465 idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
466
467 if (nr == 0 || nc == 0 || dim >= dv.ndims ())
468 return SparseMatrix (nr == 0 ? 0 : 1, nc);
469
470 octave_idx_type nel = 0;
471 if (realabs)
472 {
473 for (octave_idx_type j = 0; j < nc; j++)
474 {
475 double tmp_max = octave::numeric_limits<double>::NaN ();
476 octave_idx_type idx_j = 0;
477 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
478 {
479 if (ridx (i) != idx_j)
480 break;
481 else
482 idx_j++;
483 }
484
485 if (idx_j != nr)
486 tmp_max = 0.;
487
488 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
489 {
490 double tmp = data (i);
491
492 if (octave::math::isnan (tmp) && ! nanflag)
493 {
494 idx_j = ridx (i);
495 tmp_max = tmp;
496 break;
497 }
498 else if (octave::math::isnan (tmp))
499 continue;
500 else if (octave::math::isnan (tmp_max) || tmp < tmp_max)
501 {
502 idx_j = ridx (i);
503 tmp_max = tmp;
504 }
505
506 }
507
508 if (octave::math::isnan (tmp_max) && ! nanflag)
509 idx_arg.elem (j) = idx_j;
510 else
511 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
512 if (tmp_max != 0.)
513 nel++;
514 }
515 }
516 else
517 {
518 for (octave_idx_type j = 0; j < nc; j++)
519 {
520 double tmp_max = octave::numeric_limits<double>::NaN ();
521 double abs_max = octave::numeric_limits<double>::NaN ();
522 octave_idx_type idx_j = 0;
523 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
524 {
525 if (ridx (i) != idx_j)
526 break;
527 else
528 idx_j++;
529 }
530
531 if (idx_j != nr)
532 {
533 tmp_max = 0.;
534 abs_max = 0.;
535 }
536 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
537 {
538 double tmp = data (i);
539
540 if (octave::math::isnan (tmp) && ! nanflag)
541 {
542 idx_j = ridx (i);
543 tmp_max = tmp;
544 break;
545 }
546 else if (octave::math::isnan (tmp))
547 continue;
548
549 double abs_tmp = std::abs (tmp);
550 if (octave::math::isnan (abs_max) || abs_tmp < abs_max)
551 {
552 idx_j = ridx (i);
553 tmp_max = tmp;
554 abs_max = abs_tmp;
555 }
556 else if (abs_tmp == abs_max && tmp < tmp_max)
557 {
558 idx_j = ridx (i);
559 tmp_max = tmp;
560 abs_max = abs_tmp;
561 }
562
563 }
564
565 if (octave::math::isnan (tmp_max) && ! nanflag)
566 idx_arg.elem (j) = idx_j;
567 else
568 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
569 if (tmp_max != 0.)
570 nel++;
571 }
572 }
573
574 result = SparseMatrix (1, nc, nel);
575
576 octave_idx_type ii = 0;
577 result.xcidx (0) = 0;
578 for (octave_idx_type j = 0; j < nc; j++)
579 {
580 double tmp = elem (idx_arg(j), j);
581 if (tmp != 0.)
582 {
583 result.xdata (ii) = tmp;
584 result.xridx (ii++) = 0;
585 }
586 result.xcidx (j+1) = ii;
587
588 }
589 }
590 else
591 {
592 idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
593
594 if (nr == 0 || nc == 0 || dim >= dv.ndims ())
595 return SparseMatrix (nr, nc == 0 ? 0 : 1);
596
598
599 for (octave_idx_type i = 0; i < nr; i++)
600 found[i] = 0;
601
602 for (octave_idx_type j = 0; j < nc; j++)
603 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
604 if (found[ridx (i)] == -j)
605 found[ridx (i)] = -j - 1;
606
607 for (octave_idx_type i = 0; i < nr; i++)
608 if (found[i] > -nc && found[i] < 0)
609 idx_arg.elem (i) = -found[i];
610
611 if (realabs)
612 {
613 for (octave_idx_type j = 0; j < nc; j++)
614 {
615 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
616 {
617 octave_idx_type ir = ridx (i);
618 octave_idx_type ix = idx_arg.elem (ir);
619 double tmp = data (i);
620
621 if (octave::math::isnan (tmp) && ! nanflag)
622 {
623 idx_arg.elem (ir) = j;
624 break;
625 }
626 else if (octave::math::isnan (tmp))
627 continue;
628 else if (ix == -1 || octave::math::isnan (elem (ir, ix)) ||
629 tmp < elem (ir, ix))
630 idx_arg.elem (ir) = j;
631 }
632 }
633 }
634 else
635 {
636 for (octave_idx_type j = 0; j < nc; j++)
637 {
638 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
639 {
640 octave_idx_type ir = ridx (i);
641 octave_idx_type ix = idx_arg.elem (ir);
642 double tmp = data (i);
643
644 if (octave::math::isnan (tmp) && ! nanflag)
645 {
646 idx_arg.elem (ir) = j;
647 break;
648 }
649 else if (octave::math::isnan (tmp))
650 continue;
651 else if (ix == -1 || octave::math::isnan (elem (ir, ix)) ||
652 std::abs (tmp) < std::abs (elem (ir, ix)))
653 idx_arg.elem (ir) = j;
654 else if (std::abs (tmp) == std::abs (elem (ir, ix)) &&
655 tmp < elem (ir, ix))
656 idx_arg.elem (ir) = j;
657 }
658 }
659 }
660
661 octave_idx_type nel = 0;
662 for (octave_idx_type j = 0; j < nr; j++)
663 if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
664 nel++;
665
666 result = SparseMatrix (nr, 1, nel);
667
668 octave_idx_type ii = 0;
669 result.xcidx (0) = 0;
670 result.xcidx (1) = nel;
671 for (octave_idx_type j = 0; j < nr; j++)
672 {
673 if (idx_arg(j) == -1)
674 {
675 idx_arg(j) = 0;
676 result.xdata (ii) = octave::numeric_limits<double>::NaN ();
677 result.xridx (ii++) = j;
678 }
679 else
680 {
681 double tmp = elem (j, idx_arg(j));
682 if (tmp != 0.)
683 {
684 result.xdata (ii) = tmp;
685 result.xridx (ii++) = j;
686 }
687 }
688 }
689 }
690
691 return result;
692}
693
694/*
695
696%!assert (max (max (speye (65536))), sparse (1))
697%!assert (min (min (speye (65536))), sparse (0))
698%!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
699%!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
700%!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
701%!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
702%!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
703%!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
704%!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
705%!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
706
707*/
708
711{
712 octave_idx_type nc = columns ();
713 RowVector retval (nc, 0);
714
715 for (octave_idx_type j = 0; j < nc; j++)
716 for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
717 {
718 if (ridx (k) == i)
719 {
720 retval(j) = data (k);
721 break;
722 }
723 }
724
725 return retval;
726}
727
730{
731 octave_idx_type nr = rows ();
732 ColumnVector retval (nr, 0);
733
734 for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
735 retval(ridx (k)) = data (k);
736
737 return retval;
738}
739
743{
744 // Don't use numel to avoid all possibility of an overflow
745 if (rb.rows () > 0 && rb.cols () > 0)
746 insert (rb, ra_idx(0), ra_idx(1));
747 return *this;
748}
749
753{
754 SparseComplexMatrix retval (*this);
755 if (rb.rows () > 0 && rb.cols () > 0)
756 retval.insert (rb, ra_idx(0), ra_idx(1));
757 return retval;
758}
759
762{
763 octave_idx_type nr = a.rows ();
764 octave_idx_type nc = a.cols ();
765 octave_idx_type nz = a.nnz ();
766 SparseMatrix r (nr, nc, nz);
767
768 for (octave_idx_type i = 0; i < nc +1; i++)
769 r.cidx (i) = a.cidx (i);
770
771 for (octave_idx_type i = 0; i < nz; i++)
772 {
773 r.data (i) = std::real (a.data (i));
774 r.ridx (i) = a.ridx (i);
775 }
776
777 r.maybe_compress (true);
778 return r;
779}
780
783{
784 octave_idx_type nr = a.rows ();
785 octave_idx_type nc = a.cols ();
786 octave_idx_type nz = a.nnz ();
787 SparseMatrix r (nr, nc, nz);
788
789 for (octave_idx_type i = 0; i < nc +1; i++)
790 r.cidx (i) = a.cidx (i);
791
792 for (octave_idx_type i = 0; i < nz; i++)
793 {
794 r.data (i) = std::imag (a.data (i));
795 r.ridx (i) = a.ridx (i);
796 }
797
798 r.maybe_compress (true);
799 return r;
800}
801
802/*
803
804%!assert (nnz (real (sparse ([1i,1]))), 1)
805%!assert (nnz (real (sparse ([1i,1]))), 1)
806
807*/
808
809// Local function to check if matrix is singular based on rcond.
810static inline
811bool
812is_singular (const double rcond)
813{
814 return (std::abs (rcond) <= std::numeric_limits<double>::epsilon ());
815}
816
819{
820 octave_idx_type info;
821 double rcond;
822 MatrixType mattype (*this);
823 return inverse (mattype, info, rcond, 0, 0);
824}
825
828{
829 octave_idx_type info;
830 double rcond;
831 return inverse (mattype, info, rcond, 0, 0);
832}
833
836{
837 double rcond;
838 return inverse (mattype, info, rcond, 0, 0);
839}
840
842SparseMatrix::dinverse (MatrixType& mattype, octave_idx_type& info,
843 double& rcond, const bool,
844 const bool calccond) const
845{
846 SparseMatrix retval;
847
848 octave_idx_type nr = rows ();
849 octave_idx_type nc = cols ();
850 info = 0;
851
852 if (nr == 0 || nc == 0 || nr != nc)
853 (*current_liboctave_error_handler) ("inverse requires square matrix");
854
855 // Print spparms("spumoni") info if requested
856 int typ = mattype.type ();
857 mattype.info ();
858
860 (*current_liboctave_error_handler) ("incorrect matrix type");
861
863 retval = transpose ();
864 else
865 retval = *this;
866
867 // Force make_unique to be called
868 double *v = retval.data ();
869
870 if (calccond)
871 {
872 double dmax = 0.;
873 double dmin = octave::numeric_limits<double>::Inf ();
874 for (octave_idx_type i = 0; i < nr; i++)
875 {
876 double tmp = fabs (v[i]);
877 if (tmp > dmax)
878 dmax = tmp;
879 if (tmp < dmin)
880 dmin = tmp;
881 }
882 rcond = dmin / dmax;
883 }
884
885 for (octave_idx_type i = 0; i < nr; i++)
886 v[i] = 1.0 / v[i];
887
888 return retval;
889}
890
892SparseMatrix::tinverse (MatrixType& mattype, octave_idx_type& info,
893 double& rcond, const bool,
894 const bool calccond) const
895{
896 SparseMatrix retval;
897
898 octave_idx_type nr = rows ();
899 octave_idx_type nc = cols ();
900 info = 0;
901
902 if (nr == 0 || nc == 0 || nr != nc)
903 (*current_liboctave_error_handler) ("inverse requires square matrix");
904
905 // Print spparms("spumoni") info if requested
906 int typ = mattype.type ();
907 mattype.info ();
908
911 (*current_liboctave_error_handler) ("incorrect matrix type");
912
913 double anorm = 0.;
914 double ainvnorm = 0.;
915
916 if (calccond)
917 {
918 // Calculate the 1-norm of matrix for rcond calculation
919 for (octave_idx_type j = 0; j < nr; j++)
920 {
921 double atmp = 0.;
922 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
923 atmp += fabs (data (i));
924 if (atmp > anorm)
925 anorm = atmp;
926 }
927 }
928
929 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
930 {
931 octave_idx_type nz = nnz ();
932 octave_idx_type cx = 0;
933 octave_idx_type nz2 = nz;
934 retval = SparseMatrix (nr, nc, nz2);
935
936 for (octave_idx_type i = 0; i < nr; i++)
937 {
938 octave_quit ();
939 // place the 1 in the identity position
940 octave_idx_type cx_colstart = cx;
941
942 if (cx == nz2)
943 {
944 nz2 *= 2;
945 retval.change_capacity (nz2);
946 }
947
948 retval.xcidx (i) = cx;
949 retval.xridx (cx) = i;
950 retval.xdata (cx) = 1.0;
951 cx++;
952
953 // iterate across columns of input matrix
954 for (octave_idx_type j = i+1; j < nr; j++)
955 {
956 double v = 0.;
957 // iterate to calculate sum
958 octave_idx_type colXp = retval.xcidx (i);
959 octave_idx_type colUp = cidx (j);
960 octave_idx_type rpX, rpU;
961
962 if (cidx (j) == cidx (j+1))
963 (*current_liboctave_error_handler) ("division by zero");
964
965 do
966 {
967 octave_quit ();
968 rpX = retval.xridx (colXp);
969 rpU = ridx (colUp);
970
971 if (rpX < rpU)
972 colXp++;
973 else if (rpX > rpU)
974 colUp++;
975 else
976 {
977 v -= retval.xdata (colXp) * data (colUp);
978 colXp++;
979 colUp++;
980 }
981 }
982 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
983
984 // get A(m,m)
985 if (typ == MatrixType::Upper)
986 colUp = cidx (j+1) - 1;
987 else
988 colUp = cidx (j);
989 double pivot = data (colUp);
990 if (pivot == 0. || ridx (colUp) != j)
991 (*current_liboctave_error_handler) ("division by zero");
992
993 if (v != 0.)
994 {
995 if (cx == nz2)
996 {
997 nz2 *= 2;
998 retval.change_capacity (nz2);
999 }
1000
1001 retval.xridx (cx) = j;
1002 retval.xdata (cx) = v / pivot;
1003 cx++;
1004 }
1005 }
1006
1007 // get A(m,m)
1008 octave_idx_type colUp;
1009 if (typ == MatrixType::Upper)
1010 colUp = cidx (i+1) - 1;
1011 else
1012 colUp = cidx (i);
1013 double pivot = data (colUp);
1014 if (pivot == 0. || ridx (colUp) != i)
1015 (*current_liboctave_error_handler) ("division by zero");
1016
1017 if (pivot != 1.0)
1018 for (octave_idx_type j = cx_colstart; j < cx; j++)
1019 retval.xdata (j) /= pivot;
1020 }
1021 retval.xcidx (nr) = cx;
1022 retval.maybe_compress ();
1023 }
1024 else
1025 {
1026 octave_idx_type nz = nnz ();
1027 octave_idx_type cx = 0;
1028 octave_idx_type nz2 = nz;
1029 retval = SparseMatrix (nr, nc, nz2);
1030
1031 OCTAVE_LOCAL_BUFFER (double, work, nr);
1033
1034 octave_idx_type *perm = mattype.triangular_perm ();
1035 if (typ == MatrixType::Permuted_Upper)
1036 {
1037 for (octave_idx_type i = 0; i < nr; i++)
1038 rperm[perm[i]] = i;
1039 }
1040 else
1041 {
1042 for (octave_idx_type i = 0; i < nr; i++)
1043 rperm[i] = perm[i];
1044 for (octave_idx_type i = 0; i < nr; i++)
1045 perm[rperm[i]] = i;
1046 }
1047
1048 for (octave_idx_type i = 0; i < nr; i++)
1049 {
1050 octave_quit ();
1051 octave_idx_type iidx = rperm[i];
1052
1053 for (octave_idx_type j = 0; j < nr; j++)
1054 work[j] = 0.;
1055
1056 // place the 1 in the identity position
1057 work[iidx] = 1.0;
1058
1059 // iterate across columns of input matrix
1060 for (octave_idx_type j = iidx+1; j < nr; j++)
1061 {
1062 double v = 0.;
1063 octave_idx_type jidx = perm[j];
1064 // iterate to calculate sum
1065 for (octave_idx_type k = cidx (jidx);
1066 k < cidx (jidx+1); k++)
1067 {
1068 octave_quit ();
1069 v -= work[ridx (k)] * data (k);
1070 }
1071
1072 // get A(m,m)
1073 double pivot;
1074 if (typ == MatrixType::Permuted_Upper)
1075 pivot = data (cidx (jidx+1) - 1);
1076 else
1077 pivot = data (cidx (jidx));
1078 if (pivot == 0.)
1079 (*current_liboctave_error_handler) ("division by zero");
1080
1081 work[j] = v / pivot;
1082 }
1083
1084 // get A(m,m)
1085 octave_idx_type colUp;
1086 if (typ == MatrixType::Permuted_Upper)
1087 colUp = cidx (perm[iidx]+1) - 1;
1088 else
1089 colUp = cidx (perm[iidx]);
1090
1091 double pivot = data (colUp);
1092 if (pivot == 0.)
1093 (*current_liboctave_error_handler) ("division by zero");
1094
1095 octave_idx_type new_cx = cx;
1096 for (octave_idx_type j = iidx; j < nr; j++)
1097 if (work[j] != 0.0)
1098 {
1099 new_cx++;
1100 if (pivot != 1.0)
1101 work[j] /= pivot;
1102 }
1103
1104 if (cx < new_cx)
1105 {
1106 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1107 retval.change_capacity (nz2);
1108 }
1109
1110 retval.xcidx (i) = cx;
1111 for (octave_idx_type j = iidx; j < nr; j++)
1112 if (work[j] != 0.)
1113 {
1114 retval.xridx (cx) = j;
1115 retval.xdata (cx++) = work[j];
1116 }
1117 }
1118
1119 retval.xcidx (nr) = cx;
1120 retval.maybe_compress ();
1121 }
1122
1123 if (calccond)
1124 {
1125 // Calculate the 1-norm of inverse matrix for rcond calculation
1126 for (octave_idx_type j = 0; j < nr; j++)
1127 {
1128 double atmp = 0.;
1129 for (octave_idx_type i = retval.cidx (j);
1130 i < retval.cidx (j+1); i++)
1131 atmp += fabs (retval.data (i));
1132 if (atmp > ainvnorm)
1133 ainvnorm = atmp;
1134 }
1135
1136 rcond = 1. / ainvnorm / anorm;
1137 }
1138
1139 return retval;
1140}
1141
1144 double& rcond, bool, bool calc_cond) const
1145{
1146 if (nnz () == 0)
1147 {
1148 (*current_liboctave_error_handler)
1149 ("inverse of the null matrix not defined");
1150 }
1151
1152 int typ = mattype.type (false);
1153 SparseMatrix ret;
1154
1155 if (typ == MatrixType::Unknown)
1156 typ = mattype.type (*this);
1157
1159 ret = dinverse (mattype, info, rcond, true, calc_cond);
1160 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1161 ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
1162 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1163 {
1164 MatrixType newtype = mattype.transpose ();
1165 ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
1166 }
1167 else
1168 {
1169 if (mattype.ishermitian ())
1170 {
1171 MatrixType tmp_typ (MatrixType::Upper);
1172 octave::math::sparse_chol<SparseMatrix> fact (*this, info, false);
1173 rcond = fact.rcond ();
1174 if (info == 0)
1175 {
1176 double rcond2;
1177 SparseMatrix Q = fact.Q ();
1178 SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
1179 info, rcond2, true, false);
1180 ret = Q * InvL.transpose () * InvL * Q.transpose ();
1181 }
1182 else
1183 {
1184 // Matrix is either singular or not positive definite
1185 mattype.mark_as_unsymmetric ();
1186 }
1187 }
1188
1189 if (! mattype.ishermitian ())
1190 {
1191 octave_idx_type n = rows ();
1192 ColumnVector Qinit(n);
1193 for (octave_idx_type i = 0; i < n; i++)
1194 Qinit(i) = i;
1195
1196 MatrixType tmp_typ (MatrixType::Upper);
1197 octave::math::sparse_lu<SparseMatrix> fact (*this,
1198 Qinit, Matrix (),
1199 false, false);
1200 rcond = fact.rcond ();
1201 if (rcond == 0.0)
1202 {
1203 // Return all Inf matrix with sparsity pattern of input.
1204 octave_idx_type nz = nnz ();
1205 ret = SparseMatrix (rows (), cols (), nz);
1206 std::fill (ret.xdata (), ret.xdata () + nz,
1207 octave::numeric_limits<double>::Inf ());
1208 std::copy_n (ridx (), nz, ret.xridx ());
1209 std::copy_n (cidx (), cols () + 1, ret.xcidx ());
1210
1211 return ret;
1212 }
1213
1214 double rcond2;
1215 SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
1216 info, rcond2, true, false);
1217 SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
1218 true, false).transpose ();
1219 ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1220 }
1221 }
1222
1223 return ret;
1224}
1225
1226DET
1228{
1229 octave_idx_type info;
1230 double rcond;
1231 return determinant (info, rcond, 0);
1232}
1233
1234DET
1236{
1237 double rcond;
1238 return determinant (info, rcond, 0);
1239}
1240
1241DET
1242SparseMatrix::determinant (octave_idx_type& err, double& rcond, bool) const
1243{
1244 DET retval;
1245
1246#if defined (HAVE_UMFPACK)
1247
1248 octave_idx_type nr = rows ();
1249 octave_idx_type nc = cols ();
1250
1251 if (nr == 0 || nc == 0 || nr != nc)
1252 {
1253 retval = DET (1.0);
1254 }
1255 else
1256 {
1257 err = 0;
1258
1259 // Setup the control parameters
1260 Matrix Control (UMFPACK_CONTROL, 1);
1261 double *control = Control.rwdata ();
1262 UMFPACK_DNAME (defaults) (control);
1263
1264 double tmp = octave::sparse_params::get_key ("spumoni");
1265 if (! octave::math::isnan (tmp))
1266 Control (UMFPACK_PRL) = tmp;
1267
1268 tmp = octave::sparse_params::get_key ("piv_tol");
1269 if (! octave::math::isnan (tmp))
1270 {
1271 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1272 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1273 }
1274
1275 // Set whether we are allowed to modify Q or not
1276 tmp = octave::sparse_params::get_key ("autoamd");
1277 if (! octave::math::isnan (tmp))
1278 Control (UMFPACK_FIXQ) = tmp;
1279
1280 // Turn-off UMFPACK scaling for LU
1281 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1282
1283 UMFPACK_DNAME (report_control) (control);
1284
1285 const octave_idx_type *Ap = cidx ();
1286 const octave_idx_type *Ai = ridx ();
1287 const double *Ax = data ();
1288
1289 UMFPACK_DNAME (report_matrix) (nr, nc,
1290 octave::to_suitesparse_intptr (Ap),
1291 octave::to_suitesparse_intptr (Ai),
1292 Ax, 1, control);
1293
1294 void *Symbolic;
1295 Matrix Info (1, UMFPACK_INFO);
1296 double *info = Info.rwdata ();
1297 int status = UMFPACK_DNAME (qsymbolic) (nr, nc,
1298 octave::to_suitesparse_intptr (Ap),
1299 octave::to_suitesparse_intptr (Ai),
1300 Ax, nullptr, &Symbolic, control, info);
1301
1302 if (status < 0)
1303 {
1304 UMFPACK_DNAME (report_status) (control, status);
1305 UMFPACK_DNAME (report_info) (control, info);
1306
1307 UMFPACK_DNAME (free_symbolic) (&Symbolic);
1308
1309 (*current_liboctave_error_handler)
1310 ("SparseMatrix::determinant symbolic factorization failed");
1311 }
1312 else
1313 {
1314 UMFPACK_DNAME (report_symbolic) (Symbolic, control);
1315
1316 void *Numeric;
1317 status = UMFPACK_DNAME (numeric) (octave::to_suitesparse_intptr (Ap),
1318 octave::to_suitesparse_intptr (Ai),
1319 Ax, Symbolic,
1320 &Numeric, control, info);
1321 UMFPACK_DNAME (free_symbolic) (&Symbolic);
1322
1323 rcond = Info (UMFPACK_RCOND);
1324
1325 if (status < 0)
1326 {
1327 UMFPACK_DNAME (report_status) (control, status);
1328 UMFPACK_DNAME (report_info) (control, info);
1329
1330 UMFPACK_DNAME (free_numeric) (&Numeric);
1331 (*current_liboctave_error_handler)
1332 ("SparseMatrix::determinant numeric factorization failed");
1333 }
1334 else
1335 {
1336 UMFPACK_DNAME (report_numeric) (Numeric, control);
1337
1338 double c10, e10;
1339
1340 status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1341 info);
1342
1343 if (status < 0)
1344 {
1345 UMFPACK_DNAME (report_status) (control, status);
1346 UMFPACK_DNAME (report_info) (control, info);
1347
1348 (*current_liboctave_error_handler)
1349 ("SparseMatrix::determinant error calculating determinant");
1350 }
1351 else
1352 retval = DET (c10, e10, 10);
1353
1354 UMFPACK_DNAME (free_numeric) (&Numeric);
1355 }
1356 }
1357 }
1358
1359#else
1360
1361 octave_unused_parameter (err);
1362 octave_unused_parameter (rcond);
1363
1364 (*current_liboctave_error_handler)
1365 ("support for UMFPACK was unavailable or disabled "
1366 "when liboctave was built");
1367
1368#endif
1369
1370 return retval;
1371}
1372
1373Matrix
1374SparseMatrix::dsolve (MatrixType& mattype, const Matrix& b,
1375 octave_idx_type& err,
1376 double& rcond, solve_singularity_handler,
1377 bool calc_cond) const
1378{
1379 Matrix retval;
1380
1381 octave_idx_type nr = rows ();
1382 octave_idx_type nc = cols ();
1383 octave_idx_type nm = (nc < nr ? nc : nr);
1384 err = 0;
1385
1386 if (nr != b.rows ())
1388 ("matrix dimension mismatch solution of linear equations");
1389
1390 if (nr == 0 || nc == 0 || b.cols () == 0)
1391 retval = Matrix (nc, b.cols (), 0.0);
1392 else
1393 {
1394 // Print spparms("spumoni") info if requested
1395 int typ = mattype.type ();
1396 mattype.info ();
1397
1399 (*current_liboctave_error_handler) ("incorrect matrix type");
1400
1401 retval.resize (nc, b.cols (), 0.);
1402 if (typ == MatrixType::Diagonal)
1403 for (octave_idx_type j = 0; j < b.cols (); j++)
1404 for (octave_idx_type i = 0; i < nm; i++)
1405 retval(i, j) = b(i, j) / data (i);
1406 else
1407 for (octave_idx_type j = 0; j < b.cols (); j++)
1408 for (octave_idx_type k = 0; k < nc; k++)
1409 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1410 retval(k, j) = b(ridx (i), j) / data (i);
1411
1412 if (calc_cond)
1413 {
1414 double dmax = 0.;
1415 double dmin = octave::numeric_limits<double>::Inf ();
1416 for (octave_idx_type i = 0; i < nm; i++)
1417 {
1418 double tmp = fabs (data (i));
1419 if (tmp > dmax)
1420 dmax = tmp;
1421 if (tmp < dmin)
1422 dmin = tmp;
1423 }
1424 rcond = dmin / dmax;
1425 }
1426 else
1427 rcond = 1.;
1428 }
1429
1430 return retval;
1431}
1432
1434SparseMatrix::dsolve (MatrixType& mattype, const SparseMatrix& b,
1435 octave_idx_type& err, double& rcond,
1436 solve_singularity_handler, bool calc_cond) const
1437{
1438 SparseMatrix retval;
1439
1440 octave_idx_type nr = rows ();
1441 octave_idx_type nc = cols ();
1442 octave_idx_type nm = (nc < nr ? nc : nr);
1443 err = 0;
1444
1445 if (nr != b.rows ())
1447 ("matrix dimension mismatch solution of linear equations");
1448
1449 if (nr == 0 || nc == 0 || b.cols () == 0)
1450 retval = SparseMatrix (nc, b.cols ());
1451 else
1452 {
1453 // Print spparms("spumoni") info if requested
1454 int typ = mattype.type ();
1455 mattype.info ();
1456
1458 (*current_liboctave_error_handler) ("incorrect matrix type");
1459
1460 octave_idx_type b_nc = b.cols ();
1461 octave_idx_type b_nz = b.nnz ();
1462 retval = SparseMatrix (nc, b_nc, b_nz);
1463
1464 retval.xcidx (0) = 0;
1465 octave_idx_type ii = 0;
1466 if (typ == MatrixType::Diagonal)
1467 for (octave_idx_type j = 0; j < b_nc; j++)
1468 {
1469 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1470 {
1471 if (b.ridx (i) >= nm)
1472 break;
1473 retval.xridx (ii) = b.ridx (i);
1474 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1475 }
1476 retval.xcidx (j+1) = ii;
1477 }
1478 else
1479 for (octave_idx_type j = 0; j < b_nc; j++)
1480 {
1481 for (octave_idx_type l = 0; l < nc; l++)
1482 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1483 {
1484 bool found = false;
1486 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1487 if (ridx (i) == b.ridx (k))
1488 {
1489 found = true;
1490 break;
1491 }
1492 if (found)
1493 {
1494 retval.xridx (ii) = l;
1495 retval.xdata (ii++) = b.data (k) / data (i);
1496 }
1497 }
1498 retval.xcidx (j+1) = ii;
1499 }
1500
1501 if (calc_cond)
1502 {
1503 double dmax = 0.;
1504 double dmin = octave::numeric_limits<double>::Inf ();
1505 for (octave_idx_type i = 0; i < nm; i++)
1506 {
1507 double tmp = fabs (data (i));
1508 if (tmp > dmax)
1509 dmax = tmp;
1510 if (tmp < dmin)
1511 dmin = tmp;
1512 }
1513 rcond = dmin / dmax;
1514 }
1515 else
1516 rcond = 1.;
1517 }
1518
1519 return retval;
1520}
1521
1523SparseMatrix::dsolve (MatrixType& mattype, const ComplexMatrix& b,
1524 octave_idx_type& err, double& rcond,
1525 solve_singularity_handler, bool calc_cond) const
1526{
1527 ComplexMatrix retval;
1528
1529 octave_idx_type nr = rows ();
1530 octave_idx_type nc = cols ();
1531 octave_idx_type nm = (nc < nr ? nc : nr);
1532 err = 0;
1533
1534 if (nr != b.rows ())
1536 ("matrix dimension mismatch solution of linear equations");
1537
1538 if (nr == 0 || nc == 0 || b.cols () == 0)
1539 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1540 else
1541 {
1542 // Print spparms("spumoni") info if requested
1543 int typ = mattype.type ();
1544 mattype.info ();
1545
1547 (*current_liboctave_error_handler) ("incorrect matrix type");
1548
1549 retval.resize (nc, b.cols (), 0);
1550 if (typ == MatrixType::Diagonal)
1551 for (octave_idx_type j = 0; j < b.cols (); j++)
1552 for (octave_idx_type i = 0; i < nm; i++)
1553 retval(i, j) = b(i, j) / data (i);
1554 else
1555 for (octave_idx_type j = 0; j < b.cols (); j++)
1556 for (octave_idx_type k = 0; k < nc; k++)
1557 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1558 retval(k, j) = b(ridx (i), j) / data (i);
1559
1560 if (calc_cond)
1561 {
1562 double dmax = 0.;
1563 double dmin = octave::numeric_limits<double>::Inf ();
1564 for (octave_idx_type i = 0; i < nm; i++)
1565 {
1566 double tmp = fabs (data (i));
1567 if (tmp > dmax)
1568 dmax = tmp;
1569 if (tmp < dmin)
1570 dmin = tmp;
1571 }
1572 rcond = dmin / dmax;
1573 }
1574 else
1575 rcond = 1.;
1576 }
1577
1578 return retval;
1579}
1580
1582SparseMatrix::dsolve (MatrixType& mattype, const SparseComplexMatrix& b,
1583 octave_idx_type& err, double& rcond,
1584 solve_singularity_handler, bool calc_cond) const
1585{
1586 SparseComplexMatrix retval;
1587
1588 octave_idx_type nr = rows ();
1589 octave_idx_type nc = cols ();
1590 octave_idx_type nm = (nc < nr ? nc : nr);
1591 err = 0;
1592
1593 if (nr != b.rows ())
1595 ("matrix dimension mismatch solution of linear equations");
1596
1597 if (nr == 0 || nc == 0 || b.cols () == 0)
1598 retval = SparseComplexMatrix (nc, b.cols ());
1599 else
1600 {
1601 // Print spparms("spumoni") info if requested
1602 int typ = mattype.type ();
1603 mattype.info ();
1604
1606 (*current_liboctave_error_handler) ("incorrect matrix type");
1607
1608 octave_idx_type b_nc = b.cols ();
1609 octave_idx_type b_nz = b.nnz ();
1610 retval = SparseComplexMatrix (nc, b_nc, b_nz);
1611
1612 retval.xcidx (0) = 0;
1613 octave_idx_type ii = 0;
1614 if (typ == MatrixType::Diagonal)
1615 for (octave_idx_type j = 0; j < b.cols (); j++)
1616 {
1617 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1618 {
1619 if (b.ridx (i) >= nm)
1620 break;
1621 retval.xridx (ii) = b.ridx (i);
1622 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1623 }
1624 retval.xcidx (j+1) = ii;
1625 }
1626 else
1627 for (octave_idx_type j = 0; j < b.cols (); j++)
1628 {
1629 for (octave_idx_type l = 0; l < nc; l++)
1630 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1631 {
1632 bool found = false;
1634 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1635 if (ridx (i) == b.ridx (k))
1636 {
1637 found = true;
1638 break;
1639 }
1640 if (found)
1641 {
1642 retval.xridx (ii) = l;
1643 retval.xdata (ii++) = b.data (k) / data (i);
1644 }
1645 }
1646 retval.xcidx (j+1) = ii;
1647 }
1648
1649 if (calc_cond)
1650 {
1651 double dmax = 0.;
1652 double dmin = octave::numeric_limits<double>::Inf ();
1653 for (octave_idx_type i = 0; i < nm; i++)
1654 {
1655 double tmp = fabs (data (i));
1656 if (tmp > dmax)
1657 dmax = tmp;
1658 if (tmp < dmin)
1659 dmin = tmp;
1660 }
1661 rcond = dmin / dmax;
1662 }
1663 else
1664 rcond = 1.;
1665 }
1666
1667 return retval;
1668}
1669
1670Matrix
1671SparseMatrix::utsolve (MatrixType& mattype, const Matrix& b,
1672 octave_idx_type& err, double& rcond,
1673 solve_singularity_handler sing_handler,
1674 bool calc_cond) const
1675{
1676 Matrix retval;
1677
1678 octave_idx_type nr = rows ();
1679 octave_idx_type nc = cols ();
1680 octave_idx_type nm = std::max (nc, nr);
1681 err = 0;
1682
1683 if (nr != b.rows ())
1684 (*current_liboctave_error_handler)
1685 ("matrix dimension mismatch solution of linear equations");
1686
1687 if (nr == 0 || nc == 0 || b.cols () == 0)
1688 retval = Matrix (nc, b.cols (), 0.0);
1689 else
1690 {
1691 // Print spparms("spumoni") info if requested
1692 int typ = mattype.type ();
1693 mattype.info ();
1694
1695 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1696 (*current_liboctave_error_handler) ("incorrect matrix type");
1697
1698 double anorm = 0.;
1699 double ainvnorm = 0.;
1700 octave_idx_type b_nc = b.cols ();
1701 rcond = 1.;
1702
1703 if (calc_cond)
1704 {
1705 // Calculate the 1-norm of matrix for rcond calculation
1706 for (octave_idx_type j = 0; j < nc; j++)
1707 {
1708 double atmp = 0.;
1709 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1710 atmp += fabs (data (i));
1711 if (atmp > anorm)
1712 anorm = atmp;
1713 }
1714 }
1715
1716 if (typ == MatrixType::Permuted_Upper)
1717 {
1718 retval.resize (nc, b_nc);
1719 octave_idx_type *perm = mattype.triangular_perm ();
1720 OCTAVE_LOCAL_BUFFER (double, work, nm);
1721
1722 for (octave_idx_type j = 0; j < b_nc; j++)
1723 {
1724 for (octave_idx_type i = 0; i < nr; i++)
1725 work[i] = b(i, j);
1726 for (octave_idx_type i = nr; i < nc; i++)
1727 work[i] = 0.;
1728
1729 for (octave_idx_type k = nc-1; k >= 0; k--)
1730 {
1731 octave_idx_type kidx = perm[k];
1732
1733 if (work[k] != 0.)
1734 {
1735 if (ridx (cidx (kidx+1)-1) != k
1736 || data (cidx (kidx+1)-1) == 0.)
1737 {
1738 err = -2;
1739 goto triangular_error;
1740 }
1741
1742 double tmp = work[k] / data (cidx (kidx+1)-1);
1743 work[k] = tmp;
1744 for (octave_idx_type i = cidx (kidx);
1745 i < cidx (kidx+1)-1; i++)
1746 {
1747 octave_idx_type iidx = ridx (i);
1748 work[iidx] = work[iidx] - tmp * data (i);
1749 }
1750 }
1751 }
1752
1753 for (octave_idx_type i = 0; i < nc; i++)
1754 retval.xelem (perm[i], j) = work[i];
1755 }
1756
1757 if (calc_cond)
1758 {
1759 // Calculation of 1-norm of inv(*this)
1760 for (octave_idx_type i = 0; i < nm; i++)
1761 work[i] = 0.;
1762
1763 for (octave_idx_type j = 0; j < nr; j++)
1764 {
1765 work[j] = 1.;
1766
1767 for (octave_idx_type k = j; k >= 0; k--)
1768 {
1769 octave_idx_type iidx = perm[k];
1770
1771 if (work[k] != 0.)
1772 {
1773 double tmp = work[k] / data (cidx (iidx+1)-1);
1774 work[k] = tmp;
1775 for (octave_idx_type i = cidx (iidx);
1776 i < cidx (iidx+1)-1; i++)
1777 {
1778 octave_idx_type idx2 = ridx (i);
1779 work[idx2] = work[idx2] - tmp * data (i);
1780 }
1781 }
1782 }
1783 double atmp = 0;
1784 for (octave_idx_type i = 0; i < j+1; i++)
1785 {
1786 atmp += fabs (work[i]);
1787 work[i] = 0.;
1788 }
1789 if (atmp > ainvnorm)
1790 ainvnorm = atmp;
1791 }
1792 rcond = 1. / ainvnorm / anorm;
1793 }
1794 }
1795 else
1796 {
1797 OCTAVE_LOCAL_BUFFER (double, work, nm);
1798 retval.resize (nc, b_nc);
1799
1800 for (octave_idx_type j = 0; j < b_nc; j++)
1801 {
1802 for (octave_idx_type i = 0; i < nr; i++)
1803 work[i] = b(i, j);
1804 for (octave_idx_type i = nr; i < nc; i++)
1805 work[i] = 0.;
1806
1807 for (octave_idx_type k = nc-1; k >= 0; k--)
1808 {
1809 if (work[k] != 0.)
1810 {
1811 if (ridx (cidx (k+1)-1) != k
1812 || data (cidx (k+1)-1) == 0.)
1813 {
1814 err = -2;
1815 goto triangular_error;
1816 }
1817
1818 double tmp = work[k] / data (cidx (k+1)-1);
1819 work[k] = tmp;
1820 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1821 {
1822 octave_idx_type iidx = ridx (i);
1823 work[iidx] = work[iidx] - tmp * data (i);
1824 }
1825 }
1826 }
1827
1828 for (octave_idx_type i = 0; i < nc; i++)
1829 retval.xelem (i, j) = work[i];
1830 }
1831
1832 if (calc_cond)
1833 {
1834 // Calculation of 1-norm of inv(*this)
1835 for (octave_idx_type i = 0; i < nm; i++)
1836 work[i] = 0.;
1837
1838 for (octave_idx_type j = 0; j < nr; j++)
1839 {
1840 work[j] = 1.;
1841
1842 for (octave_idx_type k = j; k >= 0; k--)
1843 {
1844 if (work[k] != 0.)
1845 {
1846 double tmp = work[k] / data (cidx (k+1)-1);
1847 work[k] = tmp;
1848 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1849 {
1850 octave_idx_type iidx = ridx (i);
1851 work[iidx] = work[iidx] - tmp * data (i);
1852 }
1853 }
1854 }
1855 double atmp = 0;
1856 for (octave_idx_type i = 0; i < j+1; i++)
1857 {
1858 atmp += fabs (work[i]);
1859 work[i] = 0.;
1860 }
1861 if (atmp > ainvnorm)
1862 ainvnorm = atmp;
1863 }
1864 rcond = 1. / ainvnorm / anorm;
1865 }
1866 }
1867
1868 triangular_error:
1869 if (err != 0)
1870 {
1871 if (sing_handler)
1872 {
1873 sing_handler (rcond);
1874 mattype.mark_as_rectangular ();
1875 }
1876 else
1877 octave::warn_singular_matrix (rcond);
1878 }
1879
1880 if (is_singular (rcond) || octave::math::isnan (rcond))
1881 {
1882 err = -2;
1883
1884 if (sing_handler)
1885 {
1886 sing_handler (rcond);
1887 mattype.mark_as_rectangular ();
1888 }
1889 else
1890 octave::warn_singular_matrix (rcond);
1891 }
1892 }
1893
1894 return retval;
1895}
1896
1898SparseMatrix::utsolve (MatrixType& mattype, const SparseMatrix& b,
1899 octave_idx_type& err, double& rcond,
1900 solve_singularity_handler sing_handler,
1901 bool calc_cond) const
1902{
1903 SparseMatrix retval;
1904
1905 octave_idx_type nr = rows ();
1906 octave_idx_type nc = cols ();
1907 octave_idx_type nm = std::max (nc, nr);
1908 err = 0;
1909
1910 if (nr != b.rows ())
1911 (*current_liboctave_error_handler)
1912 ("matrix dimension mismatch solution of linear equations");
1913
1914 if (nr == 0 || nc == 0 || b.cols () == 0)
1915 retval = SparseMatrix (nc, b.cols ());
1916 else
1917 {
1918 // Print spparms("spumoni") info if requested
1919 int typ = mattype.type ();
1920 mattype.info ();
1921
1922 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1923 (*current_liboctave_error_handler) ("incorrect matrix type");
1924
1925 double anorm = 0.;
1926 double ainvnorm = 0.;
1927 rcond = 1.;
1928
1929 if (calc_cond)
1930 {
1931 // Calculate the 1-norm of matrix for rcond calculation
1932 for (octave_idx_type j = 0; j < nc; j++)
1933 {
1934 double atmp = 0.;
1935 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1936 atmp += fabs (data (i));
1937 if (atmp > anorm)
1938 anorm = atmp;
1939 }
1940 }
1941
1942 octave_idx_type b_nc = b.cols ();
1943 octave_idx_type b_nz = b.nnz ();
1944 retval = SparseMatrix (nc, b_nc, b_nz);
1945 retval.xcidx (0) = 0;
1946 octave_idx_type ii = 0;
1947 octave_idx_type x_nz = b_nz;
1948
1949 if (typ == MatrixType::Permuted_Upper)
1950 {
1951 octave_idx_type *perm = mattype.triangular_perm ();
1952 OCTAVE_LOCAL_BUFFER (double, work, nm);
1953
1955 for (octave_idx_type i = 0; i < nc; i++)
1956 rperm[perm[i]] = i;
1957
1958 for (octave_idx_type j = 0; j < b_nc; j++)
1959 {
1960 for (octave_idx_type i = 0; i < nm; i++)
1961 work[i] = 0.;
1962 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1963 work[b.ridx (i)] = b.data (i);
1964
1965 for (octave_idx_type k = nc-1; k >= 0; k--)
1966 {
1967 octave_idx_type kidx = perm[k];
1968
1969 if (work[k] != 0.)
1970 {
1971 if (ridx (cidx (kidx+1)-1) != k
1972 || data (cidx (kidx+1)-1) == 0.)
1973 {
1974 err = -2;
1975 goto triangular_error;
1976 }
1977
1978 double tmp = work[k] / data (cidx (kidx+1)-1);
1979 work[k] = tmp;
1980 for (octave_idx_type i = cidx (kidx);
1981 i < cidx (kidx+1)-1; i++)
1982 {
1983 octave_idx_type iidx = ridx (i);
1984 work[iidx] = work[iidx] - tmp * data (i);
1985 }
1986 }
1987 }
1988
1989 // Count nonzeros in work vector and adjust space in
1990 // retval if needed
1991 octave_idx_type new_nnz = 0;
1992 for (octave_idx_type i = 0; i < nc; i++)
1993 if (work[i] != 0.)
1994 new_nnz++;
1995
1996 if (ii + new_nnz > x_nz)
1997 {
1998 // Resize the sparse matrix
1999 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2000 retval.change_capacity (sz);
2001 x_nz = sz;
2002 }
2003
2004 for (octave_idx_type i = 0; i < nc; i++)
2005 if (work[rperm[i]] != 0.)
2006 {
2007 retval.xridx (ii) = i;
2008 retval.xdata (ii++) = work[rperm[i]];
2009 }
2010 retval.xcidx (j+1) = ii;
2011 }
2012
2013 retval.maybe_compress ();
2014
2015 if (calc_cond)
2016 {
2017 // Calculation of 1-norm of inv(*this)
2018 for (octave_idx_type i = 0; i < nm; i++)
2019 work[i] = 0.;
2020
2021 for (octave_idx_type j = 0; j < nr; j++)
2022 {
2023 work[j] = 1.;
2024
2025 for (octave_idx_type k = j; k >= 0; k--)
2026 {
2027 octave_idx_type iidx = perm[k];
2028
2029 if (work[k] != 0.)
2030 {
2031 double tmp = work[k] / data (cidx (iidx+1)-1);
2032 work[k] = tmp;
2033 for (octave_idx_type i = cidx (iidx);
2034 i < cidx (iidx+1)-1; i++)
2035 {
2036 octave_idx_type idx2 = ridx (i);
2037 work[idx2] = work[idx2] - tmp * data (i);
2038 }
2039 }
2040 }
2041 double atmp = 0;
2042 for (octave_idx_type i = 0; i < j+1; i++)
2043 {
2044 atmp += fabs (work[i]);
2045 work[i] = 0.;
2046 }
2047 if (atmp > ainvnorm)
2048 ainvnorm = atmp;
2049 }
2050 rcond = 1. / ainvnorm / anorm;
2051 }
2052 }
2053 else
2054 {
2055 OCTAVE_LOCAL_BUFFER (double, work, nm);
2056
2057 for (octave_idx_type j = 0; j < b_nc; j++)
2058 {
2059 for (octave_idx_type i = 0; i < nm; i++)
2060 work[i] = 0.;
2061 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2062 work[b.ridx (i)] = b.data (i);
2063
2064 for (octave_idx_type k = nc-1; k >= 0; k--)
2065 {
2066 if (work[k] != 0.)
2067 {
2068 if (ridx (cidx (k+1)-1) != k
2069 || data (cidx (k+1)-1) == 0.)
2070 {
2071 err = -2;
2072 goto triangular_error;
2073 }
2074
2075 double tmp = work[k] / data (cidx (k+1)-1);
2076 work[k] = tmp;
2077 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2078 {
2079 octave_idx_type iidx = ridx (i);
2080 work[iidx] = work[iidx] - tmp * data (i);
2081 }
2082 }
2083 }
2084
2085 // Count nonzeros in work vector and adjust space in
2086 // retval if needed
2087 octave_idx_type new_nnz = 0;
2088 for (octave_idx_type i = 0; i < nc; i++)
2089 if (work[i] != 0.)
2090 new_nnz++;
2091
2092 if (ii + new_nnz > x_nz)
2093 {
2094 // Resize the sparse matrix
2095 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2096 retval.change_capacity (sz);
2097 x_nz = sz;
2098 }
2099
2100 for (octave_idx_type i = 0; i < nc; i++)
2101 if (work[i] != 0.)
2102 {
2103 retval.xridx (ii) = i;
2104 retval.xdata (ii++) = work[i];
2105 }
2106 retval.xcidx (j+1) = ii;
2107 }
2108
2109 retval.maybe_compress ();
2110
2111 if (calc_cond)
2112 {
2113 // Calculation of 1-norm of inv(*this)
2114 for (octave_idx_type i = 0; i < nm; i++)
2115 work[i] = 0.;
2116
2117 for (octave_idx_type j = 0; j < nr; j++)
2118 {
2119 work[j] = 1.;
2120
2121 for (octave_idx_type k = j; k >= 0; k--)
2122 {
2123 if (work[k] != 0.)
2124 {
2125 double tmp = work[k] / data (cidx (k+1)-1);
2126 work[k] = tmp;
2127 for (octave_idx_type i = cidx (k);
2128 i < cidx (k+1)-1; i++)
2129 {
2130 octave_idx_type iidx = ridx (i);
2131 work[iidx] = work[iidx] - tmp * data (i);
2132 }
2133 }
2134 }
2135 double atmp = 0;
2136 for (octave_idx_type i = 0; i < j+1; i++)
2137 {
2138 atmp += fabs (work[i]);
2139 work[i] = 0.;
2140 }
2141 if (atmp > ainvnorm)
2142 ainvnorm = atmp;
2143 }
2144 rcond = 1. / ainvnorm / anorm;
2145 }
2146 }
2147
2148 triangular_error:
2149 if (err != 0)
2150 {
2151 if (sing_handler)
2152 {
2153 sing_handler (rcond);
2154 mattype.mark_as_rectangular ();
2155 }
2156 else
2157 octave::warn_singular_matrix (rcond);
2158 }
2159
2160 if (is_singular (rcond) || octave::math::isnan (rcond))
2161 {
2162 err = -2;
2163
2164 if (sing_handler)
2165 {
2166 sing_handler (rcond);
2167 mattype.mark_as_rectangular ();
2168 }
2169 else
2170 octave::warn_singular_matrix (rcond);
2171 }
2172 }
2173 return retval;
2174}
2175
2177SparseMatrix::utsolve (MatrixType& mattype, const ComplexMatrix& b,
2178 octave_idx_type& err, double& rcond,
2179 solve_singularity_handler sing_handler,
2180 bool calc_cond) const
2181{
2182 ComplexMatrix retval;
2183
2184 octave_idx_type nr = rows ();
2185 octave_idx_type nc = cols ();
2186 octave_idx_type nm = std::max (nc, nr);
2187 err = 0;
2188
2189 if (nr != b.rows ())
2190 (*current_liboctave_error_handler)
2191 ("matrix dimension mismatch solution of linear equations");
2192
2193 if (nr == 0 || nc == 0 || b.cols () == 0)
2194 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2195 else
2196 {
2197 // Print spparms("spumoni") info if requested
2198 int typ = mattype.type ();
2199 mattype.info ();
2200
2201 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2202 (*current_liboctave_error_handler) ("incorrect matrix type");
2203
2204 double anorm = 0.;
2205 double ainvnorm = 0.;
2206 octave_idx_type b_nc = b.cols ();
2207 rcond = 1.;
2208
2209 if (calc_cond)
2210 {
2211 // Calculate the 1-norm of matrix for rcond calculation
2212 for (octave_idx_type j = 0; j < nc; j++)
2213 {
2214 double atmp = 0.;
2215 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2216 atmp += fabs (data (i));
2217 if (atmp > anorm)
2218 anorm = atmp;
2219 }
2220 }
2221
2222 if (typ == MatrixType::Permuted_Upper)
2223 {
2224 retval.resize (nc, b_nc);
2225 octave_idx_type *perm = mattype.triangular_perm ();
2226 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2227
2228 for (octave_idx_type j = 0; j < b_nc; j++)
2229 {
2230 for (octave_idx_type i = 0; i < nr; i++)
2231 cwork[i] = b(i, j);
2232 for (octave_idx_type i = nr; i < nc; i++)
2233 cwork[i] = 0.;
2234
2235 for (octave_idx_type k = nc-1; k >= 0; k--)
2236 {
2237 octave_idx_type kidx = perm[k];
2238
2239 if (cwork[k] != 0.)
2240 {
2241 if (ridx (cidx (kidx+1)-1) != k
2242 || data (cidx (kidx+1)-1) == 0.)
2243 {
2244 err = -2;
2245 goto triangular_error;
2246 }
2247
2248 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2249 cwork[k] = tmp;
2250 for (octave_idx_type i = cidx (kidx);
2251 i < cidx (kidx+1)-1; i++)
2252 {
2253 octave_idx_type iidx = ridx (i);
2254 cwork[iidx] = cwork[iidx] - tmp * data (i);
2255 }
2256 }
2257 }
2258
2259 for (octave_idx_type i = 0; i < nc; i++)
2260 retval.xelem (perm[i], j) = cwork[i];
2261 }
2262
2263 if (calc_cond)
2264 {
2265 // Calculation of 1-norm of inv(*this)
2266 OCTAVE_LOCAL_BUFFER (double, work, nm);
2267 for (octave_idx_type i = 0; i < nm; i++)
2268 work[i] = 0.;
2269
2270 for (octave_idx_type j = 0; j < nr; j++)
2271 {
2272 work[j] = 1.;
2273
2274 for (octave_idx_type k = j; k >= 0; k--)
2275 {
2276 octave_idx_type iidx = perm[k];
2277
2278 if (work[k] != 0.)
2279 {
2280 double tmp = work[k] / data (cidx (iidx+1)-1);
2281 work[k] = tmp;
2282 for (octave_idx_type i = cidx (iidx);
2283 i < cidx (iidx+1)-1; i++)
2284 {
2285 octave_idx_type idx2 = ridx (i);
2286 work[idx2] = work[idx2] - tmp * data (i);
2287 }
2288 }
2289 }
2290 double atmp = 0;
2291 for (octave_idx_type i = 0; i < j+1; i++)
2292 {
2293 atmp += fabs (work[i]);
2294 work[i] = 0.;
2295 }
2296 if (atmp > ainvnorm)
2297 ainvnorm = atmp;
2298 }
2299 rcond = 1. / ainvnorm / anorm;
2300 }
2301 }
2302 else
2303 {
2304 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2305 retval.resize (nc, b_nc);
2306
2307 for (octave_idx_type j = 0; j < b_nc; j++)
2308 {
2309 for (octave_idx_type i = 0; i < nr; i++)
2310 cwork[i] = b(i, j);
2311 for (octave_idx_type i = nr; i < nc; i++)
2312 cwork[i] = 0.;
2313
2314 for (octave_idx_type k = nc-1; k >= 0; k--)
2315 {
2316 if (cwork[k] != 0.)
2317 {
2318 if (ridx (cidx (k+1)-1) != k
2319 || data (cidx (k+1)-1) == 0.)
2320 {
2321 err = -2;
2322 goto triangular_error;
2323 }
2324
2325 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2326 cwork[k] = tmp;
2327 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2328 {
2329 octave_idx_type iidx = ridx (i);
2330 cwork[iidx] = cwork[iidx] - tmp * data (i);
2331 }
2332 }
2333 }
2334
2335 for (octave_idx_type i = 0; i < nc; i++)
2336 retval.xelem (i, j) = cwork[i];
2337 }
2338
2339 if (calc_cond)
2340 {
2341 // Calculation of 1-norm of inv(*this)
2342 OCTAVE_LOCAL_BUFFER (double, work, nm);
2343 for (octave_idx_type i = 0; i < nm; i++)
2344 work[i] = 0.;
2345
2346 for (octave_idx_type j = 0; j < nr; j++)
2347 {
2348 work[j] = 1.;
2349
2350 for (octave_idx_type k = j; k >= 0; k--)
2351 {
2352 if (work[k] != 0.)
2353 {
2354 double tmp = work[k] / data (cidx (k+1)-1);
2355 work[k] = tmp;
2356 for (octave_idx_type i = cidx (k);
2357 i < cidx (k+1)-1; i++)
2358 {
2359 octave_idx_type iidx = ridx (i);
2360 work[iidx] = work[iidx] - tmp * data (i);
2361 }
2362 }
2363 }
2364 double atmp = 0;
2365 for (octave_idx_type i = 0; i < j+1; i++)
2366 {
2367 atmp += fabs (work[i]);
2368 work[i] = 0.;
2369 }
2370 if (atmp > ainvnorm)
2371 ainvnorm = atmp;
2372 }
2373 rcond = 1. / ainvnorm / anorm;
2374 }
2375 }
2376
2377 triangular_error:
2378 if (err != 0)
2379 {
2380 if (sing_handler)
2381 {
2382 sing_handler (rcond);
2383 mattype.mark_as_rectangular ();
2384 }
2385 else
2386 octave::warn_singular_matrix (rcond);
2387 }
2388
2389 if (is_singular (rcond) || octave::math::isnan (rcond))
2390 {
2391 err = -2;
2392
2393 if (sing_handler)
2394 {
2395 sing_handler (rcond);
2396 mattype.mark_as_rectangular ();
2397 }
2398 else
2399 octave::warn_singular_matrix (rcond);
2400 }
2401 }
2402
2403 return retval;
2404}
2405
2407SparseMatrix::utsolve (MatrixType& mattype, const SparseComplexMatrix& b,
2408 octave_idx_type& err, double& rcond,
2409 solve_singularity_handler sing_handler,
2410 bool calc_cond) const
2411{
2412 SparseComplexMatrix retval;
2413
2414 octave_idx_type nr = rows ();
2415 octave_idx_type nc = cols ();
2416 octave_idx_type nm = std::max (nc, nr);
2417 err = 0;
2418
2419 if (nr != b.rows ())
2420 (*current_liboctave_error_handler)
2421 ("matrix dimension mismatch solution of linear equations");
2422
2423 if (nr == 0 || nc == 0 || b.cols () == 0)
2424 retval = SparseComplexMatrix (nc, b.cols ());
2425 else
2426 {
2427 // Print spparms("spumoni") info if requested
2428 int typ = mattype.type ();
2429 mattype.info ();
2430
2431 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2432 (*current_liboctave_error_handler) ("incorrect matrix type");
2433
2434 double anorm = 0.;
2435 double ainvnorm = 0.;
2436 rcond = 1.;
2437
2438 if (calc_cond)
2439 {
2440 // Calculate the 1-norm of matrix for rcond calculation
2441 for (octave_idx_type j = 0; j < nc; j++)
2442 {
2443 double atmp = 0.;
2444 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2445 atmp += fabs (data (i));
2446 if (atmp > anorm)
2447 anorm = atmp;
2448 }
2449 }
2450
2451 octave_idx_type b_nc = b.cols ();
2452 octave_idx_type b_nz = b.nnz ();
2453 retval = SparseComplexMatrix (nc, b_nc, b_nz);
2454 retval.xcidx (0) = 0;
2455 octave_idx_type ii = 0;
2456 octave_idx_type x_nz = b_nz;
2457
2458 if (typ == MatrixType::Permuted_Upper)
2459 {
2460 octave_idx_type *perm = mattype.triangular_perm ();
2461 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2462
2464 for (octave_idx_type i = 0; i < nc; i++)
2465 rperm[perm[i]] = i;
2466
2467 for (octave_idx_type j = 0; j < b_nc; j++)
2468 {
2469 for (octave_idx_type i = 0; i < nm; i++)
2470 cwork[i] = 0.;
2471 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2472 cwork[b.ridx (i)] = b.data (i);
2473
2474 for (octave_idx_type k = nc-1; k >= 0; k--)
2475 {
2476 octave_idx_type kidx = perm[k];
2477
2478 if (cwork[k] != 0.)
2479 {
2480 if (ridx (cidx (kidx+1)-1) != k
2481 || data (cidx (kidx+1)-1) == 0.)
2482 {
2483 err = -2;
2484 goto triangular_error;
2485 }
2486
2487 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2488 cwork[k] = tmp;
2489 for (octave_idx_type i = cidx (kidx);
2490 i < cidx (kidx+1)-1; i++)
2491 {
2492 octave_idx_type iidx = ridx (i);
2493 cwork[iidx] = cwork[iidx] - tmp * data (i);
2494 }
2495 }
2496 }
2497
2498 // Count nonzeros in work vector and adjust space in
2499 // retval if needed
2500 octave_idx_type new_nnz = 0;
2501 for (octave_idx_type i = 0; i < nc; i++)
2502 if (cwork[i] != 0.)
2503 new_nnz++;
2504
2505 if (ii + new_nnz > x_nz)
2506 {
2507 // Resize the sparse matrix
2508 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2509 retval.change_capacity (sz);
2510 x_nz = sz;
2511 }
2512
2513 for (octave_idx_type i = 0; i < nc; i++)
2514 if (cwork[rperm[i]] != 0.)
2515 {
2516 retval.xridx (ii) = i;
2517 retval.xdata (ii++) = cwork[rperm[i]];
2518 }
2519 retval.xcidx (j+1) = ii;
2520 }
2521
2522 retval.maybe_compress ();
2523
2524 if (calc_cond)
2525 {
2526 // Calculation of 1-norm of inv(*this)
2527 OCTAVE_LOCAL_BUFFER (double, work, nm);
2528 for (octave_idx_type i = 0; i < nm; i++)
2529 work[i] = 0.;
2530
2531 for (octave_idx_type j = 0; j < nr; j++)
2532 {
2533 work[j] = 1.;
2534
2535 for (octave_idx_type k = j; k >= 0; k--)
2536 {
2537 octave_idx_type iidx = perm[k];
2538
2539 if (work[k] != 0.)
2540 {
2541 double tmp = work[k] / data (cidx (iidx+1)-1);
2542 work[k] = tmp;
2543 for (octave_idx_type i = cidx (iidx);
2544 i < cidx (iidx+1)-1; i++)
2545 {
2546 octave_idx_type idx2 = ridx (i);
2547 work[idx2] = work[idx2] - tmp * data (i);
2548 }
2549 }
2550 }
2551 double atmp = 0;
2552 for (octave_idx_type i = 0; i < j+1; i++)
2553 {
2554 atmp += fabs (work[i]);
2555 work[i] = 0.;
2556 }
2557 if (atmp > ainvnorm)
2558 ainvnorm = atmp;
2559 }
2560 rcond = 1. / ainvnorm / anorm;
2561 }
2562 }
2563 else
2564 {
2565 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2566
2567 for (octave_idx_type j = 0; j < b_nc; j++)
2568 {
2569 for (octave_idx_type i = 0; i < nm; i++)
2570 cwork[i] = 0.;
2571 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2572 cwork[b.ridx (i)] = b.data (i);
2573
2574 for (octave_idx_type k = nc-1; k >= 0; k--)
2575 {
2576 if (cwork[k] != 0.)
2577 {
2578 if (ridx (cidx (k+1)-1) != k
2579 || data (cidx (k+1)-1) == 0.)
2580 {
2581 err = -2;
2582 goto triangular_error;
2583 }
2584
2585 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2586 cwork[k] = tmp;
2587 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2588 {
2589 octave_idx_type iidx = ridx (i);
2590 cwork[iidx] = cwork[iidx] - tmp * data (i);
2591 }
2592 }
2593 }
2594
2595 // Count nonzeros in work vector and adjust space in
2596 // retval if needed
2597 octave_idx_type new_nnz = 0;
2598 for (octave_idx_type i = 0; i < nc; i++)
2599 if (cwork[i] != 0.)
2600 new_nnz++;
2601
2602 if (ii + new_nnz > x_nz)
2603 {
2604 // Resize the sparse matrix
2605 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2606 retval.change_capacity (sz);
2607 x_nz = sz;
2608 }
2609
2610 for (octave_idx_type i = 0; i < nc; i++)
2611 if (cwork[i] != 0.)
2612 {
2613 retval.xridx (ii) = i;
2614 retval.xdata (ii++) = cwork[i];
2615 }
2616 retval.xcidx (j+1) = ii;
2617 }
2618
2619 retval.maybe_compress ();
2620
2621 if (calc_cond)
2622 {
2623 // Calculation of 1-norm of inv(*this)
2624 OCTAVE_LOCAL_BUFFER (double, work, nm);
2625 for (octave_idx_type i = 0; i < nm; i++)
2626 work[i] = 0.;
2627
2628 for (octave_idx_type j = 0; j < nr; j++)
2629 {
2630 work[j] = 1.;
2631
2632 for (octave_idx_type k = j; k >= 0; k--)
2633 {
2634 if (work[k] != 0.)
2635 {
2636 double tmp = work[k] / data (cidx (k+1)-1);
2637 work[k] = tmp;
2638 for (octave_idx_type i = cidx (k);
2639 i < cidx (k+1)-1; i++)
2640 {
2641 octave_idx_type iidx = ridx (i);
2642 work[iidx] = work[iidx] - tmp * data (i);
2643 }
2644 }
2645 }
2646 double atmp = 0;
2647 for (octave_idx_type i = 0; i < j+1; i++)
2648 {
2649 atmp += fabs (work[i]);
2650 work[i] = 0.;
2651 }
2652 if (atmp > ainvnorm)
2653 ainvnorm = atmp;
2654 }
2655 rcond = 1. / ainvnorm / anorm;
2656 }
2657 }
2658
2659 triangular_error:
2660 if (err != 0)
2661 {
2662 if (sing_handler)
2663 {
2664 sing_handler (rcond);
2665 mattype.mark_as_rectangular ();
2666 }
2667 else
2668 octave::warn_singular_matrix (rcond);
2669 }
2670
2671 if (is_singular (rcond) || octave::math::isnan (rcond))
2672 {
2673 err = -2;
2674
2675 if (sing_handler)
2676 {
2677 sing_handler (rcond);
2678 mattype.mark_as_rectangular ();
2679 }
2680 else
2681 octave::warn_singular_matrix (rcond);
2682 }
2683 }
2684
2685 return retval;
2686}
2687
2688Matrix
2689SparseMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
2690 octave_idx_type& err, double& rcond,
2691 solve_singularity_handler sing_handler,
2692 bool calc_cond) const
2693{
2694 Matrix retval;
2695
2696 octave_idx_type nr = rows ();
2697 octave_idx_type nc = cols ();
2698 octave_idx_type nm = std::max (nc, nr);
2699 err = 0;
2700
2701 if (nr != b.rows ())
2702 (*current_liboctave_error_handler)
2703 ("matrix dimension mismatch solution of linear equations");
2704
2705 if (nr == 0 || nc == 0 || b.cols () == 0)
2706 retval = Matrix (nc, b.cols (), 0.0);
2707 else
2708 {
2709 // Print spparms("spumoni") info if requested
2710 int typ = mattype.type ();
2711 mattype.info ();
2712
2713 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2714 (*current_liboctave_error_handler) ("incorrect matrix type");
2715
2716 double anorm = 0.;
2717 double ainvnorm = 0.;
2718 octave_idx_type b_nc = b.cols ();
2719 rcond = 1.;
2720
2721 if (calc_cond)
2722 {
2723 // Calculate the 1-norm of matrix for rcond calculation
2724 for (octave_idx_type j = 0; j < nc; j++)
2725 {
2726 double atmp = 0.;
2727 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2728 atmp += fabs (data (i));
2729 if (atmp > anorm)
2730 anorm = atmp;
2731 }
2732 }
2733
2734 if (typ == MatrixType::Permuted_Lower)
2735 {
2736 retval.resize (nc, b_nc);
2737 OCTAVE_LOCAL_BUFFER (double, work, nm);
2738
2739 octave_idx_type *perm = mattype.triangular_perm ();
2740 OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
2741 for (octave_idx_type i = 0; i < nr; i++)
2742 inv_perm[perm[i]] = i;
2743
2744 for (octave_idx_type j = 0; j < b_nc; j++)
2745 {
2746 if (nc > nr)
2747 for (octave_idx_type i = 0; i < nm; i++)
2748 work[i] = 0.;
2749 for (octave_idx_type i = 0; i < nr; i++)
2750 work[inv_perm[i]] = b(i, j);
2751
2752 for (octave_idx_type k = 0; k < nc; k++)
2753 {
2754 if (work[k] != 0.)
2755 {
2756 octave_idx_type minr = nr;
2757 octave_idx_type mini = 0;
2758
2759 // Check that permutation leads to lower triangular form
2760 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2761 if (inv_perm[ridx (i)] < minr)
2762 {
2763 minr = inv_perm[ridx (i)];
2764 mini = i;
2765 }
2766
2767 if (minr != k || data (mini) == 0.)
2768 {
2769 err = -2;
2770 goto triangular_error;
2771 }
2772
2773 double tmp = work[k] / data (mini);
2774 work[k] = tmp;
2775 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2776 {
2777 if (i == mini)
2778 continue;
2779
2780 octave_idx_type iidx = inv_perm[ridx (i)];
2781 work[iidx] = work[iidx] - tmp * data (i);
2782 }
2783 }
2784 }
2785
2786 for (octave_idx_type i = 0; i < nc; i++)
2787 retval(i, j) = work[i];
2788 }
2789
2790 if (calc_cond)
2791 {
2792 // Calculation of 1-norm of inv(*this)
2793 for (octave_idx_type i = 0; i < nm; i++)
2794 work[i] = 0.;
2795
2796 for (octave_idx_type j = 0; j < nr; j++)
2797 {
2798 work[j] = 1.;
2799
2800 for (octave_idx_type k = 0; k < nc; k++)
2801 {
2802 if (work[k] != 0.)
2803 {
2804 octave_idx_type minr = nr;
2805 octave_idx_type mini = 0;
2806
2807 for (octave_idx_type i = cidx (k);
2808 i < cidx (k+1); i++)
2809 if (inv_perm[ridx (i)] < minr)
2810 {
2811 minr = inv_perm[ridx (i)];
2812 mini = i;
2813 }
2814
2815 double tmp = work[k] / data (mini);
2816 work[k] = tmp;
2817 for (octave_idx_type i = cidx (k);
2818 i < cidx (k+1); i++)
2819 {
2820 if (i == mini)
2821 continue;
2822
2823 octave_idx_type iidx = inv_perm[ridx (i)];
2824 work[iidx] = work[iidx] - tmp * data (i);
2825 }
2826 }
2827 }
2828
2829 double atmp = 0.;
2830 for (octave_idx_type i = j; i < nc; i++)
2831 {
2832 atmp += fabs (work[i]);
2833 work[i] = 0.;
2834 }
2835 if (atmp > ainvnorm)
2836 ainvnorm = atmp;
2837 }
2838 rcond = 1. / ainvnorm / anorm;
2839 }
2840 }
2841 else
2842 {
2843 OCTAVE_LOCAL_BUFFER (double, work, nm);
2844 retval.resize (nc, b_nc, 0.);
2845
2846 for (octave_idx_type j = 0; j < b_nc; j++)
2847 {
2848 for (octave_idx_type i = 0; i < nr; i++)
2849 work[i] = b(i, j);
2850 for (octave_idx_type i = nr; i < nc; i++)
2851 work[i] = 0.;
2852 for (octave_idx_type k = 0; k < nc; k++)
2853 {
2854 if (work[k] != 0.)
2855 {
2856 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2857 {
2858 err = -2;
2859 goto triangular_error;
2860 }
2861
2862 double tmp = work[k] / data (cidx (k));
2863 work[k] = tmp;
2864 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2865 {
2866 octave_idx_type iidx = ridx (i);
2867 work[iidx] = work[iidx] - tmp * data (i);
2868 }
2869 }
2870 }
2871
2872 for (octave_idx_type i = 0; i < nc; i++)
2873 retval.xelem (i, j) = work[i];
2874 }
2875
2876 if (calc_cond)
2877 {
2878 // Calculation of 1-norm of inv(*this)
2879 for (octave_idx_type i = 0; i < nm; i++)
2880 work[i] = 0.;
2881
2882 for (octave_idx_type j = 0; j < nr; j++)
2883 {
2884 work[j] = 1.;
2885
2886 for (octave_idx_type k = j; k < nc; k++)
2887 {
2888
2889 if (work[k] != 0.)
2890 {
2891 double tmp = work[k] / data (cidx (k));
2892 work[k] = tmp;
2893 for (octave_idx_type i = cidx (k)+1;
2894 i < cidx (k+1); i++)
2895 {
2896 octave_idx_type iidx = ridx (i);
2897 work[iidx] = work[iidx] - tmp * data (i);
2898 }
2899 }
2900 }
2901 double atmp = 0;
2902 for (octave_idx_type i = j; i < nc; i++)
2903 {
2904 atmp += fabs (work[i]);
2905 work[i] = 0.;
2906 }
2907 if (atmp > ainvnorm)
2908 ainvnorm = atmp;
2909 }
2910 rcond = 1. / ainvnorm / anorm;
2911 }
2912 }
2913
2914 triangular_error:
2915 if (err != 0)
2916 {
2917 if (sing_handler)
2918 {
2919 sing_handler (rcond);
2920 mattype.mark_as_rectangular ();
2921 }
2922 else
2923 octave::warn_singular_matrix (rcond);
2924 }
2925
2926 if (is_singular (rcond) || octave::math::isnan (rcond))
2927 {
2928 err = -2;
2929
2930 if (sing_handler)
2931 {
2932 sing_handler (rcond);
2933 mattype.mark_as_rectangular ();
2934 }
2935 else
2936 octave::warn_singular_matrix (rcond);
2937 }
2938 }
2939
2940 return retval;
2941}
2942
2944SparseMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
2945 octave_idx_type& err, double& rcond,
2946 solve_singularity_handler sing_handler,
2947 bool calc_cond) const
2948{
2949 SparseMatrix retval;
2950
2951 octave_idx_type nr = rows ();
2952 octave_idx_type nc = cols ();
2953 octave_idx_type nm = std::max (nc, nr);
2954 err = 0;
2955
2956 if (nr != b.rows ())
2957 (*current_liboctave_error_handler)
2958 ("matrix dimension mismatch solution of linear equations");
2959
2960 if (nr == 0 || nc == 0 || b.cols () == 0)
2961 retval = SparseMatrix (nc, b.cols ());
2962 else
2963 {
2964 // Print spparms("spumoni") info if requested
2965 int typ = mattype.type ();
2966 mattype.info ();
2967
2968 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2969 (*current_liboctave_error_handler) ("incorrect matrix type");
2970
2971 double anorm = 0.;
2972 double ainvnorm = 0.;
2973 rcond = 1.;
2974
2975 if (calc_cond)
2976 {
2977 // Calculate the 1-norm of matrix for rcond calculation
2978 for (octave_idx_type j = 0; j < nc; j++)
2979 {
2980 double atmp = 0.;
2981 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2982 atmp += fabs (data (i));
2983 if (atmp > anorm)
2984 anorm = atmp;
2985 }
2986 }
2987
2988 octave_idx_type b_nc = b.cols ();
2989 octave_idx_type b_nz = b.nnz ();
2990 retval = SparseMatrix (nc, b_nc, b_nz);
2991 retval.xcidx (0) = 0;
2992 octave_idx_type ii = 0;
2993 octave_idx_type x_nz = b_nz;
2994
2995 if (typ == MatrixType::Permuted_Lower)
2996 {
2997 OCTAVE_LOCAL_BUFFER (double, work, nm);
2998
2999 octave_idx_type *perm = mattype.triangular_perm ();
3000 OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
3001 for (octave_idx_type i = 0; i < nr; i++)
3002 inv_perm[perm[i]] = i;
3003
3004 for (octave_idx_type j = 0; j < b_nc; j++)
3005 {
3006 for (octave_idx_type i = 0; i < nm; i++)
3007 work[i] = 0.;
3008 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3009 work[inv_perm[b.ridx (i)]] = b.data (i);
3010
3011 for (octave_idx_type k = 0; k < nc; k++)
3012 {
3013 if (work[k] != 0.)
3014 {
3015 octave_idx_type minr = nr;
3016 octave_idx_type mini = 0;
3017
3018 // Check that permutation leads to lower triangular form
3019 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3020 if (inv_perm[ridx (i)] < minr)
3021 {
3022 minr = inv_perm[ridx (i)];
3023 mini = i;
3024 }
3025
3026 if (minr != k || data (mini) == 0.)
3027 {
3028 err = -2;
3029 goto triangular_error;
3030 }
3031
3032 double tmp = work[k] / data (mini);
3033 work[k] = tmp;
3034 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3035 {
3036 if (i == mini)
3037 continue;
3038
3039 octave_idx_type iidx = inv_perm[ridx (i)];
3040 work[iidx] = work[iidx] - tmp * data (i);
3041 }
3042 }
3043 }
3044
3045 // Count nonzeros in work vector and adjust space in
3046 // retval if needed
3047 octave_idx_type new_nnz = 0;
3048 for (octave_idx_type i = 0; i < nc; i++)
3049 if (work[i] != 0.)
3050 new_nnz++;
3051
3052 if (ii + new_nnz > x_nz)
3053 {
3054 // Resize the sparse matrix
3055 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3056 retval.change_capacity (sz);
3057 x_nz = sz;
3058 }
3059
3060 for (octave_idx_type i = 0; i < nc; i++)
3061 if (work[i] != 0.)
3062 {
3063 retval.xridx (ii) = i;
3064 retval.xdata (ii++) = work[i];
3065 }
3066 retval.xcidx (j+1) = ii;
3067 }
3068
3069 retval.maybe_compress ();
3070
3071 if (calc_cond)
3072 {
3073 // Calculation of 1-norm of inv(*this)
3074 for (octave_idx_type i = 0; i < nm; i++)
3075 work[i] = 0.;
3076
3077 for (octave_idx_type j = 0; j < nr; j++)
3078 {
3079 work[j] = 1.;
3080
3081 for (octave_idx_type k = 0; k < nc; k++)
3082 {
3083 if (work[k] != 0.)
3084 {
3085 octave_idx_type minr = nr;
3086 octave_idx_type mini = 0;
3087
3088 for (octave_idx_type i = cidx (k);
3089 i < cidx (k+1); i++)
3090 if (inv_perm[ridx (i)] < minr)
3091 {
3092 minr = inv_perm[ridx (i)];
3093 mini = i;
3094 }
3095
3096 double tmp = work[k] / data (mini);
3097 work[k] = tmp;
3098 for (octave_idx_type i = cidx (k);
3099 i < cidx (k+1); i++)
3100 {
3101 if (i == mini)
3102 continue;
3103
3104 octave_idx_type iidx = inv_perm[ridx (i)];
3105 work[iidx] = work[iidx] - tmp * data (i);
3106 }
3107 }
3108 }
3109
3110 double atmp = 0;
3111 for (octave_idx_type i = j; i < nc; i++)
3112 {
3113 atmp += fabs (work[i]);
3114 work[i] = 0.;
3115 }
3116 if (atmp > ainvnorm)
3117 ainvnorm = atmp;
3118 }
3119 rcond = 1. / ainvnorm / anorm;
3120 }
3121 }
3122 else
3123 {
3124 OCTAVE_LOCAL_BUFFER (double, work, nm);
3125
3126 for (octave_idx_type j = 0; j < b_nc; j++)
3127 {
3128 for (octave_idx_type i = 0; i < nm; i++)
3129 work[i] = 0.;
3130 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3131 work[b.ridx (i)] = b.data (i);
3132
3133 for (octave_idx_type k = 0; k < nc; k++)
3134 {
3135 if (work[k] != 0.)
3136 {
3137 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3138 {
3139 err = -2;
3140 goto triangular_error;
3141 }
3142
3143 double tmp = work[k] / data (cidx (k));
3144 work[k] = tmp;
3145 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3146 {
3147 octave_idx_type iidx = ridx (i);
3148 work[iidx] = work[iidx] - tmp * data (i);
3149 }
3150 }
3151 }
3152
3153 // Count nonzeros in work vector and adjust space in
3154 // retval if needed
3155 octave_idx_type new_nnz = 0;
3156 for (octave_idx_type i = 0; i < nc; i++)
3157 if (work[i] != 0.)
3158 new_nnz++;
3159
3160 if (ii + new_nnz > x_nz)
3161 {
3162 // Resize the sparse matrix
3163 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3164 retval.change_capacity (sz);
3165 x_nz = sz;
3166 }
3167
3168 for (octave_idx_type i = 0; i < nc; i++)
3169 if (work[i] != 0.)
3170 {
3171 retval.xridx (ii) = i;
3172 retval.xdata (ii++) = work[i];
3173 }
3174 retval.xcidx (j+1) = ii;
3175 }
3176
3177 retval.maybe_compress ();
3178
3179 if (calc_cond)
3180 {
3181 // Calculation of 1-norm of inv(*this)
3182 for (octave_idx_type i = 0; i < nm; i++)
3183 work[i] = 0.;
3184
3185 for (octave_idx_type j = 0; j < nr; j++)
3186 {
3187 work[j] = 1.;
3188
3189 for (octave_idx_type k = j; k < nc; k++)
3190 {
3191
3192 if (work[k] != 0.)
3193 {
3194 double tmp = work[k] / data (cidx (k));
3195 work[k] = tmp;
3196 for (octave_idx_type i = cidx (k)+1;
3197 i < cidx (k+1); i++)
3198 {
3199 octave_idx_type iidx = ridx (i);
3200 work[iidx] = work[iidx] - tmp * data (i);
3201 }
3202 }
3203 }
3204 double atmp = 0;
3205 for (octave_idx_type i = j; i < nc; i++)
3206 {
3207 atmp += fabs (work[i]);
3208 work[i] = 0.;
3209 }
3210 if (atmp > ainvnorm)
3211 ainvnorm = atmp;
3212 }
3213 rcond = 1. / ainvnorm / anorm;
3214 }
3215 }
3216
3217 triangular_error:
3218 if (err != 0)
3219 {
3220 if (sing_handler)
3221 {
3222 sing_handler (rcond);
3223 mattype.mark_as_rectangular ();
3224 }
3225 else
3226 octave::warn_singular_matrix (rcond);
3227 }
3228
3229 if (is_singular (rcond) || octave::math::isnan (rcond))
3230 {
3231 err = -2;
3232
3233 if (sing_handler)
3234 {
3235 sing_handler (rcond);
3236 mattype.mark_as_rectangular ();
3237 }
3238 else
3239 octave::warn_singular_matrix (rcond);
3240 }
3241 }
3242
3243 return retval;
3244}
3245
3247SparseMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
3248 octave_idx_type& err, double& rcond,
3249 solve_singularity_handler sing_handler,
3250 bool calc_cond) const
3251{
3252 ComplexMatrix retval;
3253
3254 octave_idx_type nr = rows ();
3255 octave_idx_type nc = cols ();
3256 octave_idx_type nm = std::max (nc, nr);
3257 err = 0;
3258
3259 if (nr != b.rows ())
3260 (*current_liboctave_error_handler)
3261 ("matrix dimension mismatch solution of linear equations");
3262
3263 if (nr == 0 || nc == 0 || b.cols () == 0)
3264 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3265 else
3266 {
3267 // Print spparms("spumoni") info if requested
3268 int typ = mattype.type ();
3269 mattype.info ();
3270
3271 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3272 (*current_liboctave_error_handler) ("incorrect matrix type");
3273
3274 double anorm = 0.;
3275 double ainvnorm = 0.;
3276 octave_idx_type b_nc = b.cols ();
3277 rcond = 1.;
3278
3279 if (calc_cond)
3280 {
3281 // Calculate the 1-norm of matrix for rcond calculation
3282 for (octave_idx_type j = 0; j < nc; j++)
3283 {
3284 double atmp = 0.;
3285 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3286 atmp += fabs (data (i));
3287 if (atmp > anorm)
3288 anorm = atmp;
3289 }
3290 }
3291
3292 if (typ == MatrixType::Permuted_Lower)
3293 {
3294 retval.resize (nc, b_nc);
3295 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3296
3297 octave_idx_type *perm = mattype.triangular_perm ();
3298 OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
3299 for (octave_idx_type i = 0; i < nr; i++)
3300 inv_perm[perm[i]] = i;
3301
3302 for (octave_idx_type j = 0; j < b_nc; j++)
3303 {
3304 for (octave_idx_type i = 0; i < nm; i++)
3305 cwork[i] = 0.;
3306 for (octave_idx_type i = 0; i < nr; i++)
3307 cwork[inv_perm[i]] = b(i, j);
3308
3309 for (octave_idx_type k = 0; k < nc; k++)
3310 {
3311 if (cwork[k] != 0.)
3312 {
3313 octave_idx_type minr = nr;
3314 octave_idx_type mini = 0;
3315
3316 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3317 if (inv_perm[ridx (i)] < minr)
3318 {
3319 minr = inv_perm[ridx (i)];
3320 mini = i;
3321 }
3322
3323 if (minr != k || data (mini) == 0)
3324 {
3325 err = -2;
3326 goto triangular_error;
3327 }
3328
3329 Complex tmp = cwork[k] / data (mini);
3330 cwork[k] = tmp;
3331 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3332 {
3333 if (i == mini)
3334 continue;
3335
3336 octave_idx_type iidx = inv_perm[ridx (i)];
3337 cwork[iidx] = cwork[iidx] - tmp * data (i);
3338 }
3339 }
3340 }
3341
3342 for (octave_idx_type i = 0; i < nc; i++)
3343 retval(i, j) = cwork[i];
3344 }
3345
3346 if (calc_cond)
3347 {
3348 // Calculation of 1-norm of inv(*this)
3349 OCTAVE_LOCAL_BUFFER (double, work, nm);
3350 for (octave_idx_type i = 0; i < nm; i++)
3351 work[i] = 0.;
3352
3353 for (octave_idx_type j = 0; j < nr; j++)
3354 {
3355 work[j] = 1.;
3356
3357 for (octave_idx_type k = 0; k < nc; k++)
3358 {
3359 if (work[k] != 0.)
3360 {
3361 octave_idx_type minr = nr;
3362 octave_idx_type mini = 0;
3363
3364 for (octave_idx_type i = cidx (k);
3365 i < cidx (k+1); i++)
3366 if (inv_perm[ridx (i)] < minr)
3367 {
3368 minr = inv_perm[ridx (i)];
3369 mini = i;
3370 }
3371
3372 double tmp = work[k] / data (mini);
3373 work[k] = tmp;
3374 for (octave_idx_type i = cidx (k);
3375 i < cidx (k+1); i++)
3376 {
3377 if (i == mini)
3378 continue;
3379
3380 octave_idx_type iidx = inv_perm[ridx (i)];
3381 work[iidx] = work[iidx] - tmp * data (i);
3382 }
3383 }
3384 }
3385
3386 double atmp = 0;
3387 for (octave_idx_type i = j; i < nc; i++)
3388 {
3389 atmp += fabs (work[i]);
3390 work[i] = 0.;
3391 }
3392 if (atmp > ainvnorm)
3393 ainvnorm = atmp;
3394 }
3395 rcond = 1. / ainvnorm / anorm;
3396 }
3397 }
3398 else
3399 {
3400 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3401 retval.resize (nc, b_nc, 0.);
3402
3403 for (octave_idx_type j = 0; j < b_nc; j++)
3404 {
3405 for (octave_idx_type i = 0; i < nr; i++)
3406 cwork[i] = b(i, j);
3407 for (octave_idx_type i = nr; i < nc; i++)
3408 cwork[i] = 0.;
3409
3410 for (octave_idx_type k = 0; k < nc; k++)
3411 {
3412 if (cwork[k] != 0.)
3413 {
3414 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3415 {
3416 err = -2;
3417 goto triangular_error;
3418 }
3419
3420 Complex tmp = cwork[k] / data (cidx (k));
3421 cwork[k] = tmp;
3422 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3423 {
3424 octave_idx_type iidx = ridx (i);
3425 cwork[iidx] = cwork[iidx] - tmp * data (i);
3426 }
3427 }
3428 }
3429
3430 for (octave_idx_type i = 0; i < nc; i++)
3431 retval.xelem (i, j) = cwork[i];
3432 }
3433
3434 if (calc_cond)
3435 {
3436 // Calculation of 1-norm of inv(*this)
3437 OCTAVE_LOCAL_BUFFER (double, work, nm);
3438 for (octave_idx_type i = 0; i < nm; i++)
3439 work[i] = 0.;
3440
3441 for (octave_idx_type j = 0; j < nr; j++)
3442 {
3443 work[j] = 1.;
3444
3445 for (octave_idx_type k = j; k < nc; k++)
3446 {
3447
3448 if (work[k] != 0.)
3449 {
3450 double tmp = work[k] / data (cidx (k));
3451 work[k] = tmp;
3452 for (octave_idx_type i = cidx (k)+1;
3453 i < cidx (k+1); i++)
3454 {
3455 octave_idx_type iidx = ridx (i);
3456 work[iidx] = work[iidx] - tmp * data (i);
3457 }
3458 }
3459 }
3460 double atmp = 0;
3461 for (octave_idx_type i = j; i < nc; i++)
3462 {
3463 atmp += fabs (work[i]);
3464 work[i] = 0.;
3465 }
3466 if (atmp > ainvnorm)
3467 ainvnorm = atmp;
3468 }
3469 rcond = 1. / ainvnorm / anorm;
3470 }
3471 }
3472
3473 triangular_error:
3474 if (err != 0)
3475 {
3476 if (sing_handler)
3477 {
3478 sing_handler (rcond);
3479 mattype.mark_as_rectangular ();
3480 }
3481 else
3482 octave::warn_singular_matrix (rcond);
3483 }
3484
3485 if (is_singular (rcond) || octave::math::isnan (rcond))
3486 {
3487 err = -2;
3488
3489 if (sing_handler)
3490 {
3491 sing_handler (rcond);
3492 mattype.mark_as_rectangular ();
3493 }
3494 else
3495 octave::warn_singular_matrix (rcond);
3496 }
3497 }
3498
3499 return retval;
3500}
3501
3503SparseMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
3504 octave_idx_type& err, double& rcond,
3505 solve_singularity_handler sing_handler,
3506 bool calc_cond) const
3507{
3508 SparseComplexMatrix retval;
3509
3510 octave_idx_type nr = rows ();
3511 octave_idx_type nc = cols ();
3512 octave_idx_type nm = std::max (nc, nr);
3513 err = 0;
3514
3515 if (nr != b.rows ())
3516 (*current_liboctave_error_handler)
3517 ("matrix dimension mismatch solution of linear equations");
3518
3519 if (nr == 0 || nc == 0 || b.cols () == 0)
3520 retval = SparseComplexMatrix (nc, b.cols ());
3521 else
3522 {
3523 // Print spparms("spumoni") info if requested
3524 int typ = mattype.type ();
3525 mattype.info ();
3526
3527 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3528 (*current_liboctave_error_handler) ("incorrect matrix type");
3529
3530 double anorm = 0.;
3531 double ainvnorm = 0.;
3532 rcond = 1.;
3533
3534 if (calc_cond)
3535 {
3536 // Calculate the 1-norm of matrix for rcond calculation
3537 for (octave_idx_type j = 0; j < nc; j++)
3538 {
3539 double atmp = 0.;
3540 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3541 atmp += fabs (data (i));
3542 if (atmp > anorm)
3543 anorm = atmp;
3544 }
3545 }
3546
3547 octave_idx_type b_nc = b.cols ();
3548 octave_idx_type b_nz = b.nnz ();
3549 retval = SparseComplexMatrix (nc, b_nc, b_nz);
3550 retval.xcidx (0) = 0;
3551 octave_idx_type ii = 0;
3552 octave_idx_type x_nz = b_nz;
3553
3554 if (typ == MatrixType::Permuted_Lower)
3555 {
3556 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3557
3558 octave_idx_type *perm = mattype.triangular_perm ();
3559 OCTAVE_LOCAL_BUFFER (octave_idx_type, inv_perm, nr);
3560 for (octave_idx_type i = 0; i < nr; i++)
3561 inv_perm[perm[i]] = i;
3562
3563 for (octave_idx_type j = 0; j < b_nc; j++)
3564 {
3565 for (octave_idx_type i = 0; i < nm; i++)
3566 cwork[i] = 0.;
3567 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3568 cwork[inv_perm[b.ridx (i)]] = b.data (i);
3569
3570 for (octave_idx_type k = 0; k < nc; k++)
3571 {
3572 if (cwork[k] != 0.)
3573 {
3574 octave_idx_type minr = nr;
3575 octave_idx_type mini = 0;
3576
3577 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3578 if (inv_perm[ridx (i)] < minr)
3579 {
3580 minr = inv_perm[ridx (i)];
3581 mini = i;
3582 }
3583
3584 if (minr != k || data (mini) == 0.)
3585 {
3586 err = -2;
3587 goto triangular_error;
3588 }
3589
3590 Complex tmp = cwork[k] / data (mini);
3591 cwork[k] = tmp;
3592 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3593 {
3594 if (i == mini)
3595 continue;
3596
3597 octave_idx_type iidx = inv_perm[ridx (i)];
3598 cwork[iidx] = cwork[iidx] - tmp * data (i);
3599 }
3600 }
3601 }
3602
3603 // Count nonzeros in work vector and adjust space in
3604 // retval if needed
3605 octave_idx_type new_nnz = 0;
3606 for (octave_idx_type i = 0; i < nc; i++)
3607 if (cwork[i] != 0.)
3608 new_nnz++;
3609
3610 if (ii + new_nnz > x_nz)
3611 {
3612 // Resize the sparse matrix
3613 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3614 retval.change_capacity (sz);
3615 x_nz = sz;
3616 }
3617
3618 for (octave_idx_type i = 0; i < nc; i++)
3619 if (cwork[i] != 0.)
3620 {
3621 retval.xridx (ii) = i;
3622 retval.xdata (ii++) = cwork[i];
3623 }
3624 retval.xcidx (j+1) = ii;
3625 }
3626
3627 retval.maybe_compress ();
3628
3629 if (calc_cond)
3630 {
3631 // Calculation of 1-norm of inv(*this)
3632 OCTAVE_LOCAL_BUFFER (double, work, nm);
3633 for (octave_idx_type i = 0; i < nm; i++)
3634 work[i] = 0.;
3635
3636 for (octave_idx_type j = 0; j < nr; j++)
3637 {
3638 work[j] = 1.;
3639
3640 for (octave_idx_type k = 0; k < nc; k++)
3641 {
3642 if (work[k] != 0.)
3643 {
3644 octave_idx_type minr = nr;
3645 octave_idx_type mini = 0;
3646
3647 for (octave_idx_type i = cidx (k);
3648 i < cidx (k+1); i++)
3649 if (inv_perm[ridx (i)] < minr)
3650 {
3651 minr = inv_perm[ridx (i)];
3652 mini = i;
3653 }
3654
3655 double tmp = work[k] / data (mini);
3656 work[k] = tmp;
3657 for (octave_idx_type i = cidx (k);
3658 i < cidx (k+1); i++)
3659 {
3660 if (i == mini)
3661 continue;
3662
3663 octave_idx_type iidx = inv_perm[ridx (i)];
3664 work[iidx] = work[iidx] - tmp * data (i);
3665 }
3666 }
3667 }
3668
3669 double atmp = 0;
3670 for (octave_idx_type i = j; i < nc; i++)
3671 {
3672 atmp += fabs (work[i]);
3673 work[i] = 0.;
3674 }
3675 if (atmp > ainvnorm)
3676 ainvnorm = atmp;
3677 }
3678 rcond = 1. / ainvnorm / anorm;
3679 }
3680 }
3681 else
3682 {
3683 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3684
3685 for (octave_idx_type j = 0; j < b_nc; j++)
3686 {
3687 for (octave_idx_type i = 0; i < nm; i++)
3688 cwork[i] = 0.;
3689 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3690 cwork[b.ridx (i)] = b.data (i);
3691
3692 for (octave_idx_type k = 0; k < nc; k++)
3693 {
3694 if (cwork[k] != 0.)
3695 {
3696 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3697 {
3698 err = -2;
3699 goto triangular_error;
3700 }
3701
3702 Complex tmp = cwork[k] / data (cidx (k));
3703 cwork[k] = tmp;
3704 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3705 {
3706 octave_idx_type iidx = ridx (i);
3707 cwork[iidx] = cwork[iidx] - tmp * data (i);
3708 }
3709 }
3710 }
3711
3712 // Count nonzeros in work vector and adjust space in
3713 // retval if needed
3714 octave_idx_type new_nnz = 0;
3715 for (octave_idx_type i = 0; i < nc; i++)
3716 if (cwork[i] != 0.)
3717 new_nnz++;
3718
3719 if (ii + new_nnz > x_nz)
3720 {
3721 // Resize the sparse matrix
3722 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3723 retval.change_capacity (sz);
3724 x_nz = sz;
3725 }
3726
3727 for (octave_idx_type i = 0; i < nc; i++)
3728 if (cwork[i] != 0.)
3729 {
3730 retval.xridx (ii) = i;
3731 retval.xdata (ii++) = cwork[i];
3732 }
3733 retval.xcidx (j+1) = ii;
3734 }
3735
3736 retval.maybe_compress ();
3737
3738 if (calc_cond)
3739 {
3740 // Calculation of 1-norm of inv(*this)
3741 OCTAVE_LOCAL_BUFFER (double, work, nm);
3742 for (octave_idx_type i = 0; i < nm; i++)
3743 work[i] = 0.;
3744
3745 for (octave_idx_type j = 0; j < nr; j++)
3746 {
3747 work[j] = 1.;
3748
3749 for (octave_idx_type k = j; k < nc; k++)
3750 {
3751
3752 if (work[k] != 0.)
3753 {
3754 double tmp = work[k] / data (cidx (k));
3755 work[k] = tmp;
3756 for (octave_idx_type i = cidx (k)+1;
3757 i < cidx (k+1); i++)
3758 {
3759 octave_idx_type iidx = ridx (i);
3760 work[iidx] = work[iidx] - tmp * data (i);
3761 }
3762 }
3763 }
3764 double atmp = 0;
3765 for (octave_idx_type i = j; i < nc; i++)
3766 {
3767 atmp += fabs (work[i]);
3768 work[i] = 0.;
3769 }
3770 if (atmp > ainvnorm)
3771 ainvnorm = atmp;
3772 }
3773 rcond = 1. / ainvnorm / anorm;
3774 }
3775 }
3776
3777 triangular_error:
3778 if (err != 0)
3779 {
3780 if (sing_handler)
3781 {
3782 sing_handler (rcond);
3783 mattype.mark_as_rectangular ();
3784 }
3785 else
3786 octave::warn_singular_matrix (rcond);
3787 }
3788
3789 if (is_singular (rcond) || octave::math::isnan (rcond))
3790 {
3791 err = -2;
3792
3793 if (sing_handler)
3794 {
3795 sing_handler (rcond);
3796 mattype.mark_as_rectangular ();
3797 }
3798 else
3799 octave::warn_singular_matrix (rcond);
3800 }
3801 }
3802
3803 return retval;
3804}
3805
3806Matrix
3807SparseMatrix::trisolve (MatrixType& mattype, const Matrix& b,
3808 octave_idx_type& err, double& rcond,
3809 solve_singularity_handler sing_handler,
3810 bool calc_cond) const
3811{
3812 Matrix retval;
3813
3814 octave_idx_type nr = rows ();
3815 octave_idx_type nc = cols ();
3816 err = 0;
3817
3818 if (nr != nc || nr != b.rows ())
3819 (*current_liboctave_error_handler)
3820 ("matrix dimension mismatch solution of linear equations");
3821
3822 if (nr == 0 || b.cols () == 0)
3823 retval = Matrix (nc, b.cols (), 0.0);
3824 else if (calc_cond)
3825 (*current_liboctave_error_handler)
3826 ("calculation of condition number not implemented");
3827 else
3828 {
3829 // Print spparms("spumoni") info if requested
3830 int typ = mattype.type ();
3831 mattype.info ();
3832
3834 {
3835 OCTAVE_LOCAL_BUFFER (double, D, nr);
3836 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
3837
3838 if (mattype.is_dense ())
3839 {
3840 octave_idx_type ii = 0;
3841
3842 for (octave_idx_type j = 0; j < nc-1; j++)
3843 {
3844 D[j] = data (ii++);
3845 DL[j] = data (ii);
3846 ii += 2;
3847 }
3848 D[nc-1] = data (ii);
3849 }
3850 else
3851 {
3852 D[0] = 0.;
3853 for (octave_idx_type i = 0; i < nr - 1; i++)
3854 {
3855 D[i+1] = 0.;
3856 DL[i] = 0.;
3857 }
3858
3859 for (octave_idx_type j = 0; j < nc; j++)
3860 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3861 {
3862 if (ridx (i) == j)
3863 D[j] = data (i);
3864 else if (ridx (i) == j + 1)
3865 DL[j] = data (i);
3866 }
3867 }
3868
3869 F77_INT tmp_nr = octave::to_f77_int (nr);
3870
3871 F77_INT b_nr = octave::to_f77_int (b.rows ());
3872 F77_INT b_nc = octave::to_f77_int (b.cols ());
3873
3874 retval = b;
3875 double *result = retval.rwdata ();
3876
3877 F77_INT tmp_err = 0;
3878
3879 F77_XFCN (dptsv, DPTSV, (tmp_nr, b_nc, D, DL, result,
3880 b_nr, tmp_err));
3881
3882 err = tmp_err;
3883
3884 if (err != 0)
3885 {
3886 err = 0;
3887 mattype.mark_as_unsymmetric ();
3889 }
3890 else
3891 rcond = 1.;
3892 }
3893
3894 if (typ == MatrixType::Tridiagonal)
3895 {
3896 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
3897 OCTAVE_LOCAL_BUFFER (double, D, nr);
3898 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
3899
3900 if (mattype.is_dense ())
3901 {
3902 octave_idx_type ii = 0;
3903
3904 for (octave_idx_type j = 0; j < nc-1; j++)
3905 {
3906 D[j] = data (ii++);
3907 DL[j] = data (ii++);
3908 DU[j] = data (ii++);
3909 }
3910 D[nc-1] = data (ii);
3911 }
3912 else
3913 {
3914 D[0] = 0.;
3915 for (octave_idx_type i = 0; i < nr - 1; i++)
3916 {
3917 D[i+1] = 0.;
3918 DL[i] = 0.;
3919 DU[i] = 0.;
3920 }
3921
3922 for (octave_idx_type j = 0; j < nc; j++)
3923 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3924 {
3925 if (ridx (i) == j)
3926 D[j] = data (i);
3927 else if (ridx (i) == j + 1)
3928 DL[j] = data (i);
3929 else if (ridx (i) == j - 1)
3930 DU[j-1] = data (i);
3931 }
3932 }
3933
3934 F77_INT tmp_nr = octave::to_f77_int (nr);
3935
3936 F77_INT b_nr = octave::to_f77_int (b.rows ());
3937 F77_INT b_nc = octave::to_f77_int (b.cols ());
3938
3939 retval = b;
3940 double *result = retval.rwdata ();
3941
3942 F77_INT tmp_err = 0;
3943
3944 F77_XFCN (dgtsv, DGTSV, (tmp_nr, b_nc, DL, D, DU, result,
3945 b_nr, tmp_err));
3946
3947 err = tmp_err;
3948
3949 if (err != 0)
3950 {
3951 rcond = 0.;
3952 err = -2;
3953
3954 if (sing_handler)
3955 {
3956 sing_handler (rcond);
3957 mattype.mark_as_rectangular ();
3958 }
3959 else
3960 octave::warn_singular_matrix ();
3961 }
3962 else
3963 rcond = 1.;
3964 }
3965 else if (typ != MatrixType::Tridiagonal_Hermitian)
3966 (*current_liboctave_error_handler) ("incorrect matrix type");
3967 }
3968
3969 return retval;
3970}
3971
3973SparseMatrix::trisolve (MatrixType& mattype, const SparseMatrix& b,
3974 octave_idx_type& err, double& rcond,
3975 solve_singularity_handler sing_handler,
3976 bool calc_cond) const
3977{
3978 SparseMatrix retval;
3979
3980 octave_idx_type nr = rows ();
3981 octave_idx_type nc = cols ();
3982 err = 0;
3983
3984 if (nr != nc || nr != b.rows ())
3985 (*current_liboctave_error_handler)
3986 ("matrix dimension mismatch solution of linear equations");
3987
3988 if (nr == 0 || b.cols () == 0)
3989 retval = SparseMatrix (nc, b.cols ());
3990 else if (calc_cond)
3991 (*current_liboctave_error_handler)
3992 ("calculation of condition number not implemented");
3993 else
3994 {
3995 // Print spparms("spumoni") info if requested
3996 int typ = mattype.type ();
3997 mattype.info ();
3998
3999 // Note can't treat symmetric case as there is no dpttrf function
4000 if (typ == MatrixType::Tridiagonal
4002 {
4003 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
4004 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
4005 OCTAVE_LOCAL_BUFFER (double, D, nr);
4006 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
4007 Array<F77_INT> ipvt (dim_vector (nr, 1));
4008 F77_INT *pipvt = ipvt.rwdata ();
4009
4010 if (mattype.is_dense ())
4011 {
4012 octave_idx_type ii = 0;
4013
4014 for (octave_idx_type j = 0; j < nc-1; j++)
4015 {
4016 D[j] = data (ii++);
4017 DL[j] = data (ii++);
4018 DU[j] = data (ii++);
4019 }
4020 D[nc-1] = data (ii);
4021 }
4022 else
4023 {
4024 D[0] = 0.;
4025 for (octave_idx_type i = 0; i < nr - 1; i++)
4026 {
4027 D[i+1] = 0.;
4028 DL[i] = 0.;
4029 DU[i] = 0.;
4030 }
4031
4032 for (octave_idx_type j = 0; j < nc; j++)
4033 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4034 {
4035 if (ridx (i) == j)
4036 D[j] = data (i);
4037 else if (ridx (i) == j + 1)
4038 DL[j] = data (i);
4039 else if (ridx (i) == j - 1)
4040 DU[j-1] = data (i);
4041 }
4042 }
4043
4044 F77_INT tmp_nr = octave::to_f77_int (nr);
4045
4046 F77_INT tmp_err = 0;
4047
4048 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4049
4050 if (err != 0)
4051 {
4052 rcond = 0.0;
4053 err = -2;
4054
4055 if (sing_handler)
4056 {
4057 sing_handler (rcond);
4058 mattype.mark_as_rectangular ();
4059 }
4060 else
4061 octave::warn_singular_matrix ();
4062 }
4063 else
4064 {
4065 rcond = 1.0;
4066 char job = 'N';
4067 octave_idx_type x_nz = b.nnz ();
4068 octave_idx_type b_nc = b.cols ();
4069 retval = SparseMatrix (nr, b_nc, x_nz);
4070 retval.xcidx (0) = 0;
4071 octave_idx_type ii = 0;
4072
4073 OCTAVE_LOCAL_BUFFER (double, work, nr);
4074
4075 for (octave_idx_type j = 0; j < b_nc; j++)
4076 {
4077 for (octave_idx_type i = 0; i < nr; i++)
4078 work[i] = 0.;
4079 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
4080 work[b.ridx (i)] = b.data (i);
4081
4082 F77_INT b_nr = octave::to_f77_int (b.rows ());
4083
4084 F77_XFCN (dgttrs, DGTTRS,
4085 (F77_CONST_CHAR_ARG2 (&job, 1),
4086 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4087 work, b_nr, tmp_err
4088 F77_CHAR_ARG_LEN (1)));
4089
4090 err = tmp_err;
4091
4092 // Count nonzeros in work vector and adjust
4093 // space in retval if needed
4094 octave_idx_type new_nnz = 0;
4095 for (octave_idx_type i = 0; i < nr; i++)
4096 if (work[i] != 0.)
4097 new_nnz++;
4098
4099 if (ii + new_nnz > x_nz)
4100 {
4101 // Resize the sparse matrix
4102 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4103 retval.change_capacity (sz);
4104 x_nz = sz;
4105 }
4106
4107 for (octave_idx_type i = 0; i < nr; i++)
4108 if (work[i] != 0.)
4109 {
4110 retval.xridx (ii) = i;
4111 retval.xdata (ii++) = work[i];
4112 }
4113 retval.xcidx (j+1) = ii;
4114 }
4115
4116 retval.maybe_compress ();
4117 }
4118 }
4119 else if (typ != MatrixType::Tridiagonal_Hermitian)
4120 (*current_liboctave_error_handler) ("incorrect matrix type");
4121 }
4122
4123 return retval;
4124}
4125
4127SparseMatrix::trisolve (MatrixType& mattype, const ComplexMatrix& b,
4128 octave_idx_type& err, double& rcond,
4129 solve_singularity_handler sing_handler,
4130 bool calc_cond) const
4131{
4132 ComplexMatrix retval;
4133
4134 octave_idx_type nr = rows ();
4135 octave_idx_type nc = cols ();
4136 err = 0;
4137
4138 if (nr != nc || nr != b.rows ())
4139 (*current_liboctave_error_handler)
4140 ("matrix dimension mismatch solution of linear equations");
4141
4142 if (nr == 0 || b.cols () == 0)
4143 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4144 else if (calc_cond)
4145 (*current_liboctave_error_handler)
4146 ("calculation of condition number not implemented");
4147 else
4148 {
4149 // Print spparms("spumoni") info if requested
4150 int typ = mattype.type ();
4151 mattype.info ();
4152
4154 {
4155 OCTAVE_LOCAL_BUFFER (double, D, nr);
4156 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4157
4158 if (mattype.is_dense ())
4159 {
4160 octave_idx_type ii = 0;
4161
4162 for (octave_idx_type j = 0; j < nc-1; j++)
4163 {
4164 D[j] = data (ii++);
4165 DL[j] = data (ii);
4166 ii += 2;
4167 }
4168 D[nc-1] = data (ii);
4169 }
4170 else
4171 {
4172 D[0] = 0.;
4173 for (octave_idx_type i = 0; i < nr - 1; i++)
4174 {
4175 D[i+1] = 0.;
4176 DL[i] = 0.;
4177 }
4178
4179 for (octave_idx_type j = 0; j < nc; j++)
4180 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4181 {
4182 if (ridx (i) == j)
4183 D[j] = data (i);
4184 else if (ridx (i) == j + 1)
4185 DL[j] = data (i);
4186 }
4187 }
4188
4189 F77_INT tmp_nr = octave::to_f77_int (nr);
4190
4191 F77_INT b_nr = octave::to_f77_int (b.rows ());
4192 F77_INT b_nc = octave::to_f77_int (b.cols ());
4193
4194 rcond = 1.;
4195
4196 retval = b;
4197 Complex *result = retval.rwdata ();
4198
4199 F77_INT tmp_err = 0;
4200
4201 F77_XFCN (zptsv, ZPTSV, (tmp_nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL),
4202 F77_DBLE_CMPLX_ARG (result),
4203 b_nr, tmp_err));
4204
4205 err = tmp_err;
4206
4207 if (err != 0)
4208 {
4209 err = 0;
4210 mattype.mark_as_unsymmetric ();
4212 }
4213 }
4214
4215 if (typ == MatrixType::Tridiagonal)
4216 {
4217 OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4219 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4220
4221 if (mattype.is_dense ())
4222 {
4223 octave_idx_type ii = 0;
4224
4225 for (octave_idx_type j = 0; j < nc-1; j++)
4226 {
4227 D[j] = data (ii++);
4228 DL[j] = data (ii++);
4229 DU[j] = data (ii++);
4230 }
4231 D[nc-1] = data (ii);
4232 }
4233 else
4234 {
4235 D[0] = 0.;
4236 for (octave_idx_type i = 0; i < nr - 1; i++)
4237 {
4238 D[i+1] = 0.;
4239 DL[i] = 0.;
4240 DU[i] = 0.;
4241 }
4242
4243 for (octave_idx_type j = 0; j < nc; j++)
4244 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4245 {
4246 if (ridx (i) == j)
4247 D[j] = data (i);
4248 else if (ridx (i) == j + 1)
4249 DL[j] = data (i);
4250 else if (ridx (i) == j - 1)
4251 DU[j-1] = data (i);
4252 }
4253 }
4254
4255 F77_INT tmp_nr = octave::to_f77_int (nr);
4256
4257 F77_INT b_nr = octave::to_f77_int (b.rows ());
4258 F77_INT b_nc = octave::to_f77_int (b.cols ());
4259
4260 rcond = 1.;
4261
4262 retval = b;
4263 Complex *result = retval.rwdata ();
4264
4265 F77_INT tmp_err = 0;
4266
4267 F77_XFCN (zgtsv, ZGTSV, (tmp_nr, b_nc, F77_DBLE_CMPLX_ARG (DL),
4269 F77_DBLE_CMPLX_ARG (DU),
4270 F77_DBLE_CMPLX_ARG (result),
4271 b_nr, tmp_err));
4272
4273 err = tmp_err;
4274
4275 if (err != 0)
4276 {
4277 rcond = 0.;
4278 err = -2;
4279
4280 if (sing_handler)
4281 {
4282 sing_handler (rcond);
4283 mattype.mark_as_rectangular ();
4284 }
4285 else
4286 octave::warn_singular_matrix ();
4287 }
4288 }
4289 else if (typ != MatrixType::Tridiagonal_Hermitian)
4290 (*current_liboctave_error_handler) ("incorrect matrix type");
4291 }
4292
4293 return retval;
4294}
4295
4297SparseMatrix::trisolve (MatrixType& mattype, const SparseComplexMatrix& b,
4298 octave_idx_type& err, double& rcond,
4299 solve_singularity_handler sing_handler,
4300 bool calc_cond) const
4301{
4302 SparseComplexMatrix retval;
4303
4304 octave_idx_type nr = rows ();
4305 octave_idx_type nc = cols ();
4306 err = 0;
4307
4308 if (nr != nc || nr != b.rows ())
4309 (*current_liboctave_error_handler)
4310 ("matrix dimension mismatch solution of linear equations");
4311
4312 if (nr == 0 || b.cols () == 0)
4313 retval = SparseComplexMatrix (nc, b.cols ());
4314 else if (calc_cond)
4315 (*current_liboctave_error_handler)
4316 ("calculation of condition number not implemented");
4317 else
4318 {
4319 // Print spparms("spumoni") info if requested
4320 int typ = mattype.type ();
4321 mattype.info ();
4322
4323 // Note can't treat symmetric case as there is no dpttrf function
4324 if (typ == MatrixType::Tridiagonal
4326 {
4327 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
4328 OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
4329 OCTAVE_LOCAL_BUFFER (double, D, nr);
4330 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
4331 Array<F77_INT> ipvt (dim_vector (nr, 1));
4332 F77_INT *pipvt = ipvt.rwdata ();
4333
4334 if (mattype.is_dense ())
4335 {
4336 octave_idx_type ii = 0;
4337
4338 for (octave_idx_type j = 0; j < nc-1; j++)
4339 {
4340 D[j] = data (ii++);
4341 DL[j] = data (ii++);
4342 DU[j] = data (ii++);
4343 }
4344 D[nc-1] = data (ii);
4345 }
4346 else
4347 {
4348 D[0] = 0.;
4349 for (octave_idx_type i = 0; i < nr - 1; i++)
4350 {
4351 D[i+1] = 0.;
4352 DL[i] = 0.;
4353 DU[i] = 0.;
4354 }
4355
4356 for (octave_idx_type j = 0; j < nc; j++)
4357 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4358 {
4359 if (ridx (i) == j)
4360 D[j] = data (i);
4361 else if (ridx (i) == j + 1)
4362 DL[j] = data (i);
4363 else if (ridx (i) == j - 1)
4364 DU[j-1] = data (i);
4365 }
4366 }
4367
4368 F77_INT tmp_nr = octave::to_f77_int (nr);
4369
4370 F77_INT tmp_err = 0;
4371
4372 F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4373
4374 err = tmp_err;
4375
4376 if (err != 0)
4377 {
4378 rcond = 0.0;
4379 err = -2;
4380
4381 if (sing_handler)
4382 {
4383 sing_handler (rcond);
4384 mattype.mark_as_rectangular ();
4385 }
4386 else
4387 octave::warn_singular_matrix ();
4388 }
4389 else
4390 {
4391 rcond = 1.;
4392 char job = 'N';
4393 F77_INT b_nr = octave::to_f77_int (b.rows ());
4394 octave_idx_type b_nc = b.cols ();
4395 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4396 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
4397
4398 // Take a first guess that the number of nonzero terms
4399 // will be as many as in b
4400 octave_idx_type x_nz = b.nnz ();
4401 octave_idx_type ii = 0;
4402 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4403
4404 retval.xcidx (0) = 0;
4405 for (octave_idx_type j = 0; j < b_nc; j++)
4406 {
4407
4408 for (F77_INT i = 0; i < b_nr; i++)
4409 {
4410 Complex c = b(i, j);
4411 Bx[i] = c.real ();
4412 Bz[i] = c.imag ();
4413 }
4414
4415 F77_XFCN (dgttrs, DGTTRS,
4416 (F77_CONST_CHAR_ARG2 (&job, 1),
4417 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4418 Bx, b_nr, tmp_err
4419 F77_CHAR_ARG_LEN (1)));
4420
4421 err = tmp_err;
4422
4423 if (err != 0)
4424 {
4425 // FIXME: Should this be a warning?
4426 (*current_liboctave_error_handler)
4427 ("SparseMatrix::solve solve failed");
4428
4429 err = -1;
4430 break;
4431 }
4432
4433 F77_XFCN (dgttrs, DGTTRS,
4434 (F77_CONST_CHAR_ARG2 (&job, 1),
4435 tmp_nr, 1, DL, D, DU, DU2, pipvt,
4436 Bz, b_nr, tmp_err
4437 F77_CHAR_ARG_LEN (1)));
4438
4439 err = tmp_err;
4440
4441 if (err != 0)
4442 {
4443 // FIXME: Should this be a warning?
4444 (*current_liboctave_error_handler)
4445 ("SparseMatrix::solve solve failed");
4446
4447 err = -1;
4448 break;
4449 }
4450
4451 // Count nonzeros in work vector and adjust
4452 // space in retval if needed
4453 octave_idx_type new_nnz = 0;
4454 for (octave_idx_type i = 0; i < nr; i++)
4455 if (Bx[i] != 0. || Bz[i] != 0.)
4456 new_nnz++;
4457
4458 if (ii + new_nnz > x_nz)
4459 {
4460 // Resize the sparse matrix
4461 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4462 retval.change_capacity (sz);
4463 x_nz = sz;
4464 }
4465
4466 for (octave_idx_type i = 0; i < nr; i++)
4467 if (Bx[i] != 0. || Bz[i] != 0.)
4468 {
4469 retval.xridx (ii) = i;
4470 retval.xdata (ii++) = Complex (Bx[i], Bz[i]);
4471 }
4472
4473 retval.xcidx (j+1) = ii;
4474 }
4475
4476 retval.maybe_compress ();
4477 }
4478 }
4479 else if (typ != MatrixType::Tridiagonal_Hermitian)
4480 (*current_liboctave_error_handler) ("incorrect matrix type");
4481 }
4482
4483 return retval;
4484}
4485
4486Matrix
4487SparseMatrix::bsolve (MatrixType& mattype, const Matrix& b,
4488 octave_idx_type& err, double& rcond,
4489 solve_singularity_handler sing_handler,
4490 bool calc_cond) const
4491{
4492 Matrix retval;
4493
4494 octave_idx_type nr = rows ();
4495 octave_idx_type nc = cols ();
4496 err = 0;
4497
4498 if (nr != nc || nr != b.rows ())
4499 (*current_liboctave_error_handler)
4500 ("matrix dimension mismatch solution of linear equations");
4501
4502 if (nr == 0 || b.cols () == 0)
4503 retval = Matrix (nc, b.cols (), 0.0);
4504 else
4505 {
4506 // Print spparms("spumoni") info if requested
4507 int typ = mattype.type ();
4508 mattype.info ();
4509
4511 {
4512 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4513 F77_INT ldm = n_lower + 1;
4514 Matrix m_band (ldm, nc);
4515 double *tmp_data = m_band.rwdata ();
4516
4517 if (! mattype.is_dense ())
4518 {
4519 octave_idx_type ii = 0;
4520
4521 for (octave_idx_type j = 0; j < ldm; j++)
4522 for (octave_idx_type i = 0; i < nc; i++)
4523 tmp_data[ii++] = 0.;
4524 }
4525
4526 for (octave_idx_type j = 0; j < nc; j++)
4527 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4528 {
4529 octave_idx_type ri = ridx (i);
4530 if (ri >= j)
4531 m_band(ri - j, j) = data (i);
4532 }
4533
4534 // Calculate the norm of the matrix, for later use.
4535 double anorm;
4536 if (calc_cond)
4537 anorm = m_band.abs ().sum ().row (0).max ();
4538
4539 F77_INT tmp_nr = octave::to_f77_int (nr);
4540
4541 F77_INT tmp_err = 0;
4542
4543 char job = 'L';
4544 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4545 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4546 F77_CHAR_ARG_LEN (1)));
4547
4548 err = tmp_err;
4549
4550 if (err != 0)
4551 {
4552 // Matrix is not positive definite!! Fall through to
4553 // unsymmetric banded solver.
4554 mattype.mark_as_unsymmetric ();
4555 typ = MatrixType::Banded;
4556 rcond = 0.0;
4557 err = 0;
4558 }
4559 else
4560 {
4561 if (calc_cond)
4562 {
4563 Array<double> z (dim_vector (3 * nr, 1));
4564 double *pz = z.rwdata ();
4565 Array<F77_INT> iz (dim_vector (nr, 1));
4566 F77_INT *piz = iz.rwdata ();
4567
4568 F77_XFCN (dpbcon, DPBCON,
4569 (F77_CONST_CHAR_ARG2 (&job, 1),
4570 tmp_nr, n_lower, tmp_data, ldm,
4571 anorm, rcond, pz, piz, tmp_err
4572 F77_CHAR_ARG_LEN (1)));
4573
4574 err = tmp_err;
4575
4576 if (err != 0)
4577 err = -2;
4578
4579 if (is_singular (rcond) || octave::math::isnan (rcond))
4580 {
4581 err = -2;
4582
4583 if (sing_handler)
4584 {
4585 sing_handler (rcond);
4586 mattype.mark_as_rectangular ();
4587 }
4588 else
4589 octave::warn_singular_matrix (rcond);
4590 }
4591 }
4592 else
4593 rcond = 1.;
4594
4595 if (err == 0)
4596 {
4597 retval = b;
4598 double *result = retval.rwdata ();
4599
4600 F77_INT b_nr = octave::to_f77_int (b.rows ());
4601 F77_INT b_nc = octave::to_f77_int (b.cols ());
4602
4603 F77_XFCN (dpbtrs, DPBTRS,
4604 (F77_CONST_CHAR_ARG2 (&job, 1),
4605 tmp_nr, n_lower, b_nc, tmp_data,
4606 ldm, result, b_nr, tmp_err
4607 F77_CHAR_ARG_LEN (1)));
4608
4609 err = tmp_err;
4610
4611 if (err != 0)
4612 {
4613 // FIXME: Should this be a warning?
4614 (*current_liboctave_error_handler)
4615 ("SparseMatrix::solve solve failed");
4616 err = -1;
4617 }
4618 }
4619 }
4620 }
4621
4622 if (typ == MatrixType::Banded)
4623 {
4624 // Create the storage for the banded form of the sparse matrix
4625 F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
4626 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4627 F77_INT ldm = n_upper + 2 * n_lower + 1;
4628
4629 Matrix m_band (ldm, nc);
4630 double *tmp_data = m_band.rwdata ();
4631
4632 if (! mattype.is_dense ())
4633 {
4634 octave_idx_type ii = 0;
4635
4636 for (F77_INT j = 0; j < ldm; j++)
4637 for (octave_idx_type i = 0; i < nc; i++)
4638 tmp_data[ii++] = 0.;
4639 }
4640
4641 for (octave_idx_type j = 0; j < nc; j++)
4642 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4643 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4644
4645 // Calculate the norm of the matrix, for later use.
4646 double anorm = 0.0;
4647 if (calc_cond)
4648 {
4649 for (octave_idx_type j = 0; j < nr; j++)
4650 {
4651 double atmp = 0.;
4652 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4653 atmp += fabs (data (i));
4654 if (atmp > anorm)
4655 anorm = atmp;
4656 }
4657 }
4658
4659 F77_INT tmp_nr = octave::to_f77_int (nr);
4660
4661 Array<F77_INT> ipvt (dim_vector (nr, 1));
4662 F77_INT *pipvt = ipvt.rwdata ();
4663
4664 F77_INT tmp_err = 0;
4665
4666 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4667 tmp_data, ldm, pipvt, tmp_err));
4668
4669 err = tmp_err;
4670
4671 // Throw away extra info LAPACK gives so as to not
4672 // change output.
4673 if (err != 0)
4674 {
4675 err = -2;
4676 rcond = 0.0;
4677
4678 if (sing_handler)
4679 {
4680 sing_handler (rcond);
4681 mattype.mark_as_rectangular ();
4682 }
4683 else
4684 octave::warn_singular_matrix ();
4685 }
4686 else
4687 {
4688 if (calc_cond)
4689 {
4690 char job = '1';
4691 Array<double> z (dim_vector (3 * nr, 1));
4692 double *pz = z.rwdata ();
4693 Array<F77_INT> iz (dim_vector (nr, 1));
4694 F77_INT *piz = iz.rwdata ();
4695
4696 F77_INT tmp_nc = octave::to_f77_int (nc);
4697
4698 F77_XFCN (dgbcon, DGBCON,
4699 (F77_CONST_CHAR_ARG2 (&job, 1),
4700 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4701 anorm, rcond, pz, piz, tmp_err
4702 F77_CHAR_ARG_LEN (1)));
4703
4704 err = tmp_err;
4705
4706 if (err != 0)
4707 err = -2;
4708
4709 if (is_singular (rcond) || octave::math::isnan (rcond))
4710 {
4711 err = -2;
4712
4713 if (sing_handler)
4714 {
4715 sing_handler (rcond);
4716 mattype.mark_as_rectangular ();
4717 }
4718 else
4719 octave::warn_singular_matrix (rcond);
4720 }
4721 }
4722 else
4723 rcond = 1.;
4724
4725 if (err == 0)
4726 {
4727 retval = b;
4728 double *result = retval.rwdata ();
4729
4730 F77_INT b_nr = octave::to_f77_int (b.rows ());
4731 F77_INT b_nc = octave::to_f77_int (b.cols ());
4732
4733 char job = 'N';
4734 F77_XFCN (dgbtrs, DGBTRS,
4735 (F77_CONST_CHAR_ARG2 (&job, 1),
4736 tmp_nr, n_lower, n_upper, b_nc, tmp_data,
4737 ldm, pipvt, result, b_nr, tmp_err
4738 F77_CHAR_ARG_LEN (1)));
4739
4740 err = tmp_err;
4741 }
4742 }
4743 }
4744 else if (typ != MatrixType::Banded_Hermitian)
4745 (*current_liboctave_error_handler) ("incorrect matrix type");
4746 }
4747
4748 return retval;
4749}
4750
4752SparseMatrix::bsolve (MatrixType& mattype, const SparseMatrix& b,
4753 octave_idx_type& err, double& rcond,
4754 solve_singularity_handler sing_handler,
4755 bool calc_cond) const
4756{
4757 SparseMatrix retval;
4758
4759 octave_idx_type nr = rows ();
4760 octave_idx_type nc = cols ();
4761 err = 0;
4762
4763 if (nr != nc || nr != b.rows ())
4764 (*current_liboctave_error_handler)
4765 ("matrix dimension mismatch solution of linear equations");
4766
4767 if (nr == 0 || b.cols () == 0)
4768 retval = SparseMatrix (nc, b.cols ());
4769 else
4770 {
4771 // Print spparms("spumoni") info if requested
4772 int typ = mattype.type ();
4773 mattype.info ();
4774
4776 {
4777 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4778 F77_INT ldm = octave::to_f77_int (n_lower + 1);
4779
4780 Matrix m_band (ldm, nc);
4781 double *tmp_data = m_band.rwdata ();
4782
4783 if (! mattype.is_dense ())
4784 {
4785 octave_idx_type ii = 0;
4786
4787 for (F77_INT j = 0; j < ldm; j++)
4788 for (octave_idx_type i = 0; i < nc; i++)
4789 tmp_data[ii++] = 0.;
4790 }
4791
4792 for (octave_idx_type j = 0; j < nc; j++)
4793 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4794 {
4795 octave_idx_type ri = ridx (i);
4796 if (ri >= j)
4797 m_band(ri - j, j) = data (i);
4798 }
4799
4800 // Calculate the norm of the matrix, for later use.
4801 double anorm;
4802 if (calc_cond)
4803 anorm = m_band.abs ().sum ().row (0).max ();
4804
4805 F77_INT tmp_nr = octave::to_f77_int (nr);
4806
4807 F77_INT tmp_err = 0;
4808
4809 char job = 'L';
4810 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4811 tmp_nr, n_lower, tmp_data, ldm, tmp_err
4812 F77_CHAR_ARG_LEN (1)));
4813
4814 err = tmp_err;
4815
4816 if (err != 0)
4817 {
4818 mattype.mark_as_unsymmetric ();
4819 typ = MatrixType::Banded;
4820 rcond = 0.0;
4821 err = 0;
4822 }
4823 else
4824 {
4825 if (calc_cond)
4826 {
4827 Array<double> z (dim_vector (3 * nr, 1));
4828 double *pz = z.rwdata ();
4829 Array<F77_INT> iz (dim_vector (nr, 1));
4830 F77_INT *piz = iz.rwdata ();
4831
4832 F77_XFCN (dpbcon, DPBCON,
4833 (F77_CONST_CHAR_ARG2 (&job, 1),
4834 tmp_nr, n_lower, tmp_data, ldm,
4835 anorm, rcond, pz, piz, tmp_err
4836 F77_CHAR_ARG_LEN (1)));
4837
4838 err = tmp_err;
4839
4840 if (err != 0)
4841 err = -2;
4842
4843 if (is_singular (rcond) || octave::math::isnan (rcond))
4844 {
4845 err = -2;
4846
4847 if (sing_handler)
4848 {
4849 sing_handler (rcond);
4850 mattype.mark_as_rectangular ();
4851 }
4852 else
4853 octave::warn_singular_matrix (rcond);
4854 }
4855 }
4856 else
4857 rcond = 1.;
4858
4859 if (err == 0)
4860 {
4861 F77_INT b_nr = octave::to_f77_int (b.rows ());
4862 octave_idx_type b_nc = b.cols ();
4863 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4864
4865 // Take a first guess that the number of nonzero terms
4866 // will be as many as in b
4867 octave_idx_type x_nz = b.nnz ();
4868 octave_idx_type ii = 0;
4869 retval = SparseMatrix (b_nr, b_nc, x_nz);
4870
4871 retval.xcidx (0) = 0;
4872 for (octave_idx_type j = 0; j < b_nc; j++)
4873 {
4874 for (F77_INT i = 0; i < b_nr; i++)
4875 Bx[i] = b.elem (i, j);
4876
4877 F77_XFCN (dpbtrs, DPBTRS,
4878 (F77_CONST_CHAR_ARG2 (&job, 1),
4879 tmp_nr, n_lower, 1, tmp_data,
4880 ldm, Bx, b_nr, tmp_err
4881 F77_CHAR_ARG_LEN (1)));
4882
4883 err = tmp_err;
4884
4885 if (err != 0)
4886 {
4887 // FIXME: Should this be a warning?
4888 (*current_liboctave_error_handler)
4889 ("SparseMatrix::solve solve failed");
4890 err = -1;
4891 break;
4892 }
4893
4894 for (F77_INT i = 0; i < b_nr; i++)
4895 {
4896 double tmp = Bx[i];
4897 if (tmp != 0.0)
4898 {
4899 if (ii == x_nz)
4900 {
4901 // Resize the sparse matrix
4902 octave_idx_type sz;
4903 sz = (static_cast<double> (b_nc) - j) / b_nc
4904 * x_nz;
4905 sz = x_nz + (sz > 100 ? sz : 100);
4906 retval.change_capacity (sz);
4907 x_nz = sz;
4908 }
4909 retval.xdata (ii) = tmp;
4910 retval.xridx (ii++) = i;
4911 }
4912 }
4913 retval.xcidx (j+1) = ii;
4914 }
4915
4916 retval.maybe_compress ();
4917 }
4918 }
4919 }
4920
4921 if (typ == MatrixType::Banded)
4922 {
4923 // Create the storage for the banded form of the sparse matrix
4924 F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
4925 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4926 F77_INT ldm = octave::to_f77_int (n_upper + 2 * n_lower + 1);
4927
4928 Matrix m_band (ldm, nc);
4929 double *tmp_data = m_band.rwdata ();
4930
4931 if (! mattype.is_dense ())
4932 {
4933 octave_idx_type ii = 0;
4934
4935 for (octave_idx_type j = 0; j < ldm; j++)
4936 for (octave_idx_type i = 0; i < nc; i++)
4937 tmp_data[ii++] = 0.;
4938 }
4939
4940 for (octave_idx_type j = 0; j < nc; j++)
4941 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4942 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4943
4944 // Calculate the norm of the matrix, for later use.
4945 double anorm;
4946 if (calc_cond)
4947 {
4948 anorm = 0.0;
4949 for (octave_idx_type j = 0; j < nr; j++)
4950 {
4951 double atmp = 0.0;
4952 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4953 atmp += fabs (data (i));
4954 if (atmp > anorm)
4955 anorm = atmp;
4956 }
4957 }
4958
4959 F77_INT tmp_nr = octave::to_f77_int (nr);
4960
4961 Array<F77_INT> ipvt (dim_vector (nr, 1));
4962 F77_INT *pipvt = ipvt.rwdata ();
4963
4964 F77_INT tmp_err = 0;
4965
4966 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4967 tmp_data, ldm, pipvt, tmp_err));
4968
4969 err = tmp_err;
4970
4971 if (err != 0)
4972 {
4973 err = -2;
4974 rcond = 0.0;
4975
4976 if (sing_handler)
4977 {
4978 sing_handler (rcond);
4979 mattype.mark_as_rectangular ();
4980 }
4981 else
4982 octave::warn_singular_matrix ();
4983 }
4984 else
4985 {
4986 if (calc_cond)
4987 {
4988 char job = '1';
4989 Array<double> z (dim_vector (3 * nr, 1));
4990 double *pz = z.rwdata ();
4991 Array<F77_INT> iz (dim_vector (nr, 1));
4992 F77_INT *piz = iz.rwdata ();
4993
4994 F77_INT tmp_nc = octave::to_f77_int (nc);
4995
4996 F77_XFCN (dgbcon, DGBCON,
4997 (F77_CONST_CHAR_ARG2 (&job, 1),
4998 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4999 anorm, rcond, pz, piz, tmp_err
5000 F77_CHAR_ARG_LEN (1)));
5001
5002 err = tmp_err;
5003
5004 if (err != 0)
5005 err = -2;
5006
5007 if (is_singular (rcond) || octave::math::isnan (rcond))
5008 {
5009 err = -2;
5010
5011 if (sing_handler)
5012 {
5013 sing_handler (rcond);
5014 mattype.mark_as_rectangular ();
5015 }
5016 else
5017 octave::warn_singular_matrix (rcond);
5018 }
5019 }
5020 else
5021 rcond = 1.;
5022
5023 if (err == 0)
5024 {
5025 char job = 'N';
5026 octave_idx_type x_nz = b.nnz ();
5027 octave_idx_type b_nc = b.cols ();
5028 retval = SparseMatrix (nr, b_nc, x_nz);
5029 retval.xcidx (0) = 0;
5030 octave_idx_type ii = 0;
5031
5032 OCTAVE_LOCAL_BUFFER (double, work, nr);
5033
5034 for (octave_idx_type j = 0; j < b_nc; j++)
5035 {
5036 for (octave_idx_type i = 0; i < nr; i++)
5037 work[i] = 0.;
5038 for (octave_idx_type i = b.cidx (j);
5039 i < b.cidx (j+1); i++)
5040 work[b.ridx (i)] = b.data (i);
5041
5042 F77_INT b_nr = octave::to_f77_int (b.rows ());
5043
5044 F77_XFCN (dgbtrs, DGBTRS,
5045 (F77_CONST_CHAR_ARG2 (&job, 1),
5046 tmp_nr, n_lower, n_upper, 1, tmp_data,
5047 ldm, pipvt, work, b_nr, tmp_err
5048 F77_CHAR_ARG_LEN (1)));
5049
5050 err = tmp_err;
5051
5052 // Count nonzeros in work vector and adjust
5053 // space in retval if needed
5054 octave_idx_type new_nnz = 0;
5055 for (octave_idx_type i = 0; i < nr; i++)
5056 if (work[i] != 0.)
5057 new_nnz++;
5058
5059 if (ii + new_nnz > x_nz)
5060 {
5061 // Resize the sparse matrix
5062 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5063 retval.change_capacity (sz);
5064 x_nz = sz;
5065 }
5066
5067 for (octave_idx_type i = 0; i < nr; i++)
5068 if (work[i] != 0.)
5069 {
5070 retval.xridx (ii) = i;
5071 retval.xdata (ii++) = work[i];
5072 }
5073 retval.xcidx (j+1) = ii;
5074 }
5075
5076 retval.maybe_compress ();
5077 }
5078 }
5079 }
5080 else if (typ != MatrixType::Banded_Hermitian)
5081 (*current_liboctave_error_handler) ("incorrect matrix type");
5082 }
5083
5084 return retval;
5085}
5086
5088SparseMatrix::bsolve (MatrixType& mattype, const ComplexMatrix& b,
5089 octave_idx_type& err, double& rcond,
5090 solve_singularity_handler sing_handler,
5091 bool calc_cond) const
5092{
5093 ComplexMatrix retval;
5094
5095 octave_idx_type nr = rows ();
5096 octave_idx_type nc = cols ();
5097 err = 0;
5098
5099 if (nr != nc || nr != b.rows ())
5100 (*current_liboctave_error_handler)
5101 ("matrix dimension mismatch solution of linear equations");
5102
5103 if (nr == 0 || b.cols () == 0)
5104 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
5105 else
5106 {
5107 // Print spparms("spumoni") info if requested
5108 int typ = mattype.type ();
5109 mattype.info ();
5110
5112 {
5113 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
5114 F77_INT ldm = n_lower + 1;
5115
5116 Matrix m_band (ldm, nc);
5117 double *tmp_data = m_band.rwdata ();
5118
5119 if (! mattype.is_dense ())
5120 {
5121 octave_idx_type ii = 0;
5122
5123 for (F77_INT j = 0; j < ldm; j++)
5124 for (octave_idx_type i = 0; i < nc; i++)
5125 tmp_data[ii++] = 0.;
5126 }
5127
5128 for (octave_idx_type j = 0; j < nc; j++)
5129 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5130 {
5131 octave_idx_type ri = ridx (i);
5132 if (ri >= j)
5133 m_band(ri - j, j) = data (i);
5134 }
5135
5136 // Calculate the norm of the matrix, for later use.
5137 double anorm;
5138 if (calc_cond)
5139 anorm = m_band.abs ().sum ().row (0).max ();
5140
5141 F77_INT tmp_nr = octave::to_f77_int (nr);
5142
5143 F77_INT tmp_err = 0;
5144
5145 char job = 'L';
5146 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5147 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5148 F77_CHAR_ARG_LEN (1)));
5149
5150 err = tmp_err;
5151
5152 if (err != 0)
5153 {
5154 // Matrix is not positive definite!! Fall through to
5155 // unsymmetric banded solver.
5156 mattype.mark_as_unsymmetric ();
5157 typ = MatrixType::Banded;
5158 rcond = 0.0;
5159 err = 0;
5160 }
5161 else
5162 {
5163 if (calc_cond)
5164 {
5165 Array<double> z (dim_vector (3 * nr, 1));
5166 double *pz = z.rwdata ();
5167 Array<F77_INT> iz (dim_vector (nr, 1));
5168 F77_INT *piz = iz.rwdata ();
5169
5170 F77_XFCN (dpbcon, DPBCON,
5171 (F77_CONST_CHAR_ARG2 (&job, 1),
5172 tmp_nr, n_lower, tmp_data, ldm,
5173 anorm, rcond, pz, piz, tmp_err
5174 F77_CHAR_ARG_LEN (1)));
5175
5176 err = tmp_err;
5177
5178 if (err != 0)
5179 err = -2;
5180
5181 if (is_singular (rcond) || octave::math::isnan (rcond))
5182 {
5183 err = -2;
5184
5185 if (sing_handler)
5186 {
5187 sing_handler (rcond);
5188 mattype.mark_as_rectangular ();
5189 }
5190 else
5191 octave::warn_singular_matrix (rcond);
5192 }
5193 }
5194 else
5195 rcond = 1.;
5196
5197 if (err == 0)
5198 {
5199 F77_INT b_nr = octave::to_f77_int (b.rows ());
5200 octave_idx_type b_nc = b.cols ();
5201
5202 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
5203 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
5204
5205 retval.resize (b_nr, b_nc);
5206
5207 for (octave_idx_type j = 0; j < b_nc; j++)
5208 {
5209 for (F77_INT i = 0; i < b_nr; i++)
5210 {
5211 Complex c = b(i, j);
5212 Bx[i] = c.real ();
5213 Bz[i] = c.imag ();
5214 }
5215
5216 F77_XFCN (dpbtrs, DPBTRS,
5217 (F77_CONST_CHAR_ARG2 (&job, 1),
5218 tmp_nr, n_lower, 1, tmp_data,
5219 ldm, Bx, b_nr, tmp_err
5220 F77_CHAR_ARG_LEN (1)));
5221
5222 err = tmp_err;
5223
5224 if (err != 0)
5225 {
5226 // FIXME: Should this be a warning?
5227 (*current_liboctave_error_handler)
5228 ("SparseMatrix::solve solve failed");
5229 err = -1;
5230 break;
5231 }
5232
5233 F77_XFCN (dpbtrs, DPBTRS,
5234 (F77_CONST_CHAR_ARG2 (&job, 1),
5235 tmp_nr, n_lower, 1, tmp_data,
5236 ldm, Bz, b_nr, tmp_err
5237 F77_CHAR_ARG_LEN (1)));
5238
5239 err = tmp_err;
5240
5241 if (err != 0)
5242 {
5243 // FIXME: Should this be a warning?
5244 (*current_liboctave_error_handler)
5245 ("SparseMatrix::solve solve failed");
5246 err = -1;
5247 break;
5248 }
5249
5250 for (octave_idx_type i = 0; i < b_nr; i++)
5251 retval(i, j) = Complex (Bx[i], Bz[i]);
5252 }
5253 }
5254 }
5255 }
5256
5257 if (typ == MatrixType::Banded)
5258 {
5259 // Create the storage for the banded form of the sparse matrix
5260 F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
5261 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
5262 F77_INT ldm = n_upper + 2 * n_lower + 1;
5263
5264 Matrix m_band (ldm, nc);
5265 double *tmp_data = m_band.rwdata ();
5266
5267 if (! mattype.is_dense ())
5268 {
5269 octave_idx_type ii = 0;
5270
5271 for (F77_INT j = 0; j < ldm; j++)
5272 for (octave_idx_type i = 0; i < nc; i++)
5273 tmp_data[ii++] = 0.;
5274 }
5275
5276 for (octave_idx_type j = 0; j < nc; j++)
5277 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5278 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5279
5280 // Calculate the norm of the matrix, for later use.
5281 double anorm;
5282 if (calc_cond)
5283 {
5284 anorm = 0.0;
5285 for (octave_idx_type j = 0; j < nr; j++)
5286 {
5287 double atmp = 0.0;
5288 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5289 atmp += fabs (data (i));
5290 if (atmp > anorm)
5291 anorm = atmp;
5292 }
5293 }
5294
5295 F77_INT tmp_nr = octave::to_f77_int (nr);
5296
5297 Array<F77_INT> ipvt (dim_vector (nr, 1));
5298 F77_INT *pipvt = ipvt.rwdata ();
5299
5300 F77_INT tmp_err = 0;
5301
5302 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5303 tmp_data, ldm, pipvt, tmp_err));
5304
5305 err = tmp_err;
5306
5307 if (err != 0)
5308 {
5309 err = -2;
5310 rcond = 0.0;
5311
5312 if (sing_handler)
5313 {
5314 sing_handler (rcond);
5315 mattype.mark_as_rectangular ();
5316 }
5317 else
5318 octave::warn_singular_matrix ();
5319 }
5320 else
5321 {
5322 if (calc_cond)
5323 {
5324 char job = '1';
5325 Array<double> z (dim_vector (3 * nr, 1));
5326 double *pz = z.rwdata ();
5327 Array<F77_INT> iz (dim_vector (nr, 1));
5328 F77_INT *piz = iz.rwdata ();
5329
5330 F77_INT tmp_nc = octave::to_f77_int (nc);
5331
5332 F77_XFCN (dgbcon, DGBCON,
5333 (F77_CONST_CHAR_ARG2 (&job, 1),
5334 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5335 anorm, rcond, pz, piz, tmp_err
5336 F77_CHAR_ARG_LEN (1)));
5337
5338 err = tmp_err;
5339
5340 if (err != 0)
5341 err = -2;
5342
5343 if (is_singular (rcond) || octave::math::isnan (rcond))
5344 {
5345 err = -2;
5346
5347 if (sing_handler)
5348 {
5349 sing_handler (rcond);
5350 mattype.mark_as_rectangular ();
5351 }
5352 else
5353 octave::warn_singular_matrix (rcond);
5354 }
5355 }
5356 else
5357 rcond = 1.;
5358
5359 if (err == 0)
5360 {
5361 char job = 'N';
5362 F77_INT b_nr = octave::to_f77_int (b.rows ());
5363 octave_idx_type b_nc = b.cols ();
5364 retval.resize (nr, b_nc);
5365
5366 OCTAVE_LOCAL_BUFFER (double, Bz, nr);
5367 OCTAVE_LOCAL_BUFFER (double, Bx, nr);
5368
5369 for (octave_idx_type j = 0; j < b_nc; j++)
5370 {
5371 for (octave_idx_type i = 0; i < nr; i++)
5372 {
5373 Complex c = b(i, j);
5374 Bx[i] = c.real ();
5375 Bz[i] = c.imag ();
5376 }
5377
5378 F77_XFCN (dgbtrs, DGBTRS,
5379 (F77_CONST_CHAR_ARG2 (&job, 1),
5380 tmp_nr, n_lower, n_upper, 1, tmp_data,
5381 ldm, pipvt, Bx, b_nr, tmp_err
5382 F77_CHAR_ARG_LEN (1)));
5383
5384 err = tmp_err;
5385
5386 F77_XFCN (dgbtrs, DGBTRS,
5387 (F77_CONST_CHAR_ARG2 (&job, 1),
5388 tmp_nr, n_lower, n_upper, 1, tmp_data,
5389 ldm, pipvt, Bz, b_nr, tmp_err
5390 F77_CHAR_ARG_LEN (1)));
5391
5392 err = tmp_err;
5393
5394 for (octave_idx_type i = 0; i < nr; i++)
5395 retval(i, j) = Complex (Bx[i], Bz[i]);
5396 }
5397 }
5398 }
5399 }
5400 else if (typ != MatrixType::Banded_Hermitian)
5401 (*current_liboctave_error_handler) ("incorrect matrix type");
5402 }
5403
5404 return retval;
5405}
5406
5408SparseMatrix::bsolve (MatrixType& mattype, const SparseComplexMatrix& b,
5409 octave_idx_type& err, double& rcond,
5410 solve_singularity_handler sing_handler,
5411 bool calc_cond) const
5412{
5413 SparseComplexMatrix retval;
5414
5415 octave_idx_type nr = rows ();
5416 octave_idx_type nc = cols ();
5417 err = 0;
5418
5419 if (nr != nc || nr != b.rows ())
5420 (*current_liboctave_error_handler)
5421 ("matrix dimension mismatch solution of linear equations");
5422
5423 if (nr == 0 || b.cols () == 0)
5424 retval = SparseComplexMatrix (nc, b.cols ());
5425 else
5426 {
5427 // Print spparms("spumoni") info if requested
5428 int typ = mattype.type ();
5429 mattype.info ();
5430
5432 {
5433 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
5434 F77_INT ldm = n_lower + 1;
5435
5436 Matrix m_band (ldm, nc);
5437 double *tmp_data = m_band.rwdata ();
5438
5439 if (! mattype.is_dense ())
5440 {
5441 octave_idx_type ii = 0;
5442
5443 for (F77_INT j = 0; j < ldm; j++)
5444 for (octave_idx_type i = 0; i < nc; i++)
5445 tmp_data[ii++] = 0.;
5446 }
5447
5448 for (octave_idx_type j = 0; j < nc; j++)
5449 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5450 {
5451 octave_idx_type ri = ridx (i);
5452 if (ri >= j)
5453 m_band(ri - j, j) = data (i);
5454 }
5455
5456 // Calculate the norm of the matrix, for later use.
5457 double anorm;
5458 if (calc_cond)
5459 anorm = m_band.abs ().sum ().row (0).max ();
5460
5461 F77_INT tmp_nr = octave::to_f77_int (nr);
5462
5463 F77_INT tmp_err = 0;
5464
5465 char job = 'L';
5466 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5467 tmp_nr, n_lower, tmp_data, ldm, tmp_err
5468 F77_CHAR_ARG_LEN (1)));
5469
5470 err = tmp_err;
5471
5472 if (err != 0)
5473 {
5474 // Matrix is not positive definite!! Fall through to
5475 // unsymmetric banded solver.
5476 mattype.mark_as_unsymmetric ();
5477 typ = MatrixType::Banded;
5478
5479 rcond = 0.0;
5480 err = 0;
5481 }
5482 else
5483 {
5484 if (calc_cond)
5485 {
5486 Array<double> z (dim_vector (3 * nr, 1));
5487 double *pz = z.rwdata ();
5488 Array<F77_INT> iz (dim_vector (nr, 1));
5489 F77_INT *piz = iz.rwdata ();
5490
5491 F77_XFCN (dpbcon, DPBCON,
5492 (F77_CONST_CHAR_ARG2 (&job, 1),
5493 tmp_nr, n_lower, tmp_data, ldm,
5494 anorm, rcond, pz, piz, tmp_err
5495 F77_CHAR_ARG_LEN (1)));
5496
5497 err = tmp_err;
5498
5499 if (err != 0)
5500 err = -2;
5501
5502 if (is_singular (rcond) || octave::math::isnan (rcond))
5503 {
5504 err = -2;
5505
5506 if (sing_handler)
5507 {
5508 sing_handler (rcond);
5509 mattype.mark_as_rectangular ();
5510 }
5511 else
5512 octave::warn_singular_matrix (rcond);
5513 }
5514 }
5515 else
5516 rcond = 1.;
5517
5518 if (err == 0)
5519 {
5520 F77_INT b_nr = octave::to_f77_int (b.rows ());
5521 octave_idx_type b_nc = b.cols ();
5522 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
5523 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
5524
5525 // Take a first guess that the number of nonzero terms
5526 // will be as many as in b
5527 octave_idx_type x_nz = b.nnz ();
5528 octave_idx_type ii = 0;
5529 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5530
5531 retval.xcidx (0) = 0;
5532 for (octave_idx_type j = 0; j < b_nc; j++)
5533 {
5534
5535 for (F77_INT i = 0; i < b_nr; i++)
5536 {
5537 Complex c = b(i, j);
5538 Bx[i] = c.real ();
5539 Bz[i] = c.imag ();
5540 }
5541
5542 F77_XFCN (dpbtrs, DPBTRS,
5543 (F77_CONST_CHAR_ARG2 (&job, 1),
5544 tmp_nr, n_lower, 1, tmp_data,
5545 ldm, Bx, b_nr, tmp_err
5546 F77_CHAR_ARG_LEN (1)));
5547
5548 err = tmp_err;
5549
5550 if (err != 0)
5551 {
5552 // FIXME: Should this be a warning?
5553 (*current_liboctave_error_handler)
5554 ("SparseMatrix::solve solve failed");
5555 err = -1;
5556 break;
5557 }
5558
5559 F77_XFCN (dpbtrs, DPBTRS,
5560 (F77_CONST_CHAR_ARG2 (&job, 1),
5561 tmp_nr, n_lower, 1, tmp_data,
5562 ldm, Bz, b_nr, tmp_err
5563 F77_CHAR_ARG_LEN (1)));
5564
5565 err = tmp_err;
5566
5567 if (err != 0)
5568 {
5569 // FIXME: Should this be a warning?
5570 (*current_liboctave_error_handler)
5571 ("SparseMatrix::solve solve failed");
5572
5573 err = -1;
5574 break;
5575 }
5576
5577 // Count nonzeros in work vector and adjust
5578 // space in retval if needed
5579 octave_idx_type new_nnz = 0;
5580 for (octave_idx_type i = 0; i < nr; i++)
5581 if (Bx[i] != 0. || Bz[i] != 0.)
5582 new_nnz++;
5583
5584 if (ii + new_nnz > x_nz)
5585 {
5586 // Resize the sparse matrix
5587 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5588 retval.change_capacity (sz);
5589 x_nz = sz;
5590 }
5591
5592 for (octave_idx_type i = 0; i < nr; i++)
5593 if (Bx[i] != 0. || Bz[i] != 0.)
5594 {
5595 retval.xridx (ii) = i;
5596 retval.xdata (ii++) = Complex (Bx[i], Bz[i]);
5597 }
5598
5599 retval.xcidx (j+1) = ii;
5600 }
5601
5602 retval.maybe_compress ();
5603 }
5604 }
5605 }
5606
5607 if (typ == MatrixType::Banded)
5608 {
5609 // Create the storage for the banded form of the sparse matrix
5610 F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
5611 F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
5612 F77_INT ldm = n_upper + 2 * n_lower + 1;
5613
5614 Matrix m_band (ldm, nc);
5615 double *tmp_data = m_band.rwdata ();
5616
5617 if (! mattype.is_dense ())
5618 {
5619 octave_idx_type ii = 0;
5620
5621 for (F77_INT j = 0; j < ldm; j++)
5622 for (octave_idx_type i = 0; i < nc; i++)
5623 tmp_data[ii++] = 0.;
5624 }
5625
5626 for (octave_idx_type j = 0; j < nc; j++)
5627 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5628 m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5629
5630 // Calculate the norm of the matrix, for later use.
5631 double anorm;
5632 if (calc_cond)
5633 {
5634 anorm = 0.0;
5635 for (octave_idx_type j = 0; j < nr; j++)
5636 {
5637 double atmp = 0.0;
5638 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5639 atmp += fabs (data (i));
5640 if (atmp > anorm)
5641 anorm = atmp;
5642 }
5643 }
5644
5645 F77_INT tmp_nr = octave::to_f77_int (nr);
5646
5647 Array<F77_INT> ipvt (dim_vector (nr, 1));
5648 F77_INT *pipvt = ipvt.rwdata ();
5649
5650 F77_INT tmp_err = 0;
5651
5652 F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
5653 tmp_data, ldm, pipvt, tmp_err));
5654
5655 err = tmp_err;
5656
5657 if (err != 0)
5658 {
5659 err = -2;
5660 rcond = 0.0;
5661
5662 if (sing_handler)
5663 {
5664 sing_handler (rcond);
5665 mattype.mark_as_rectangular ();
5666 }
5667 else
5668 octave::warn_singular_matrix ();
5669 }
5670 else
5671 {
5672 if (calc_cond)
5673 {
5674 char job = '1';
5675 Array<double> z (dim_vector (3 * nr, 1));
5676 double *pz = z.rwdata ();
5677 Array<F77_INT> iz (dim_vector (nr, 1));
5678 F77_INT *piz = iz.rwdata ();
5679
5680 F77_INT tmp_nc = octave::to_f77_int (nc);
5681
5682 F77_XFCN (dgbcon, DGBCON,
5683 (F77_CONST_CHAR_ARG2 (&job, 1),
5684 tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5685 anorm, rcond, pz, piz, tmp_err
5686 F77_CHAR_ARG_LEN (1)));
5687
5688 err = tmp_err;
5689
5690 if (err != 0)
5691 err = -2;
5692
5693 if (is_singular (rcond) || octave::math::isnan (rcond))
5694 {
5695 err = -2;
5696
5697 if (sing_handler)
5698 {
5699 sing_handler (rcond);
5700 mattype.mark_as_rectangular ();
5701 }
5702 else
5703 octave::warn_singular_matrix (rcond);
5704 }
5705 }
5706 else
5707 rcond = 1.;
5708
5709 if (err == 0)
5710 {
5711 char job = 'N';
5712 octave_idx_type x_nz = b.nnz ();
5713 F77_INT b_nr = octave::to_f77_int (b.rows ());
5714 octave_idx_type b_nc = b.cols ();
5715 retval = SparseComplexMatrix (nr, b_nc, x_nz);
5716 retval.xcidx (0) = 0;
5717 octave_idx_type ii = 0;
5718
5719 OCTAVE_LOCAL_BUFFER (double, Bx, nr);
5720 OCTAVE_LOCAL_BUFFER (double, Bz, nr);
5721
5722 for (octave_idx_type j = 0; j < b_nc; j++)
5723 {
5724 for (octave_idx_type i = 0; i < nr; i++)
5725 {
5726 Bx[i] = 0.;
5727 Bz[i] = 0.;
5728 }
5729 for (octave_idx_type i = b.cidx (j);
5730 i < b.cidx (j+1); i++)
5731 {
5732 Complex c = b.data (i);
5733 Bx[b.ridx (i)] = c.real ();
5734 Bz[b.ridx (i)] = c.imag ();
5735 }
5736
5737 F77_XFCN (dgbtrs, DGBTRS,
5738 (F77_CONST_CHAR_ARG2 (&job, 1),
5739 tmp_nr, n_lower, n_upper, 1, tmp_data,
5740 ldm, pipvt, Bx, b_nr, tmp_err
5741 F77_CHAR_ARG_LEN (1)));
5742
5743 err = tmp_err;
5744
5745 F77_XFCN (dgbtrs, DGBTRS,
5746 (F77_CONST_CHAR_ARG2 (&job, 1),
5747 tmp_nr, n_lower, n_upper, 1, tmp_data,
5748 ldm, pipvt, Bz, b_nr, tmp_err
5749 F77_CHAR_ARG_LEN (1)));
5750
5751 err = tmp_err;
5752
5753 // Count nonzeros in work vector and adjust
5754 // space in retval if needed
5755 octave_idx_type new_nnz = 0;
5756 for (octave_idx_type i = 0; i < nr; i++)
5757 if (Bx[i] != 0. || Bz[i] != 0.)
5758 new_nnz++;
5759
5760 if (ii + new_nnz > x_nz)
5761 {
5762 // Resize the sparse matrix
5763 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5764 retval.change_capacity (sz);
5765 x_nz = sz;
5766 }
5767
5768 for (octave_idx_type i = 0; i < nr; i++)
5769 if (Bx[i] != 0. || Bz[i] != 0.)
5770 {
5771 retval.xridx (ii) = i;
5772 retval.xdata (ii++) = Complex (Bx[i], Bz[i]);
5773 }
5774 retval.xcidx (j+1) = ii;
5775 }
5776
5777 retval.maybe_compress ();
5778 }
5779 }
5780 }
5781 else if (typ != MatrixType::Banded_Hermitian)
5782 (*current_liboctave_error_handler) ("incorrect matrix type");
5783 }
5784
5785 return retval;
5786}
5787
5788void *
5789SparseMatrix::factorize (octave_idx_type& err, double& rcond, Matrix& Control,
5790 Matrix& Info, solve_singularity_handler sing_handler,
5791 bool calc_cond) const
5792{
5793 // The return values
5794 void *Numeric = nullptr;
5795 err = 0;
5796
5797#if defined (HAVE_UMFPACK)
5798
5799 // Setup the control parameters
5800 Control = Matrix (UMFPACK_CONTROL, 1);
5801 double *control = Control.rwdata ();
5802 UMFPACK_DNAME (defaults) (control);
5803
5804 double tmp = octave::sparse_params::get_key ("spumoni");
5805 if (! octave::math::isnan (tmp))
5806 Control (UMFPACK_PRL) = tmp;
5807 tmp = octave::sparse_params::get_key ("piv_tol");
5808 if (! octave::math::isnan (tmp))
5809 {
5810 Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5811 Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5812 }
5813
5814 // Set whether we are allowed to modify Q or not
5815 tmp = octave::sparse_params::get_key ("autoamd");
5816 if (! octave::math::isnan (tmp))
5817 Control (UMFPACK_FIXQ) = tmp;
5818
5819 UMFPACK_DNAME (report_control) (control);
5820
5821 const octave_idx_type *Ap = cidx ();
5822 const octave_idx_type *Ai = ridx ();
5823 const double *Ax = data ();
5824 octave_idx_type nr = rows ();
5825 octave_idx_type nc = cols ();
5826
5827 UMFPACK_DNAME (report_matrix) (nr, nc,
5828 octave::to_suitesparse_intptr (Ap),
5829 octave::to_suitesparse_intptr (Ai),
5830 Ax, 1, control);
5831
5832 void *Symbolic;
5833 Info = Matrix (1, UMFPACK_INFO);
5834 double *info = Info.rwdata ();
5835 int status = UMFPACK_DNAME (qsymbolic) (nr, nc,
5836 octave::to_suitesparse_intptr (Ap),
5837 octave::to_suitesparse_intptr (Ai),
5838 Ax, nullptr, &Symbolic, control, info);
5839
5840 if (status < 0)
5841 {
5842 UMFPACK_DNAME (report_status) (control, status);
5843 UMFPACK_DNAME (report_info) (control, info);
5844
5845 UMFPACK_DNAME (free_symbolic) (&Symbolic);
5846
5847 // FIXME: Should this be a warning?
5848 (*current_liboctave_error_handler)
5849 ("SparseMatrix::solve symbolic factorization failed");
5850 err = -1;
5851 }
5852 else
5853 {
5854 UMFPACK_DNAME (report_symbolic) (Symbolic, control);
5855
5856 status = UMFPACK_DNAME (numeric) (octave::to_suitesparse_intptr (Ap),
5857 octave::to_suitesparse_intptr (Ai),
5858 Ax, Symbolic, &Numeric, control, info);
5859 UMFPACK_DNAME (free_symbolic) (&Symbolic);
5860
5861 if (calc_cond)
5862 rcond = Info (UMFPACK_RCOND);
5863 else
5864 rcond = 1.;
5865
5866 if (status == UMFPACK_WARNING_singular_matrix
5867 || is_singular (rcond) || octave::math::isnan (rcond))
5868 {
5869 UMFPACK_DNAME (report_numeric) (Numeric, control);
5870
5871 err = -2;
5872
5873 if (sing_handler)
5874 sing_handler (rcond);
5875 else
5876 octave::warn_singular_matrix (rcond);
5877 }
5878 else if (status < 0)
5879 {
5880 UMFPACK_DNAME (report_status) (control, status);
5881 UMFPACK_DNAME (report_info) (control, info);
5882
5883 // FIXME: Should this be a warning?
5884 (*current_liboctave_error_handler)
5885 ("SparseMatrix::solve numeric factorization failed");
5886
5887 err = -1;
5888 }
5889 else
5890 {
5891 UMFPACK_DNAME (report_numeric) (Numeric, control);
5892 }
5893 }
5894
5895 if (err != 0)
5896 UMFPACK_DNAME (free_numeric) (&Numeric);
5897
5898#else
5899
5900 octave_unused_parameter (rcond);
5901 octave_unused_parameter (Control);
5902 octave_unused_parameter (Info);
5903 octave_unused_parameter (sing_handler);
5904 octave_unused_parameter (calc_cond);
5905
5906 (*current_liboctave_error_handler)
5907 ("support for UMFPACK was unavailable or disabled "
5908 "when liboctave was built");
5909
5910#endif
5911
5912 return Numeric;
5913}
5914
5915Matrix
5916SparseMatrix::fsolve (MatrixType& mattype, const Matrix& b,
5917 octave_idx_type& err, double& rcond,
5918 solve_singularity_handler sing_handler,
5919 bool calc_cond) const
5920{
5921 Matrix retval;
5922
5923 octave_idx_type nr = rows ();
5924 octave_idx_type nc = cols ();
5925 err = 0;
5926
5927 if (nr != nc || nr != b.rows ())
5928 (*current_liboctave_error_handler)
5929 ("matrix dimension mismatch solution of linear equations");
5930
5931 if (nr == 0 || b.cols () == 0)
5932 retval = Matrix (nc, b.cols (), 0.0);
5933 else
5934 {
5935 // Print spparms("spumoni") info if requested
5936 int typ = mattype.type ();
5937 mattype.info ();
5938
5939 if (typ == MatrixType::Hermitian)
5940 {
5941#if defined (HAVE_CHOLMOD)
5942 cholmod_common Common;
5943 cholmod_common *cm = &Common;
5944
5945 // Setup initial parameters
5946 CHOLMOD_NAME(start) (cm);
5947 cm->prefer_zomplex = false;
5948
5949 double spu = octave::sparse_params::get_key ("spumoni");
5950 if (spu == 0.)
5951 {
5952 cm->print = -1;
5953 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
5954 }
5955 else
5956 {
5957 cm->print = static_cast<int> (spu) + 2;
5958 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
5959 }
5960
5961 cm->error_handler = &SparseCholError;
5962 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5963 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5964
5965 cm->final_ll = true;
5966
5967 cholmod_sparse Astore;
5968 cholmod_sparse *A = &Astore;
5969 A->nrow = nr;
5970 A->ncol = nc;
5971
5972 A->p = cidx ();
5973 A->i = ridx ();
5974 A->nzmax = nnz ();
5975 A->packed = true;
5976 A->sorted = true;
5977 A->nz = nullptr;
5978#if defined (OCTAVE_ENABLE_64)
5979 A->itype = CHOLMOD_LONG;
5980#else
5981 A->itype = CHOLMOD_INT;
5982#endif
5983 A->dtype = CHOLMOD_DOUBLE;
5984 A->stype = 1;
5985 A->xtype = CHOLMOD_REAL;
5986
5987 A->x = data ();
5988
5989 cholmod_dense Bstore;
5990 cholmod_dense *B = &Bstore;
5991 B->nrow = b.rows ();
5992 B->ncol = b.cols ();
5993 B->d = B->nrow;
5994 B->nzmax = B->nrow * B->ncol;
5995 B->dtype = CHOLMOD_DOUBLE;
5996 B->xtype = CHOLMOD_REAL;
5997
5998 B->x = const_cast<double *> (b.data ());
5999
6000 cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6001 CHOLMOD_NAME(factorize) (A, L, cm);
6002 if (calc_cond)
6003 rcond = CHOLMOD_NAME(rcond)(L, cm);
6004 else
6005 rcond = 1.0;
6006
6007 if (rcond == 0.0)
6008 {
6009 // Either its indefinite or singular. Try UMFPACK
6010 mattype.mark_as_unsymmetric ();
6011 typ = MatrixType::Full;
6012 }
6013 else
6014 {
6015 if (is_singular (rcond) || octave::math::isnan (rcond))
6016 {
6017 err = -2;
6018
6019 if (sing_handler)
6020 {
6021 sing_handler (rcond);
6022 mattype.mark_as_rectangular ();
6023 }
6024 else
6025 octave::warn_singular_matrix (rcond);
6026
6027 return retval;
6028 }
6029
6030 cholmod_dense *X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
6031
6032 retval.resize (b.rows (), b.cols ());
6033 for (octave_idx_type j = 0; j < b.cols (); j++)
6034 {
6035 octave_idx_type jr = j * b.rows ();
6036 for (octave_idx_type i = 0; i < b.rows (); i++)
6037 retval.xelem (i, j) = static_cast<double *> (X->x)[jr + i];
6038 }
6039
6040 CHOLMOD_NAME(free_dense) (&X, cm);
6041 CHOLMOD_NAME(free_factor) (&L, cm);
6042 CHOLMOD_NAME(finish) (cm);
6043 static char blank_name[] = " ";
6044 CHOLMOD_NAME(print_common) (blank_name, cm);
6045 }
6046#else
6047 (*current_liboctave_warning_with_id_handler)
6048 ("Octave:missing-dependency",
6049 "support for CHOLMOD was unavailable or disabled "
6050 "when liboctave was built");
6051
6052 mattype.mark_as_unsymmetric ();
6053 typ = MatrixType::Full;
6054#endif
6055 }
6056
6057 if (typ == MatrixType::Full)
6058 {
6059#if defined (HAVE_UMFPACK)
6060 Matrix Control, Info;
6061 void *Numeric
6062 = factorize (err, rcond, Control, Info, sing_handler, calc_cond);
6063
6064 if (err == 0)
6065 {
6066 // one iterative refinement instead of the default two in UMFPACK
6067 Control (UMFPACK_IRSTEP) = 1;
6068 const double *Bx = b.data ();
6069 retval.resize (b.rows (), b.cols ());
6070 double *result = retval.rwdata ();
6071 octave_idx_type b_nr = b.rows ();
6072 octave_idx_type b_nc = b.cols ();
6073 int status = 0;
6074 double *control = Control.rwdata ();
6075 double *info = Info.rwdata ();
6076 const octave_idx_type *Ap = cidx ();
6077 const octave_idx_type *Ai = ridx ();
6078 const double *Ax = data ();
6079
6080 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
6081 {
6082 status = UMFPACK_DNAME (solve) (UMFPACK_A,
6083 octave::to_suitesparse_intptr (Ap),
6084 octave::to_suitesparse_intptr (Ai),
6085 Ax, &result[iidx], &Bx[iidx],
6086 Numeric, control, info);
6087 if (status < 0)
6088 {
6089 UMFPACK_DNAME (report_status) (control, status);
6090
6091 // FIXME: Should this be a warning?
6092 (*current_liboctave_error_handler)
6093 ("SparseMatrix::solve solve failed");
6094
6095 err = -1;
6096 break;
6097 }
6098 }
6099
6100 UMFPACK_DNAME (report_info) (control, info);
6101
6102 UMFPACK_DNAME (free_numeric) (&Numeric);
6103 }
6104 else
6105 mattype.mark_as_rectangular ();
6106
6107#else
6108 octave_unused_parameter (rcond);
6109 octave_unused_parameter (sing_handler);
6110 octave_unused_parameter (calc_cond);
6111
6112 (*current_liboctave_error_handler)
6113 ("support for UMFPACK was unavailable or disabled "
6114 "when liboctave was built");
6115#endif
6116 }
6117 else if (typ != MatrixType::Hermitian)
6118 (*current_liboctave_error_handler) ("incorrect matrix type");
6119 }
6120
6121 return retval;
6122}
6123
6125SparseMatrix::fsolve (MatrixType& mattype, const SparseMatrix& b,
6126 octave_idx_type& err, double& rcond,
6127 solve_singularity_handler sing_handler,
6128 bool calc_cond) const
6129{
6130 SparseMatrix retval;
6131
6132 octave_idx_type nr = rows ();
6133 octave_idx_type nc = cols ();
6134 err = 0;
6135
6136 if (nr != nc || nr != b.rows ())
6137 (*current_liboctave_error_handler)
6138 ("matrix dimension mismatch solution of linear equations");
6139
6140 if (nr == 0 || b.cols () == 0)
6141 retval = SparseMatrix (nc, b.cols ());
6142 else
6143 {
6144 // Print spparms("spumoni") info if requested
6145 int typ = mattype.type ();
6146 mattype.info ();
6147
6148 if (typ == MatrixType::Hermitian)
6149 {
6150#if defined (HAVE_CHOLMOD)
6151 cholmod_common Common;
6152 cholmod_common *cm = &Common;
6153
6154 // Setup initial parameters
6155 CHOLMOD_NAME(start) (cm);
6156 cm->prefer_zomplex = false;
6157
6158 double spu = octave::sparse_params::get_key ("spumoni");
6159 if (spu == 0.)
6160 {
6161 cm->print = -1;
6162 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6163 }
6164 else
6165 {
6166 cm->print = static_cast<int> (spu) + 2;
6167 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6168 }
6169
6170 cm->error_handler = &SparseCholError;
6171 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6172 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6173
6174 cm->final_ll = true;
6175
6176 cholmod_sparse Astore;
6177 cholmod_sparse *A = &Astore;
6178 A->nrow = nr;
6179 A->ncol = nc;
6180
6181 A->p = cidx ();
6182 A->i = ridx ();
6183 A->nzmax = nnz ();
6184 A->packed = true;
6185 A->sorted = true;
6186 A->nz = nullptr;
6187#if defined (OCTAVE_ENABLE_64)
6188 A->itype = CHOLMOD_LONG;
6189#else
6190 A->itype = CHOLMOD_INT;
6191#endif
6192 A->dtype = CHOLMOD_DOUBLE;
6193 A->stype = 1;
6194 A->xtype = CHOLMOD_REAL;
6195
6196 A->x = data ();
6197
6198 cholmod_sparse Bstore;
6199 cholmod_sparse *B = &Bstore;
6200 B->nrow = b.rows ();
6201 B->ncol = b.cols ();
6202 B->p = b.cidx ();
6203 B->i = b.ridx ();
6204 B->nzmax = b.nnz ();
6205 B->packed = true;
6206 B->sorted = true;
6207 B->nz = nullptr;
6208#if defined (OCTAVE_ENABLE_64)
6209 B->itype = CHOLMOD_LONG;
6210#else
6211 B->itype = CHOLMOD_INT;
6212#endif
6213 B->dtype = CHOLMOD_DOUBLE;
6214 B->stype = 0;
6215 B->xtype = CHOLMOD_REAL;
6216
6217 B->x = b.data ();
6218
6219 cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6220 CHOLMOD_NAME(factorize) (A, L, cm);
6221 if (calc_cond)
6222 rcond = CHOLMOD_NAME(rcond)(L, cm);
6223 else
6224 rcond = 1.;
6225
6226 if (rcond == 0.0)
6227 {
6228 // Either its indefinite or singular. Try UMFPACK
6229 mattype.mark_as_unsymmetric ();
6230 typ = MatrixType::Full;
6231 }
6232 else
6233 {
6234 if (is_singular (rcond) || octave::math::isnan (rcond))
6235 {
6236 err = -2;
6237
6238 if (sing_handler)
6239 {
6240 sing_handler (rcond);
6241 mattype.mark_as_rectangular ();
6242 }
6243 else
6244 octave::warn_singular_matrix (rcond);
6245
6246 return retval;
6247 }
6248
6249 cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6250
6251 retval = SparseMatrix
6252 (octave::from_size_t (X->nrow),
6253 octave::from_size_t (X->ncol),
6254 octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
6255 (X, cm)));
6256 for (octave_idx_type j = 0;
6257 j <= static_cast<octave_idx_type> (X->ncol); j++)
6258 retval.xcidx (j) = static_cast<octave_idx_type *> (X->p)[j];
6259 for (octave_idx_type j = 0;
6260 j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
6261 j++)
6262 {
6263 retval.xridx (j) = static_cast<octave_idx_type *> (X->i)[j];
6264 retval.xdata (j) = static_cast<double *> (X->x)[j];
6265 }
6266
6267 CHOLMOD_NAME(free_sparse) (&X, cm);
6268 CHOLMOD_NAME(free_factor) (&L, cm);
6269 CHOLMOD_NAME(finish) (cm);
6270 static char blank_name[] = " ";
6271 CHOLMOD_NAME(print_common) (blank_name, cm);
6272 }
6273#else
6274 (*current_liboctave_warning_with_id_handler)
6275 ("Octave:missing-dependency",
6276 "support for CHOLMOD was unavailable or disabled "
6277 "when liboctave was built");
6278
6279 mattype.mark_as_unsymmetric ();
6280 typ = MatrixType::Full;
6281#endif
6282 }
6283
6284 if (typ == MatrixType::Full)
6285 {
6286#if defined (HAVE_UMFPACK)
6287 Matrix Control, Info;
6288 void *Numeric = factorize (err, rcond, Control, Info,
6289 sing_handler, calc_cond);
6290
6291 if (err == 0)
6292 {
6293 // one iterative refinement instead of the default two in UMFPACK
6294 Control (UMFPACK_IRSTEP) = 1;
6295 octave_idx_type b_nr = b.rows ();
6296 octave_idx_type b_nc = b.cols ();
6297 int status = 0;
6298 double *control = Control.rwdata ();
6299 double *info = Info.rwdata ();
6300 const octave_idx_type *Ap = cidx ();
6301 const octave_idx_type *Ai = ridx ();
6302 const double *Ax = data ();
6303
6304 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6305 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6306
6307 // Take a first guess that the number of nonzero terms
6308 // will be as many as in b
6309 octave_idx_type x_nz = b.nnz ();
6310 octave_idx_type ii = 0;
6311 retval = SparseMatrix (b_nr, b_nc, x_nz);
6312
6313 retval.xcidx (0) = 0;
6314 for (octave_idx_type j = 0; j < b_nc; j++)
6315 {
6316
6317 for (octave_idx_type i = 0; i < b_nr; i++)
6318 Bx[i] = b.elem (i, j);
6319
6320 status = UMFPACK_DNAME (solve) (UMFPACK_A,
6321 octave::to_suitesparse_intptr (Ap),
6322 octave::to_suitesparse_intptr (Ai),
6323 Ax, Xx, Bx, Numeric,
6324 control, info);
6325 if (status < 0)
6326 {
6327 UMFPACK_DNAME (report_status) (control, status);
6328
6329 // FIXME: Should this be a warning?
6330 (*current_liboctave_error_handler)
6331 ("SparseMatrix::solve solve failed");
6332
6333 err = -1;
6334 break;
6335 }
6336
6337 for (octave_idx_type i = 0; i < b_nr; i++)
6338 {
6339 double tmp = Xx[i];
6340 if (tmp != 0.0)
6341 {
6342 if (ii == x_nz)
6343 {
6344 // Resize the sparse matrix
6345 octave_idx_type sz;
6346 sz = (static_cast<double> (b_nc) - j) / b_nc
6347 * x_nz;
6348 sz = x_nz + (sz > 100 ? sz : 100);
6349 retval.change_capacity (sz);
6350 x_nz = sz;
6351 }
6352 retval.xdata (ii) = tmp;
6353 retval.xridx (ii++) = i;
6354 }
6355 }
6356 retval.xcidx (j+1) = ii;
6357 }
6358
6359 retval.maybe_compress ();
6360
6361 UMFPACK_DNAME (report_info) (control, info);
6362
6363 UMFPACK_DNAME (free_numeric) (&Numeric);
6364 }
6365 else
6366 mattype.mark_as_rectangular ();
6367
6368#else
6369 octave_unused_parameter (rcond);
6370 octave_unused_parameter (sing_handler);
6371 octave_unused_parameter (calc_cond);
6372
6373 (*current_liboctave_error_handler)
6374 ("support for UMFPACK was unavailable or disabled "
6375 "when liboctave was built");
6376#endif
6377 }
6378 else if (typ != MatrixType::Hermitian)
6379 (*current_liboctave_error_handler) ("incorrect matrix type");
6380 }
6381
6382 return retval;
6383}
6384
6386SparseMatrix::fsolve (MatrixType& mattype, const ComplexMatrix& b,
6387 octave_idx_type& err, double& rcond,
6388 solve_singularity_handler sing_handler,
6389 bool calc_cond) const
6390{
6391 ComplexMatrix retval;
6392
6393 octave_idx_type nr = rows ();
6394 octave_idx_type nc = cols ();
6395 err = 0;
6396
6397 if (nr != nc || nr != b.rows ())
6398 (*current_liboctave_error_handler)
6399 ("matrix dimension mismatch solution of linear equations");
6400
6401 if (nr == 0 || b.cols () == 0)
6402 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
6403 else
6404 {
6405 // Print spparms("spumoni") info if requested
6406 int typ = mattype.type ();
6407 mattype.info ();
6408
6409 if (typ == MatrixType::Hermitian)
6410 {
6411#if defined (HAVE_CHOLMOD)
6412 cholmod_common Common;
6413 cholmod_common *cm = &Common;
6414
6415 // Setup initial parameters
6416 CHOLMOD_NAME(start) (cm);
6417 cm->prefer_zomplex = false;
6418
6419 double spu = octave::sparse_params::get_key ("spumoni");
6420 if (spu == 0.)
6421 {
6422 cm->print = -1;
6423 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6424 }
6425 else
6426 {
6427 cm->print = static_cast<int> (spu) + 2;
6428 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6429 }
6430
6431 cm->error_handler = &SparseCholError;
6432 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6433 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6434
6435 cm->final_ll = true;
6436
6437 cholmod_sparse Astore;
6438 cholmod_sparse *A = &Astore;
6439 A->nrow = nr;
6440 A->ncol = nc;
6441
6442 A->p = cidx ();
6443 A->i = ridx ();
6444 A->nzmax = nnz ();
6445 A->packed = true;
6446 A->sorted = true;
6447 A->nz = nullptr;
6448#if defined (OCTAVE_ENABLE_64)
6449 A->itype = CHOLMOD_LONG;
6450#else
6451 A->itype = CHOLMOD_INT;
6452#endif
6453 A->dtype = CHOLMOD_DOUBLE;
6454 A->stype = 1;
6455 A->xtype = CHOLMOD_REAL;
6456
6457 A->x = data ();
6458
6459 cholmod_dense Bstore;
6460 cholmod_dense *B = &Bstore;
6461 B->nrow = b.rows ();
6462 B->ncol = b.cols ();
6463 B->d = B->nrow;
6464 B->nzmax = B->nrow * B->ncol;
6465 B->dtype = CHOLMOD_DOUBLE;
6466 B->xtype = CHOLMOD_COMPLEX;
6467
6468 B->x = const_cast<Complex *> (b.data ());
6469
6470 cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6471 CHOLMOD_NAME(factorize) (A, L, cm);
6472 if (calc_cond)
6473 rcond = CHOLMOD_NAME(rcond)(L, cm);
6474 else
6475 rcond = 1.0;
6476
6477 if (rcond == 0.0)
6478 {
6479 // Either its indefinite or singular. Try UMFPACK
6480 mattype.mark_as_unsymmetric ();
6481 typ = MatrixType::Full;
6482 }
6483 else
6484 {
6485 if (is_singular (rcond) || octave::math::isnan (rcond))
6486 {
6487 err = -2;
6488
6489 if (sing_handler)
6490 {
6491 sing_handler (rcond);
6492 mattype.mark_as_rectangular ();
6493 }
6494 else
6495 octave::warn_singular_matrix (rcond);
6496
6497 return retval;
6498 }
6499
6500 cholmod_dense *X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
6501
6502 retval.resize (b.rows (), b.cols ());
6503 for (octave_idx_type j = 0; j < b.cols (); j++)
6504 {
6505 octave_idx_type jr = j * b.rows ();
6506 for (octave_idx_type i = 0; i < b.rows (); i++)
6507 retval.xelem (i, j) = static_cast<Complex *> (X->x)[jr + i];
6508 }
6509
6510 CHOLMOD_NAME(free_dense) (&X, cm);
6511 CHOLMOD_NAME(free_factor) (&L, cm);
6512 CHOLMOD_NAME(finish) (cm);
6513 static char blank_name[] = " ";
6514 CHOLMOD_NAME(print_common) (blank_name, cm);
6515 }
6516#else
6517 (*current_liboctave_warning_with_id_handler)
6518 ("Octave:missing-dependency",
6519 "support for CHOLMOD was unavailable or disabled "
6520 "when liboctave was built");
6521
6522 mattype.mark_as_unsymmetric ();
6523 typ = MatrixType::Full;
6524#endif
6525 }
6526
6527 if (typ == MatrixType::Full)
6528 {
6529#if defined (HAVE_UMFPACK)
6530 Matrix Control, Info;
6531 void *Numeric = factorize (err, rcond, Control, Info,
6532 sing_handler, calc_cond);
6533
6534 if (err == 0)
6535 {
6536 // one iterative refinement instead of the default two in UMFPACK
6537 Control (UMFPACK_IRSTEP) = 1;
6538 octave_idx_type b_nr = b.rows ();
6539 octave_idx_type b_nc = b.cols ();
6540 int status = 0;
6541 double *control = Control.rwdata ();
6542 double *info = Info.rwdata ();
6543 const octave_idx_type *Ap = cidx ();
6544 const octave_idx_type *Ai = ridx ();
6545 const double *Ax = data ();
6546
6547 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6548 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6549
6550 retval.resize (b_nr, b_nc);
6551
6552 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6553 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6554
6555 for (octave_idx_type j = 0; j < b_nc; j++)
6556 {
6557 for (octave_idx_type i = 0; i < b_nr; i++)
6558 {
6559 Complex c = b(i, j);
6560 Bx[i] = c.real ();
6561 Bz[i] = c.imag ();
6562 }
6563
6564 status = UMFPACK_DNAME (solve) (UMFPACK_A,
6565 octave::to_suitesparse_intptr (Ap),
6566 octave::to_suitesparse_intptr (Ai),
6567 Ax, Xx, Bx, Numeric,
6568 control, info);
6569 int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6570 octave::to_suitesparse_intptr (Ap),
6571 octave::to_suitesparse_intptr (Ai),
6572 Ax, Xz, Bz,
6573 Numeric, control, info);
6574
6575 if (status < 0 || status2 < 0)
6576 {
6577 UMFPACK_DNAME (report_status) (control, status);
6578
6579 // FIXME: Should this be a warning?
6580 (*current_liboctave_error_handler)
6581 ("SparseMatrix::solve solve failed");
6582
6583 err = -1;
6584 break;
6585 }
6586
6587 for (octave_idx_type i = 0; i < b_nr; i++)
6588 retval(i, j) = Complex (Xx[i], Xz[i]);
6589 }
6590
6591 UMFPACK_DNAME (report_info) (control, info);
6592
6593 UMFPACK_DNAME (free_numeric) (&Numeric);
6594 }
6595 else
6596 mattype.mark_as_rectangular ();
6597
6598#else
6599 octave_unused_parameter (rcond);
6600 octave_unused_parameter (sing_handler);
6601 octave_unused_parameter (calc_cond);
6602
6603 (*current_liboctave_error_handler)
6604 ("support for UMFPACK was unavailable or disabled "
6605 "when liboctave was built");
6606#endif
6607 }
6608 else if (typ != MatrixType::Hermitian)
6609 (*current_liboctave_error_handler) ("incorrect matrix type");
6610 }
6611
6612 return retval;
6613}
6614
6616SparseMatrix::fsolve (MatrixType& mattype, const SparseComplexMatrix& b,
6617 octave_idx_type& err, double& rcond,
6618 solve_singularity_handler sing_handler,
6619 bool calc_cond) const
6620{
6621 SparseComplexMatrix retval;
6622
6623 octave_idx_type nr = rows ();
6624 octave_idx_type nc = cols ();
6625 err = 0;
6626
6627 if (nr != nc || nr != b.rows ())
6628 (*current_liboctave_error_handler)
6629 ("matrix dimension mismatch solution of linear equations");
6630
6631 if (nr == 0 || b.cols () == 0)
6632 retval = SparseComplexMatrix (nc, b.cols ());
6633 else
6634 {
6635 // Print spparms("spumoni") info if requested
6636 int typ = mattype.type ();
6637 mattype.info ();
6638
6639 if (typ == MatrixType::Hermitian)
6640 {
6641#if defined (HAVE_CHOLMOD)
6642 cholmod_common Common;
6643 cholmod_common *cm = &Common;
6644
6645 // Setup initial parameters
6646 CHOLMOD_NAME(start) (cm);
6647 cm->prefer_zomplex = false;
6648
6649 double spu = octave::sparse_params::get_key ("spumoni");
6650 if (spu == 0.)
6651 {
6652 cm->print = -1;
6653 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6654 }
6655 else
6656 {
6657 cm->print = static_cast<int> (spu) + 2;
6658 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6659 }
6660
6661 cm->error_handler = &SparseCholError;
6662 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6663 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6664
6665 cm->final_ll = true;
6666
6667 cholmod_sparse Astore;
6668 cholmod_sparse *A = &Astore;
6669 A->nrow = nr;
6670 A->ncol = nc;
6671
6672 A->p = cidx ();
6673 A->i = ridx ();
6674 A->nzmax = nnz ();
6675 A->packed = true;
6676 A->sorted = true;
6677 A->nz = nullptr;
6678#if defined (OCTAVE_ENABLE_64)
6679 A->itype = CHOLMOD_LONG;
6680#else
6681 A->itype = CHOLMOD_INT;
6682#endif
6683 A->dtype = CHOLMOD_DOUBLE;
6684 A->stype = 1;
6685 A->xtype = CHOLMOD_REAL;
6686
6687 A->x = data ();
6688
6689 cholmod_sparse Bstore;
6690 cholmod_sparse *B = &Bstore;
6691 B->nrow = b.rows ();
6692 B->ncol = b.cols ();
6693 B->p = b.cidx ();
6694 B->i = b.ridx ();
6695 B->nzmax = b.nnz ();
6696 B->packed = true;
6697 B->sorted = true;
6698 B->nz = nullptr;
6699#if defined (OCTAVE_ENABLE_64)
6700 B->itype = CHOLMOD_LONG;
6701#else
6702 B->itype = CHOLMOD_INT;
6703#endif
6704 B->dtype = CHOLMOD_DOUBLE;
6705 B->stype = 0;
6706 B->xtype = CHOLMOD_COMPLEX;
6707
6708 B->x = b.data ();
6709
6710 cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6711 CHOLMOD_NAME(factorize) (A, L, cm);
6712 if (calc_cond)
6713 rcond = CHOLMOD_NAME(rcond)(L, cm);
6714 else
6715 rcond = 1.0;
6716
6717 if (rcond == 0.0)
6718 {
6719 // Either its indefinite or singular. Try UMFPACK
6720 mattype.mark_as_unsymmetric ();
6721 typ = MatrixType::Full;
6722 }
6723 else
6724 {
6725 if (is_singular (rcond) || octave::math::isnan (rcond))
6726 {
6727 err = -2;
6728
6729 if (sing_handler)
6730 {
6731 sing_handler (rcond);
6732 mattype.mark_as_rectangular ();
6733 }
6734 else
6735 octave::warn_singular_matrix (rcond);
6736
6737 return retval;
6738 }
6739
6740 cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6741
6742 retval = SparseComplexMatrix
6743 (octave::from_size_t (X->nrow),
6744 octave::from_size_t (X->ncol),
6745 octave::from_suitesparse_long (CHOLMOD_NAME(nnz)
6746 (X, cm)));
6747 for (octave_idx_type j = 0;
6748 j <= static_cast<octave_idx_type> (X->ncol); j++)
6749 retval.xcidx (j) = static_cast<octave_idx_type *> (X->p)[j];
6750 for (octave_idx_type j = 0;
6751 j < octave::from_suitesparse_long (CHOLMOD_NAME(nnz) (X, cm));
6752 j++)
6753 {
6754 retval.xridx (j) = static_cast<octave_idx_type *> (X->i)[j];
6755 retval.xdata (j) = static_cast<Complex *> (X->x)[j];
6756 }
6757
6758 CHOLMOD_NAME(free_sparse) (&X, cm);
6759 CHOLMOD_NAME(free_factor) (&L, cm);
6760 CHOLMOD_NAME(finish) (cm);
6761 static char blank_name[] = " ";
6762 CHOLMOD_NAME(print_common) (blank_name, cm);
6763 }
6764#else
6765 (*current_liboctave_warning_with_id_handler)
6766 ("Octave:missing-dependency",
6767 "support for CHOLMOD was unavailable or disabled "
6768 "when liboctave was built");
6769
6770 mattype.mark_as_unsymmetric ();
6771 typ = MatrixType::Full;
6772#endif
6773 }
6774
6775 if (typ == MatrixType::Full)
6776 {
6777#if defined (HAVE_UMFPACK)
6778 Matrix Control, Info;
6779 void *Numeric = factorize (err, rcond, Control, Info,
6780 sing_handler, calc_cond);
6781
6782 if (err == 0)
6783 {
6784 // one iterative refinement instead of the default two in UMFPACK
6785 Control (UMFPACK_IRSTEP) = 1;
6786 octave_idx_type b_nr = b.rows ();
6787 octave_idx_type b_nc = b.cols ();
6788 int status = 0;
6789 double *control = Control.rwdata ();
6790 double *info = Info.rwdata ();
6791 const octave_idx_type *Ap = cidx ();
6792 const octave_idx_type *Ai = ridx ();
6793 const double *Ax = data ();
6794
6795 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6796 OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6797
6798 // Take a first guess that the number of nonzero terms
6799 // will be as many as in b
6800 octave_idx_type x_nz = b.nnz ();
6801 octave_idx_type ii = 0;
6802 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6803
6804 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6805 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6806
6807 retval.xcidx (0) = 0;
6808 for (octave_idx_type j = 0; j < b_nc; j++)
6809 {
6810 for (octave_idx_type i = 0; i < b_nr; i++)
6811 {
6812 Complex c = b(i, j);
6813 Bx[i] = c.real ();
6814 Bz[i] = c.imag ();
6815 }
6816
6817 status = UMFPACK_DNAME (solve) (UMFPACK_A,
6818 octave::to_suitesparse_intptr (Ap),
6819 octave::to_suitesparse_intptr (Ai),
6820 Ax, Xx, Bx, Numeric,
6821 control, info);
6822 int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6823 octave::to_suitesparse_intptr (Ap),
6824 octave::to_suitesparse_intptr (Ai),
6825 Ax, Xz, Bz,
6826 Numeric, control, info);
6827
6828 if (status < 0 || status2 < 0)
6829 {
6830 UMFPACK_DNAME (report_status) (control, status);
6831
6832 // FIXME: Should this be a warning?
6833 (*current_liboctave_error_handler)
6834 ("SparseMatrix::solve solve failed");
6835
6836 err = -1;
6837 break;
6838 }
6839
6840 for (octave_idx_type i = 0; i < b_nr; i++)
6841 {
6842 Complex tmp = Complex (Xx[i], Xz[i]);
6843 if (tmp != 0.0)
6844 {
6845 if (ii == x_nz)
6846 {
6847 // Resize the sparse matrix
6848 octave_idx_type sz;
6849 sz = (static_cast<double> (b_nc) - j) / b_nc
6850 * x_nz;
6851 sz = x_nz + (sz > 100 ? sz : 100);
6852 retval.change_capacity (sz);
6853 x_nz = sz;
6854 }
6855 retval.xdata (ii) = tmp;
6856 retval.xridx (ii++) = i;
6857 }
6858 }
6859 retval.xcidx (j+1) = ii;
6860 }
6861
6862 retval.maybe_compress ();
6863
6864 UMFPACK_DNAME (report_info) (control, info);
6865
6866 UMFPACK_DNAME (free_numeric) (&Numeric);
6867 }
6868 else
6869 mattype.mark_as_rectangular ();
6870#else
6871 octave_unused_parameter (rcond);
6872 octave_unused_parameter (sing_handler);
6873 octave_unused_parameter (calc_cond);
6874
6875 (*current_liboctave_error_handler)
6876 ("support for UMFPACK was unavailable or disabled "
6877 "when liboctave was built");
6878#endif
6879 }
6880 else if (typ != MatrixType::Hermitian)
6881 (*current_liboctave_error_handler) ("incorrect matrix type");
6882 }
6883
6884 return retval;
6885}
6886
6887Matrix
6888SparseMatrix::solve (MatrixType& mattype, const Matrix& b) const
6889{
6890 octave_idx_type info;
6891 double rcond;
6892 return solve (mattype, b, info, rcond, nullptr);
6893}
6894
6895Matrix
6897 octave_idx_type& info) const
6898{
6899 double rcond;
6900 return solve (mattype, b, info, rcond, nullptr);
6901}
6902
6903Matrix
6905 octave_idx_type& info, double& rcond) const
6906{
6907 return solve (mattype, b, info, rcond, nullptr);
6908}
6909
6910Matrix
6912 double& rcond, solve_singularity_handler sing_handler,
6913 bool singular_fallback) const
6914{
6915 Matrix retval;
6916 int typ = mattype.type (false);
6917
6918 if (typ == MatrixType::Unknown)
6919 typ = mattype.type (*this);
6920
6921 // Only calculate the condition number for CHOLMOD/UMFPACK
6923 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6924 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6925 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6926 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6927 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6928 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6929 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6930 else if (typ == MatrixType::Tridiagonal
6932 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6933 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6934 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6935 else if (typ != MatrixType::Rectangular)
6936 (*current_liboctave_error_handler) ("unknown matrix type");
6937
6938 // Rectangular or one of the above solvers flags a singular matrix
6939 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6940 {
6941 rcond = 1.;
6942 retval = dmsolve<Matrix, SparseMatrix, Matrix> (*this, b, err);
6943 }
6944
6945 return retval;
6946}
6947
6950{
6951 octave_idx_type info;
6952 double rcond;
6953 return solve (mattype, b, info, rcond, nullptr);
6954}
6955
6958 octave_idx_type& info) const
6959{
6960 double rcond;
6961 return solve (mattype, b, info, rcond, nullptr);
6962}
6963
6966 octave_idx_type& info, double& rcond) const
6967{
6968 return solve (mattype, b, info, rcond, nullptr);
6969}
6970
6973 octave_idx_type& err, double& rcond,
6974 solve_singularity_handler sing_handler,
6975 bool singular_fallback) const
6976{
6977 SparseMatrix retval;
6978 int typ = mattype.type (false);
6979
6980 if (typ == MatrixType::Unknown)
6981 typ = mattype.type (*this);
6982
6984 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6985 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6986 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6987 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6988 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6989 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6990 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6991 else if (typ == MatrixType::Tridiagonal
6993 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6994 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6995 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6996 else if (typ != MatrixType::Rectangular)
6997 (*current_liboctave_error_handler) ("unknown matrix type");
6998
6999 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
7000 {
7001 rcond = 1.;
7003 (*this, b, err);
7004 }
7005
7006 return retval;
7007}
7008
7011{
7012 octave_idx_type info;
7013 double rcond;
7014 return solve (mattype, b, info, rcond, nullptr);
7015}
7016
7019 octave_idx_type& info) const
7020{
7021 double rcond;
7022 return solve (mattype, b, info, rcond, nullptr);
7023}
7024
7027 octave_idx_type& info, double& rcond) const
7028{
7029 return solve (mattype, b, info, rcond, nullptr);
7030}
7031
7034 octave_idx_type& err, double& rcond,
7035 solve_singularity_handler sing_handler,
7036 bool singular_fallback) const
7037{
7038 ComplexMatrix retval;
7039 int typ = mattype.type (false);
7040
7041 if (typ == MatrixType::Unknown)
7042 typ = mattype.type (*this);
7043
7045 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
7046 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
7047 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
7048 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
7049 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
7050 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
7051 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
7052 else if (typ == MatrixType::Tridiagonal
7054 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
7055 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
7056 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
7057 else if (typ != MatrixType::Rectangular)
7058 (*current_liboctave_error_handler) ("unknown matrix type");
7059
7060 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
7061 {
7062 rcond = 1.;
7064 (*this, b, err);
7065 }
7066
7067 return retval;
7068}
7069
7072{
7073 octave_idx_type info;
7074 double rcond;
7075 return solve (mattype, b, info, rcond, nullptr);
7076}
7077
7080 octave_idx_type& info) const
7081{
7082 double rcond;
7083 return solve (mattype, b, info, rcond, nullptr);
7084}
7085
7088 octave_idx_type& info, double& rcond) const
7089{
7090 return solve (mattype, b, info, rcond, nullptr);
7091}
7092
7095 octave_idx_type& err, double& rcond,
7096 solve_singularity_handler sing_handler,
7097 bool singular_fallback) const
7098{
7099 SparseComplexMatrix retval;
7100 int typ = mattype.type (false);
7101
7102 if (typ == MatrixType::Unknown)
7103 typ = mattype.type (*this);
7104
7106 retval = dsolve (mattype, b, err, rcond, sing_handler, false);
7107 else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
7108 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
7109 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
7110 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
7111 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
7112 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
7113 else if (typ == MatrixType::Tridiagonal
7115 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
7116 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
7117 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
7118 else if (typ != MatrixType::Rectangular)
7119 (*current_liboctave_error_handler) ("unknown matrix type");
7120
7121 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
7122 {
7123 rcond = 1.;
7125 (*this, b, err);
7126 }
7127
7128 return retval;
7129}
7130
7133{
7134 octave_idx_type info; double rcond;
7135 return solve (mattype, b, info, rcond);
7136}
7137
7140 octave_idx_type& info) const
7141{
7142 double rcond;
7143 return solve (mattype, b, info, rcond);
7144}
7145
7148 octave_idx_type& info, double& rcond) const
7149{
7150 return solve (mattype, b, info, rcond, nullptr);
7151}
7152
7155 octave_idx_type& info, double& rcond,
7156 solve_singularity_handler sing_handler) const
7157{
7158 Matrix tmp (b);
7159 return solve (mattype, tmp, info, rcond,
7160 sing_handler).column (static_cast<octave_idx_type> (0));
7161}
7162
7165{
7166 octave_idx_type info;
7167 double rcond;
7168 return solve (mattype, b, info, rcond, nullptr);
7169}
7170
7173 octave_idx_type& info) const
7174{
7175 double rcond;
7176 return solve (mattype, b, info, rcond, nullptr);
7177}
7178
7181 octave_idx_type& info,
7182 double& rcond) const
7183{
7184 return solve (mattype, b, info, rcond, nullptr);
7185}
7186
7189 octave_idx_type& info, double& rcond,
7190 solve_singularity_handler sing_handler) const
7191{
7192 ComplexMatrix tmp (b);
7193 return solve (mattype, tmp, info, rcond,
7194 sing_handler).column (static_cast<octave_idx_type> (0));
7195}
7196
7197Matrix
7199{
7200 octave_idx_type info;
7201 double rcond;
7202 return solve (b, info, rcond, nullptr);
7203}
7204
7205Matrix
7207{
7208 double rcond;
7209 return solve (b, info, rcond, nullptr);
7210}
7211
7212Matrix
7214 double& rcond) const
7215{
7216 return solve (b, info, rcond, nullptr);
7217}
7218
7219Matrix
7220SparseMatrix::solve (const Matrix& b, octave_idx_type& err, double& rcond,
7221 solve_singularity_handler sing_handler) const
7222{
7223 MatrixType mattype (*this);
7224 return solve (mattype, b, err, rcond, sing_handler);
7225}
7226
7229{
7230 octave_idx_type info;
7231 double rcond;
7232 return solve (b, info, rcond, nullptr);
7233}
7234
7237 octave_idx_type& info) const
7238{
7239 double rcond;
7240 return solve (b, info, rcond, nullptr);
7241}
7242
7245 octave_idx_type& info, double& rcond) const
7246{
7247 return solve (b, info, rcond, nullptr);
7248}
7249
7251SparseMatrix::solve (const SparseMatrix& b, octave_idx_type& err, double& rcond,
7252 solve_singularity_handler sing_handler) const
7253{
7254 MatrixType mattype (*this);
7255 return solve (mattype, b, err, rcond, sing_handler);
7256}
7257
7260{
7261 double rcond;
7262 return solve (b, info, rcond, nullptr);
7263}
7264
7267 double& rcond) const
7268{
7269 return solve (b, info, rcond, nullptr);
7270}
7271
7274 double& rcond,
7275 solve_singularity_handler sing_handler) const
7276{
7277 MatrixType mattype (*this);
7278 return solve (mattype, b, err, rcond, sing_handler);
7279}
7280
7283{
7284 octave_idx_type info;
7285 double rcond;
7286 return solve (b, info, rcond, nullptr);
7287}
7288
7291{
7292 double rcond;
7293 return solve (b, info, rcond, nullptr);
7294}
7295
7298 double& rcond) const
7299{
7300 return solve (b, info, rcond, nullptr);
7301}
7302
7305 octave_idx_type& err, double& rcond,
7306 solve_singularity_handler sing_handler) const
7307{
7308 MatrixType mattype (*this);
7309 return solve (mattype, b, err, rcond, sing_handler);
7310}
7311
7314{
7315 octave_idx_type info; double rcond;
7316 return solve (b, info, rcond);
7317}
7318
7321{
7322 double rcond;
7323 return solve (b, info, rcond);
7324}
7325
7328 double& rcond) const
7329{
7330 return solve (b, info, rcond, nullptr);
7331}
7332
7335 double& rcond,
7336 solve_singularity_handler sing_handler) const
7337{
7338 Matrix tmp (b);
7339 return solve (tmp, info, rcond,
7340 sing_handler).column (static_cast<octave_idx_type> (0));
7341}
7342
7345{
7346 octave_idx_type info;
7347 double rcond;
7348 return solve (b, info, rcond, nullptr);
7349}
7350
7353{
7354 double rcond;
7355 return solve (b, info, rcond, nullptr);
7356}
7357
7360 double& rcond) const
7361{
7362 return solve (b, info, rcond, nullptr);
7363}
7364
7367 double& rcond,
7368 solve_singularity_handler sing_handler) const
7369{
7370 ComplexMatrix tmp (b);
7371 return solve (tmp, info, rcond,
7372 sing_handler).column (static_cast<octave_idx_type> (0));
7373}
7374
7375// other operations.
7376
7377bool
7379{
7380 octave_idx_type nel = nnz ();
7381
7382 if (neg_zero)
7383 {
7384 for (octave_idx_type i = 0; i < nel; i++)
7385 if (octave::math::signbit (data (i)))
7386 return true;
7387 }
7388 else
7389 {
7390 for (octave_idx_type i = 0; i < nel; i++)
7391 if (data (i) < 0)
7392 return true;
7393 }
7394
7395 return false;
7396}
7397
7398bool
7400{
7401 octave_idx_type nel = nnz ();
7402
7403 for (octave_idx_type i = 0; i < nel; i++)
7404 {
7405 double val = data (i);
7406 if (octave::math::isnan (val))
7407 return true;
7408 }
7409
7410 return false;
7411}
7412
7413bool
7415{
7416 octave_idx_type nel = nnz ();
7417
7418 for (octave_idx_type i = 0; i < nel; i++)
7419 {
7420 double val = data (i);
7421 if (octave::math::isinf (val) || octave::math::isnan (val))
7422 return true;
7423 }
7424
7425 return false;
7426}
7427
7428bool
7430{
7431 octave_idx_type nel = nnz ();
7432
7433 for (octave_idx_type i = 0; i < nel; i++)
7434 {
7435 double val = data (i);
7436 if (val != 0.0 && val != 1.0)
7437 return true;
7438 }
7439
7440 return false;
7441}
7442
7443bool
7445{
7446 octave_idx_type nel = nnz ();
7447
7448 for (octave_idx_type i = 0; i < nel; i++)
7449 if (data (i) != 0)
7450 return false;
7451
7452 return true;
7453}
7454
7455bool
7457{
7458 octave_idx_type nel = nnz ();
7459
7460 for (octave_idx_type i = 0; i < nel; i++)
7461 {
7462 double val = data (i);
7463 if (octave::math::isnan (val) || octave::math::is_integer (val))
7464 continue;
7465 else
7466 return false;
7467 }
7468
7469 return true;
7470}
7471
7472// Return nonzero if any element of M is not an integer. Also extract
7473// the largest and smallest values and return them in MAX_VAL and MIN_VAL.
7474
7475bool
7476SparseMatrix::all_integers (double& max_val, double& min_val) const
7477{
7478 octave_idx_type nel = nnz ();
7479
7480 if (nel == 0)
7481 return false;
7482
7483 max_val = data (0);
7484 min_val = data (0);
7485
7486 for (octave_idx_type i = 0; i < nel; i++)
7487 {
7488 double val = data (i);
7489
7490 if (val > max_val)
7491 max_val = val;
7492
7493 if (val < min_val)
7494 min_val = val;
7495
7496 if (! octave::math::is_integer (val))
7497 return false;
7498 }
7499
7500 return true;
7501}
7502
7503bool
7505{
7506 return test_any (octave::too_large_for_float);
7507}
7508
7511{
7512 if (any_element_is_nan ())
7513 octave::err_nan_to_logical_conversion ();
7514
7515 octave_idx_type nr = rows ();
7516 octave_idx_type nc = cols ();
7517 octave_idx_type nz1 = nnz ();
7518 octave_idx_type nz2 = nr*nc - nz1;
7519
7520 SparseBoolMatrix r (nr, nc, nz2);
7521
7522 octave_idx_type ii = 0;
7523 octave_idx_type jj = 0;
7524 r.cidx (0) = 0;
7525 for (octave_idx_type i = 0; i < nc; i++)
7526 {
7527 for (octave_idx_type j = 0; j < nr; j++)
7528 {
7529 if (jj < cidx (i+1) && ridx (jj) == j)
7530 jj++;
7531 else
7532 {
7533 r.data (ii) = true;
7534 r.ridx (ii++) = j;
7535 }
7536 }
7537 r.cidx (i+1) = ii;
7538 }
7539
7540 return r;
7541}
7542
7543// FIXME: Do these really belong here? Maybe they should be in a base class?
7544
7546SparseMatrix::all (int dim) const
7547{
7549 SPARSE_ALL_OP (dim);
7550}
7551
7553SparseMatrix::any (int dim) const
7554{
7556 SPARSE_ANY_OP (dim);
7557}
7558
7560SparseMatrix::cumprod (int dim, bool nanflag) const
7561{
7562 if (dim > 1)
7563 return *this;
7565}
7566
7568SparseMatrix::cumsum (int dim, bool nanflag) const
7569{
7570#define NAN_EXPR \
7571 if (! octave::math::isnan (data (j))) \
7572 t += data (j); \
7573 else \
7574 t += 0.0
7575
7576#define EXPR \
7577 t += data (j)
7578
7579 if (dim > 1)
7580 return *this;
7582
7583#undef NAN_EXPR
7584#undef EXPR
7585}
7586
7588SparseMatrix::prod (int dim, bool nanflag) const
7589{
7590 if (dim > 1)
7591 return *this;
7592 else if ((rows () == 1 && dim == -1) || dim == 1)
7593 return transpose ().prod (0, nanflag).transpose ();
7594 else
7595 {
7596 #define ROW_EXPR \
7597 double d = data (i); \
7598 if (nanflag && octave::math::isnan (d)) \
7599 tmp[ridx (i)] *= 1.0; \
7600 else \
7601 tmp[ridx (i)] *= d
7602
7603 #define COL_EXPR \
7604 double d = data (i); \
7605 if (nanflag && octave::math::isnan (d)) \
7606 tmp[j] *= 1.0; \
7607 else \
7608 tmp[j] *= d
7609
7611 (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
7612
7613 #undef ROW_EXPR
7614 #undef COL_EXPR
7615 }
7616}
7617
7619SparseMatrix::sum (int dim, bool nanflag) const
7620{
7621#define ROW_EXPR \
7622 double d = data (i); \
7623 if (nanflag && octave::math::isnan (d)) \
7624 tmp[ridx (i)] += 0.0; \
7625 else \
7626 tmp[ridx (i)] += d
7627
7628#define COL_EXPR \
7629 double d = data (i); \
7630 if (nanflag && octave::math::isnan (d)) \
7631 tmp[j] += 0.0; \
7632 else \
7633 tmp[j] += d
7634
7635 if (dim > 1)
7636 return *this;
7638 0.0, 0.0);
7639
7640#undef ROW_EXPR
7641#undef COL_EXPR
7642}
7643
7645SparseMatrix::xsum (int dim, bool nanflag) const
7646{
7647 if (dim > 1)
7648 return *this;
7650}
7651
7653SparseMatrix::sumsq (int dim, bool nanflag) const
7654{
7655#define EXPR r.data (nel++) = data (i) * data (i);
7656
7657#define ROW_EXPR \
7658 double d = data (i); \
7659 if (nanflag && octave::math::isnan (d)) \
7660 tmp[ridx (i)] += 0.0; \
7661 else \
7662 tmp[ridx (i)] += d * d
7663
7664#define COL_EXPR \
7665 double d = data (i); \
7666 if (nanflag && octave::math::isnan (d)) \
7667 tmp[j] += 0.0; \
7668 else \
7669 tmp[j] += d * d
7670
7673 0.0, 0.0);
7674
7675#undef EXPR
7676#undef ROW_EXPR
7677#undef COL_EXPR
7678}
7679
7682{
7683 octave_idx_type nz = nnz ();
7684
7685 SparseMatrix retval (*this);
7686
7687 for (octave_idx_type i = 0; i < nz; i++)
7688 retval.data (i) = fabs (retval.data (i));
7689
7690 return retval;
7691}
7692
7695{
7696 return MSparse<double>::diag (k);
7697}
7698
7699Matrix
7701{
7703}
7704
7705std::ostream&
7706operator << (std::ostream& os, const SparseMatrix& a)
7707{
7708 octave_idx_type nc = a.cols ();
7709
7710 // add one to the printed indices to go from
7711 // zero-based to one-based arrays
7712 for (octave_idx_type j = 0; j < nc; j++)
7713 {
7714 octave_quit ();
7715 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7716 {
7717 os << a.ridx (i) + 1 << ' ' << j + 1 << ' ';
7718 octave::write_value<double> (os, a.data (i));
7719 os << "\n";
7720 }
7721 }
7722
7723 return os;
7724}
7725
7726std::istream&
7727operator >> (std::istream& is, SparseMatrix& a)
7728{
7729 typedef SparseMatrix::element_type elt_type;
7730
7731 return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7732}
7733
7736{
7737 return MSparse<double>::squeeze ();
7738}
7739
7741SparseMatrix::reshape (const dim_vector& new_dims) const
7742{
7743 return MSparse<double>::reshape (new_dims);
7744}
7745
7748{
7749 return MSparse<double>::permute (vec, inv);
7750}
7751
7754{
7755 return MSparse<double>::ipermute (vec);
7756}
7757
7758// matrix by matrix -> matrix operations
7759
7762{
7763 SPARSE_SPARSE_MUL (SparseMatrix, double, double);
7764}
7765
7766Matrix
7768{
7769 FULL_SPARSE_MUL (Matrix, double);
7770}
7771
7772Matrix
7773mul_trans (const Matrix& m, const SparseMatrix& a)
7774{
7775 FULL_SPARSE_MUL_TRANS (Matrix, double, );
7776}
7777
7778Matrix
7780{
7781 SPARSE_FULL_MUL (Matrix, double);
7782}
7783
7784Matrix
7785trans_mul (const SparseMatrix& m, const Matrix& a)
7786{
7787 SPARSE_FULL_TRANS_MUL (Matrix, double, );
7788}
7789
7790// diag * sparse and sparse * diag
7791
7794{
7795 return do_mul_dm_sm<SparseMatrix> (d, a);
7796}
7797
7800{
7801 return do_mul_sm_dm<SparseMatrix> (a, d);
7802}
7803
7806{
7807 return do_add_dm_sm<SparseMatrix> (d, a);
7808}
7809
7812{
7813 return do_sub_dm_sm<SparseMatrix> (d, a);
7814}
7815
7818{
7819 return do_add_sm_dm<SparseMatrix> (a, d);
7820}
7821
7824{
7825 return do_sub_sm_dm<SparseMatrix> (a, d);
7826}
7827
7828// perm * sparse and sparse * perm
7829
7832{
7833 return octinternal_do_mul_pm_sm (p, a);
7834}
7835
7838{
7839 return octinternal_do_mul_sm_pm (a, p);
7840}
7841
7842// FIXME: it would be nice to share code among the min/max functions below.
7843
7844#define EMPTY_RETURN_CHECK(T) \
7845 if (nr == 0 || nc == 0) \
7846 return T (nr, nc);
7847
7849min (double d, const SparseMatrix& m)
7850{
7851 const bool nanflag = true;
7852 const bool realabs = true;
7853 return min (d, m, nanflag, realabs);
7854}
7855
7857min (double d, const SparseMatrix& m, const bool nanflag)
7858{
7859 const bool realabs = true;
7860 return min (d, m, nanflag, realabs);
7861}
7862
7864min (double d, const SparseMatrix& m, const bool nanflag, const bool realabs)
7865{
7866 SparseMatrix result;
7867
7868 octave_idx_type nr = m.rows ();
7869 octave_idx_type nc = m.columns ();
7870
7872
7873 // Count the number of nonzero elements
7874 if (d < 0.)
7875 {
7876 result = SparseMatrix (nr, nc, d);
7877 for (octave_idx_type j = 0; j < nc; j++)
7878 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7879 {
7880 double tmp = octave::math::min (d, m.data (i), nanflag, realabs);
7881 if (tmp != 0.)
7882 {
7883 octave_idx_type idx = m.ridx (i) + j * nr;
7884 result.xdata (idx) = tmp;
7885 result.xridx (idx) = m.ridx (i);
7886 }
7887 }
7888 }
7889 else
7890 {
7891 octave_idx_type nel = 0;
7892 for (octave_idx_type j = 0; j < nc; j++)
7893 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7894 if (octave::math::min (d, m.data (i), nanflag, realabs) != 0.)
7895 nel++;
7896
7897 result = SparseMatrix (nr, nc, nel);
7898
7899 octave_idx_type ii = 0;
7900 result.xcidx (0) = 0;
7901 for (octave_idx_type j = 0; j < nc; j++)
7902 {
7903 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7904 {
7905 double tmp = octave::math::min (d, m.data (i), nanflag, realabs);
7906
7907 if (tmp != 0.)
7908 {
7909 result.xdata (ii) = tmp;
7910 result.xridx (ii++) = m.ridx (i);
7911 }
7912 }
7913 result.xcidx (j+1) = ii;
7914 }
7915 }
7916
7917 return result;
7918}
7919
7921min (const SparseMatrix& m, double d)
7922{
7923 const bool nanflag = true;
7924 const bool realabs = true;
7925 return min (d, m, nanflag, realabs);
7926}
7927
7929min (const SparseMatrix& m, double d, const bool nanflag)
7930{
7931 const bool realabs = true;
7932 return min (d, m, nanflag, realabs);
7933}
7934
7936min (const SparseMatrix& m, double d, const bool nanflag, const bool realabs)
7937{
7938 return min (d, m, nanflag, realabs);
7939}
7940
7942min (const SparseMatrix& a, const SparseMatrix& b)
7943{
7944 const bool nanflag = true;
7945 const bool realabs = true;
7946 return min (a, b, nanflag, realabs);
7947}
7948
7950min (const SparseMatrix& a, const SparseMatrix& b, const bool nanflag)
7951{
7952 const bool realabs = true;
7953 return min (a, b, nanflag, realabs);
7954}
7955
7957min (const SparseMatrix& a, const SparseMatrix& b, const bool nanflag,
7958 const bool realabs)
7959{
7960 SparseMatrix r;
7961
7962 octave_idx_type a_nr = a.rows ();
7963 octave_idx_type a_nc = a.cols ();
7964 octave_idx_type b_nr = b.rows ();
7965 octave_idx_type b_nc = b.cols ();
7966
7967 if (a_nr == b_nr && a_nc == b_nc)
7968 {
7969 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7970
7971 octave_idx_type jx = 0;
7972 r.cidx (0) = 0;
7973 for (octave_idx_type i = 0 ; i < a_nc ; i++)
7974 {
7975 octave_idx_type ja = a.cidx (i);
7976 octave_idx_type ja_max = a.cidx (i+1);
7977 bool ja_lt_max = ja < ja_max;
7978
7979 octave_idx_type jb = b.cidx (i);
7980 octave_idx_type jb_max = b.cidx (i+1);
7981 bool jb_lt_max = jb < jb_max;
7982
7983 while (ja_lt_max || jb_lt_max)
7984 {
7985 octave_quit ();
7986 if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7987 {
7988 double tmp = octave::math::min (a.data (ja), 0.,
7989 nanflag, realabs);
7990 if (tmp != 0.)
7991 {
7992 r.ridx (jx) = a.ridx (ja);
7993 r.data (jx) = tmp;
7994 jx++;
7995 }
7996 ja++;
7997 ja_lt_max= ja < ja_max;
7998 }
7999 else if ((! ja_lt_max)
8000 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
8001 {
8002 double tmp = octave::math::min (0., b.data (jb),
8003 nanflag, realabs);
8004 if (tmp != 0.)
8005 {
8006 r.ridx (jx) = b.ridx (jb);
8007 r.data (jx) = tmp;
8008 jx++;
8009 }
8010 jb++;
8011 jb_lt_max= jb < jb_max;
8012 }
8013 else
8014 {
8015 double tmp = octave::math::min (a.data (ja), b.data (jb),
8016 nanflag, realabs);
8017 if (tmp != 0.)
8018 {
8019 r.data (jx) = tmp;
8020 r.ridx (jx) = a.ridx (ja);
8021 jx++;
8022 }
8023 ja++;
8024 ja_lt_max= ja < ja_max;
8025 jb++;
8026 jb_lt_max= jb < jb_max;
8027 }
8028 }
8029 r.cidx (i+1) = jx;
8030 }
8031
8032 r.maybe_compress ();
8033 }
8034 else
8035 {
8036 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8037 {
8038 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8039 r.resize (a_nr, std::max (a_nc, b_nc));
8040 else
8041 octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
8042 }
8043 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8044 {
8045 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8046 r.resize (std::max (a_nr, b_nr), a_nc);
8047 else
8048 octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
8049 }
8050 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8051 {
8052 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8053 r.resize (b_nr, std::max (a_nc, b_nc));
8054 else
8055 octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
8056 }
8057 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8058 {
8059 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8060 r.resize (std::max (a_nr, b_nr), b_nc);
8061 else
8062 octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
8063 }
8064 else
8065 octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
8066 }
8067
8068 return r;
8069}
8070
8072max (double d, const SparseMatrix& m)
8073{
8074 const bool nanflag = true;
8075 const bool realabs = true;
8076 return max (d, m, nanflag, realabs);
8077}
8078
8080max (double d, const SparseMatrix& m, const bool nanflag)
8081{
8082 const bool realabs = true;
8083 return max (d, m, nanflag, realabs);
8084}
8085
8087max (double d, const SparseMatrix& m, const bool nanflag, const bool realabs)
8088{
8089 SparseMatrix result;
8090
8091 octave_idx_type nr = m.rows ();
8092 octave_idx_type nc = m.columns ();
8093
8095
8096 // Count the number of nonzero elements
8097 if (d > 0.)
8098 {
8099 result = SparseMatrix (nr, nc, d);
8100 for (octave_idx_type j = 0; j < nc; j++)
8101 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
8102 {
8103 double tmp = octave::math::max (d, m.data (i), nanflag, realabs);
8104
8105 if (tmp != 0.)
8106 {
8107 octave_idx_type idx = m.ridx (i) + j * nr;
8108 result.xdata (idx) = tmp;
8109 result.xridx (idx) = m.ridx (i);
8110 }
8111 }
8112 }
8113 else
8114 {
8115 octave_idx_type nel = 0;
8116 for (octave_idx_type j = 0; j < nc; j++)
8117 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
8118 if (octave::math::max (d, m.data (i), nanflag, realabs) != 0.)
8119 nel++;
8120
8121 result = SparseMatrix (nr, nc, nel);
8122
8123 octave_idx_type ii = 0;
8124 result.xcidx (0) = 0;
8125 for (octave_idx_type j = 0; j < nc; j++)
8126 {
8127 for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
8128 {
8129 double tmp = octave::math::max (d, m.data (i), nanflag, realabs);
8130 if (tmp != 0.)
8131 {
8132 result.xdata (ii) = tmp;
8133 result.xridx (ii++) = m.ridx (i);
8134 }
8135 }
8136 result.xcidx (j+1) = ii;
8137 }
8138 }
8139
8140 return result;
8141}
8142
8144max (const SparseMatrix& m, double d)
8145{
8146 const bool nanflag = true;
8147 const bool realabs = true;
8148 return max (d, m, nanflag, realabs);
8149}
8150
8152max (const SparseMatrix& m, double d, const bool nanflag)
8153{
8154 const bool realabs = true;
8155 return max (d, m, nanflag, realabs);
8156}
8157
8159max (const SparseMatrix& m, double d, const bool nanflag, const bool realabs)
8160{
8161 return max (d, m, nanflag, realabs);
8162}
8163
8165max (const SparseMatrix& a, const SparseMatrix& b)
8166{
8167 const bool nanflag = true;
8168 const bool realabs = true;
8169 return max (a, b, nanflag, realabs);
8170}
8171
8173max (const SparseMatrix& a, const SparseMatrix& b, const bool nanflag)
8174{
8175 const bool realabs = true;
8176 return max (a, b, nanflag, realabs);
8177}
8178
8180max (const SparseMatrix& a, const SparseMatrix& b, const bool nanflag,
8181 const bool realabs)
8182{
8183 SparseMatrix r;
8184
8185 octave_idx_type a_nr = a.rows ();
8186 octave_idx_type a_nc = a.cols ();
8187 octave_idx_type b_nr = b.rows ();
8188 octave_idx_type b_nc = b.cols ();
8189
8190 if (a_nr == b_nr && a_nc == b_nc)
8191 {
8192 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
8193
8194 octave_idx_type jx = 0;
8195 r.cidx (0) = 0;
8196 for (octave_idx_type i = 0 ; i < a_nc ; i++)
8197 {
8198 octave_idx_type ja = a.cidx (i);
8199 octave_idx_type ja_max = a.cidx (i+1);
8200 bool ja_lt_max = ja < ja_max;
8201
8202 octave_idx_type jb = b.cidx (i);
8203 octave_idx_type jb_max = b.cidx (i+1);
8204 bool jb_lt_max = jb < jb_max;
8205
8206 while (ja_lt_max || jb_lt_max)
8207 {
8208 octave_quit ();
8209 if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
8210 {
8211 double tmp = octave::math::max (a.data (ja), 0.,
8212 nanflag, realabs);
8213 if (tmp != 0.)
8214 {
8215 r.ridx (jx) = a.ridx (ja);
8216 r.data (jx) = tmp;
8217 jx++;
8218 }
8219 ja++;
8220 ja_lt_max= ja < ja_max;
8221 }
8222 else if ((! ja_lt_max)
8223 || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
8224 {
8225 double tmp = octave::math::max (0., b.data (jb),
8226 nanflag, realabs);
8227 if (tmp != 0.)
8228 {
8229 r.ridx (jx) = b.ridx (jb);
8230 r.data (jx) = tmp;
8231 jx++;
8232 }
8233 jb++;
8234 jb_lt_max= jb < jb_max;
8235 }
8236 else
8237 {
8238 double tmp = octave::math::max (a.data (ja), b.data (jb),
8239 nanflag, realabs);
8240 if (tmp != 0.)
8241 {
8242 r.data (jx) = tmp;
8243 r.ridx (jx) = a.ridx (ja);
8244 jx++;
8245 }
8246 ja++;
8247 ja_lt_max= ja < ja_max;
8248 jb++;
8249 jb_lt_max= jb < jb_max;
8250 }
8251 }
8252 r.cidx (i+1) = jx;
8253 }
8254
8255 r.maybe_compress ();
8256 }
8257 else
8258 {
8259 if (a_nr == 0 && (b_nr == 0 || b_nr == 1))
8260 {
8261 if (a_nc == 1 || b_nc == 1 || a_nc == b_nc)
8262 r.resize (a_nr, std::max (a_nc, b_nc));
8263 else
8264 octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
8265 }
8266 else if (a_nc == 0 && (b_nc == 0 || b_nc == 1))
8267 {
8268 if (a_nr == 1 || b_nr == 1 || a_nr == b_nr)
8269 r.resize (std::max (a_nr, b_nr), a_nc);
8270 else
8271 octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
8272 }
8273 else if (b_nr == 0 && (a_nr == 0 || a_nr == 1))
8274 {
8275 if (b_nc == 1 || a_nc == 1 || b_nc == a_nc)
8276 r.resize (b_nr, std::max (a_nc, b_nc));
8277 else
8278 octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
8279 }
8280 else if (b_nc == 0 && (a_nc == 0 || a_nc == 1))
8281 {
8282 if (b_nr == 1 || a_nr == 1 || b_nr == a_nr)
8283 r.resize (std::max (a_nr, b_nr), b_nc);
8284 else
8285 octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
8286 }
8287 else
8288 octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
8289 }
8290
8291 return r;
8292}
8293
8296
8299
#define EXPR
#define NAN_EXPR
#define COL_EXPR
#define ROW_EXPR
base_det< double > DET
Definition DET.h:85
#define SPARSE_SSM_BOOL_OPS(S, M)
#define SPARSE_ANY_ALL_HEADER
#define SPARSE_ANY_OP(DIM)
#define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, INIT_VAL, MT_RESULT)
#define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE)
#define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, CONJ_OP)
#define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN, NAN_EXPR, EXPR)
#define SPARSE_XSUM_REDUCTION_OP(RET_TYPE, EL_TYPE)
#define SPARSE_SUMSQ_HEADER(RET_TYPE, EXPR)
#define SPARSE_SSM_CMP_OPS(S, M)
#define SPARSE_SMSM_CMP_OPS(M1, M2)
#define SPARSE_SMS_CMP_OPS(M, S)
#define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE)
#define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN)
#define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, CONJ_OP)
#define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE)
#define SPARSE_SMSM_BOOL_OPS(M1, M2)
#define SPARSE_SMS_BOOL_OPS(M, S)
#define SPARSE_ALL_OP(DIM)
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
N Dimensional Array with copy-on-write semantics.
Definition Array-base.h:130
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array-base.h:547
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition Array-base.h:585
octave_idx_type rows() const
Definition Array-base.h:485
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
octave_idx_type cols() const
Definition Array-base.h:495
const T * data() const
Size of the specified dimension.
Definition Array-base.h:687
T * rwdata()
Size of the specified dimension.
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition CMatrix.h:191
octave_idx_type length() const
Definition DiagArray2.h:92
octave_idx_type cols() const
Definition DiagArray2.h:87
MSparse< T > diag(octave_idx_type k=0) const
Definition MSparse.h:109
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
Definition MSparse.h:84
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
Definition MSparse.h:103
MSparse< T > reshape(const dim_vector &new_dims) const
Definition MSparse.h:100
MSparse< T > squeeze() const
Definition MSparse.h:98
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
Definition MSparse.h:106
bool ishermitian() const
Definition MatrixType.h:121
void mark_as_unsymmetric()
@ Tridiagonal_Hermitian
Definition MatrixType.h:52
@ Banded_Hermitian
Definition MatrixType.h:50
@ Permuted_Diagonal
Definition MatrixType.h:43
int nupper() const
Definition MatrixType.h:100
int nlower() const
Definition MatrixType.h:102
void info() const
octave_idx_type * triangular_perm() const
Definition MatrixType.h:135
MatrixType transpose() const
int type(bool quiet=true)
void mark_as_rectangular()
Definition MatrixType.h:154
bool is_dense() const
Definition MatrixType.h:104
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition dMatrix.h:156
ColumnVector column(octave_idx_type i) const
Definition dMatrix.cc:423
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
Definition CSparse.cc:803
SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
Definition dSparse.cc:168
SparseMatrix squeeze() const
Definition dSparse.cc:7735
SparseMatrix min(int dim=-1, bool nanflag=true, bool realabs=true) const
Definition dSparse.cc:439
bool issymmetric() const
Definition dSparse.cc:128
bool any_element_not_one_or_zero() const
Definition dSparse.cc:7429
SparseMatrix reshape(const dim_vector &new_dims) const
Definition dSparse.cc:7741
SparseBoolMatrix any(int dim=-1) const
Definition dSparse.cc:7553
bool any_element_is_inf_or_nan() const
Definition dSparse.cc:7414
SparseMatrix cumprod(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7560
bool all_integers(double &max_val, double &min_val) const
Definition dSparse.cc:7476
bool any_element_is_nan() const
Definition dSparse.cc:7399
bool too_large_for_float() const
Definition dSparse.cc:7504
bool operator==(const SparseMatrix &a) const
Definition dSparse.cc:98
SparseMatrix max(int dim=-1, bool nanflag=true, bool realabs=true) const
Definition dSparse.cc:183
SparseMatrix transpose() const
Definition dSparse.h:131
DET determinant() const
Definition dSparse.cc:1227
SparseMatrix abs() const
Definition dSparse.cc:7681
SparseMatrix prod(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7588
SparseMatrix inverse() const
Definition dSparse.cc:818
bool operator!=(const SparseMatrix &a) const
Definition dSparse.cc:122
RowVector row(octave_idx_type i) const
Definition dSparse.cc:710
SparseMatrix sum(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7619
SparseMatrix diag(octave_idx_type k=0) const
Definition dSparse.cc:7694
bool all_elements_are_zero() const
Definition dSparse.cc:7444
SparseBoolMatrix all(int dim=-1) const
Definition dSparse.cc:7546
SparseMatrix xsum(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7645
SparseMatrix sumsq(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7653
SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
Definition dSparse.cc:7747
ColumnVector column(octave_idx_type i) const
Definition dSparse.cc:729
SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
Definition dSparse.cc:741
SparseMatrix cumsum(int dim=-1, bool nanflag=false) const
Definition dSparse.cc:7568
SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
Definition dSparse.cc:7753
SparseBoolMatrix operator!() const
Definition dSparse.cc:7510
bool all_elements_are_int_or_inf_or_nan() const
Definition dSparse.cc:7456
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition dSparse.cc:6888
bool any_element_is_negative(bool=false) const
Definition dSparse.cc:7378
Matrix matrix_value() const
Definition dSparse.cc:7700
octave_idx_type cols() const
Definition Sparse.h:351
Array< T > array_value() const
Definition Sparse.cc:2780
octave_idx_type * cidx()
Definition Sparse.h:595
octave_idx_type columns() const
Definition Sparse.h:352
T * data()
Definition Sparse.h:573
void resize(octave_idx_type r, octave_idx_type c)
Definition Sparse.cc:1006
bool test_any(F fcn) const
Definition Sparse.h:662
T element_type
Definition Sparse.h:49
T * xdata()
Definition Sparse.h:575
T & elem(octave_idx_type n)
Definition Sparse.h:455
octave_idx_type * ridx()
Definition Sparse.h:582
Sparse< T, Alloc > maybe_compress(bool remove_zeros=false)
Definition Sparse.h:530
octave_idx_type nnz() const
Actual number of nonzero terms.
Definition Sparse.h:338
octave_idx_type rows() const
Definition Sparse.h:350
void change_capacity(octave_idx_type nz)
Definition Sparse.h:555
dim_vector dims() const
Definition Sparse.h:370
octave_idx_type * xcidx()
Definition Sparse.h:601
octave_idx_type * xridx()
Definition Sparse.h:588
Definition DET.h:38
Vector representing the dimensions (size) of an Array.
Definition dim-vector.h:92
octave_idx_type ndims() const
Number of dimensions.
Definition dim-vector.h:263
int first_non_singleton(int def=0) const
Definition dim-vector.h:481
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
Definition dSparse.cc:7785
SparseMatrix real(const SparseComplexMatrix &a)
Definition dSparse.cc:761
#define EMPTY_RETURN_CHECK(T)
Definition dSparse.cc:7844
SparseMatrix min(double d, const SparseMatrix &m)
Definition dSparse.cc:7849
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
Definition dSparse.cc:7773
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
Definition dSparse.cc:7811
SparseMatrix max(double d, const SparseMatrix &m)
Definition dSparse.cc:8072
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
Definition dSparse.cc:7761
SparseMatrix imag(const SparseComplexMatrix &a)
Definition dSparse.cc:782
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)
Definition dSparse.cc:7706
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
Definition dSparse.cc:7805
std::istream & operator>>(std::istream &is, SparseMatrix &a)
Definition dSparse.cc:7727
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
F77_RET_T const F77_INT F77_CMPLX * A
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
std::complex< double > Complex
Definition oct-cmplx.h:33
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
#define CHOLMOD_NAME(name)
Definition oct-sparse.h:140
#define UMFPACK_DNAME(name)
Definition oct-sparse.h:169
const octave_base_value const Array< octave_idx_type > & ra_idx
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
template SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
template Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
template ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
void SparseCholError(int status, char *file, int line, char *message)
int SparseCholPrint(const char *fmt,...)