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