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