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