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