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