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