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