GNU Octave  8.1.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-2023 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 {
631  SparseMatrix retval;
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 {
681  SparseMatrix retval;
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 {
1223  SparseMatrix retval;
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 {
1312  ComplexMatrix retval;
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 {
1371  SparseComplexMatrix retval;
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 {
1690  SparseMatrix retval;
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 {
1971  ComplexMatrix retval;
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 {
2203  SparseComplexMatrix retval;
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 {
2740  SparseMatrix retval;
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 {
3040  ComplexMatrix retval;
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 {
3294  SparseComplexMatrix retval;
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 {
3762  SparseMatrix retval;
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 {
3916  ComplexMatrix retval;
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 {
4086  SparseComplexMatrix retval;
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 {
4545  SparseMatrix retval;
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 {
4885  ComplexMatrix retval;
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 {
5209  SparseComplexMatrix retval;
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.data ());
5800 
5801  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
5802  CHOLMOD_NAME(factorize) (A, L, cm);
5803  if (calc_cond)
5804  rcond = CHOLMOD_NAME(rcond)(L, cm);
5805  else
5806  rcond = 1.0;
5807 
5808  if (rcond == 0.0)
5809  {
5810  // Either its indefinite or singular. Try UMFPACK
5811  mattype.mark_as_unsymmetric ();
5812  typ = MatrixType::Full;
5813  }
5814  else
5815  {
5816  volatile double rcond_plus_one = rcond + 1.0;
5817 
5818  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
5819  {
5820  err = -2;
5821 
5822  if (sing_handler)
5823  {
5824  sing_handler (rcond);
5825  mattype.mark_as_rectangular ();
5826  }
5827  else
5829 
5830  return retval;
5831  }
5832 
5833  cholmod_dense *X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5834 
5835  retval.resize (b.rows (), b.cols ());
5836  for (octave_idx_type j = 0; j < b.cols (); j++)
5837  {
5838  octave_idx_type jr = j * b.rows ();
5839  for (octave_idx_type i = 0; i < b.rows (); i++)
5840  retval.xelem (i, j) = static_cast<double *>(X->x)[jr + i];
5841  }
5842 
5843  CHOLMOD_NAME(free_dense) (&X, cm);
5844  CHOLMOD_NAME(free_factor) (&L, cm);
5845  CHOLMOD_NAME(finish) (cm);
5846  static char blank_name[] = " ";
5847  CHOLMOD_NAME(print_common) (blank_name, cm);
5848  }
5849 #else
5850  (*current_liboctave_warning_with_id_handler)
5851  ("Octave:missing-dependency",
5852  "support for CHOLMOD was unavailable or disabled "
5853  "when liboctave was built");
5854 
5855  mattype.mark_as_unsymmetric ();
5856  typ = MatrixType::Full;
5857 #endif
5858  }
5859 
5860  if (typ == MatrixType::Full)
5861  {
5862 #if defined (HAVE_UMFPACK)
5863  Matrix Control, Info;
5864  void *Numeric
5865  = factorize (err, rcond, Control, Info, sing_handler, calc_cond);
5866 
5867  if (err == 0)
5868  {
5869  // one iterative refinement instead of the default two in UMFPACK
5870  Control (UMFPACK_IRSTEP) = 1;
5871  const double *Bx = b.data ();
5872  retval.resize (b.rows (), b.cols ());
5873  double *result = retval.fortran_vec ();
5874  octave_idx_type b_nr = b.rows ();
5875  octave_idx_type b_nc = b.cols ();
5876  int status = 0;
5877  double *control = Control.fortran_vec ();
5878  double *info = Info.fortran_vec ();
5879  const octave_idx_type *Ap = cidx ();
5880  const octave_idx_type *Ai = ridx ();
5881  const double *Ax = data ();
5882 
5883  for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
5884  {
5885  status = UMFPACK_DNAME (solve) (UMFPACK_A,
5888  Ax, &result[iidx], &Bx[iidx],
5889  Numeric, control, info);
5890  if (status < 0)
5891  {
5892  UMFPACK_DNAME (report_status) (control, status);
5893 
5894  // FIXME: Should this be a warning?
5895  (*current_liboctave_error_handler)
5896  ("SparseMatrix::solve solve failed");
5897 
5898  err = -1;
5899  break;
5900  }
5901  }
5902 
5903  UMFPACK_DNAME (report_info) (control, info);
5904 
5905  UMFPACK_DNAME (free_numeric) (&Numeric);
5906  }
5907  else
5908  mattype.mark_as_rectangular ();
5909 
5910 #else
5911  octave_unused_parameter (rcond);
5912  octave_unused_parameter (sing_handler);
5913  octave_unused_parameter (calc_cond);
5914 
5915  (*current_liboctave_error_handler)
5916  ("support for UMFPACK was unavailable or disabled "
5917  "when liboctave was built");
5918 #endif
5919  }
5920  else if (typ != MatrixType::Hermitian)
5921  (*current_liboctave_error_handler) ("incorrect matrix type");
5922  }
5923 
5924  return retval;
5925 }
5926 
5929  octave_idx_type& err, double& rcond,
5930  solve_singularity_handler sing_handler,
5931  bool calc_cond) const
5932 {
5933  SparseMatrix retval;
5934 
5935  octave_idx_type nr = rows ();
5936  octave_idx_type nc = cols ();
5937  err = 0;
5938 
5939  if (nr != nc || nr != b.rows ())
5940  (*current_liboctave_error_handler)
5941  ("matrix dimension mismatch solution of linear equations");
5942 
5943  if (nr == 0 || b.cols () == 0)
5944  retval = SparseMatrix (nc, b.cols ());
5945  else
5946  {
5947  // Print spparms("spumoni") info if requested
5948  volatile int typ = mattype.type ();
5949  mattype.info ();
5950 
5951  if (typ == MatrixType::Hermitian)
5952  {
5953 #if defined (HAVE_CHOLMOD)
5954  cholmod_common Common;
5955  cholmod_common *cm = &Common;
5956 
5957  // Setup initial parameters
5958  CHOLMOD_NAME(start) (cm);
5959  cm->prefer_zomplex = false;
5960 
5961  double spu = octave::sparse_params::get_key ("spumoni");
5962  if (spu == 0.)
5963  {
5964  cm->print = -1;
5965  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
5966  }
5967  else
5968  {
5969  cm->print = static_cast<int> (spu) + 2;
5970  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
5971  }
5972 
5973  cm->error_handler = &SparseCholError;
5974  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5975  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5976 
5977  cm->final_ll = true;
5978 
5979  cholmod_sparse Astore;
5980  cholmod_sparse *A = &Astore;
5981  A->nrow = nr;
5982  A->ncol = nc;
5983 
5984  A->p = cidx ();
5985  A->i = ridx ();
5986  A->nzmax = nnz ();
5987  A->packed = true;
5988  A->sorted = true;
5989  A->nz = nullptr;
5990 #if defined (OCTAVE_ENABLE_64)
5991  A->itype = CHOLMOD_LONG;
5992 #else
5993  A->itype = CHOLMOD_INT;
5994 #endif
5995  A->dtype = CHOLMOD_DOUBLE;
5996  A->stype = 1;
5997  A->xtype = CHOLMOD_REAL;
5998 
5999  A->x = data ();
6000 
6001  cholmod_sparse Bstore;
6002  cholmod_sparse *B = &Bstore;
6003  B->nrow = b.rows ();
6004  B->ncol = b.cols ();
6005  B->p = b.cidx ();
6006  B->i = b.ridx ();
6007  B->nzmax = b.nnz ();
6008  B->packed = true;
6009  B->sorted = true;
6010  B->nz = nullptr;
6011 #if defined (OCTAVE_ENABLE_64)
6012  B->itype = CHOLMOD_LONG;
6013 #else
6014  B->itype = CHOLMOD_INT;
6015 #endif
6016  B->dtype = CHOLMOD_DOUBLE;
6017  B->stype = 0;
6018  B->xtype = CHOLMOD_REAL;
6019 
6020  B->x = b.data ();
6021 
6022  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6023  CHOLMOD_NAME(factorize) (A, L, cm);
6024  if (calc_cond)
6025  rcond = CHOLMOD_NAME(rcond)(L, cm);
6026  else
6027  rcond = 1.;
6028 
6029  if (rcond == 0.0)
6030  {
6031  // Either its indefinite or singular. Try UMFPACK
6032  mattype.mark_as_unsymmetric ();
6033  typ = MatrixType::Full;
6034  }
6035  else
6036  {
6037  volatile double rcond_plus_one = rcond + 1.0;
6038 
6039  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
6040  {
6041  err = -2;
6042 
6043  if (sing_handler)
6044  {
6045  sing_handler (rcond);
6046  mattype.mark_as_rectangular ();
6047  }
6048  else
6050 
6051  return retval;
6052  }
6053 
6054  cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6055 
6056  retval = SparseMatrix (static_cast<octave_idx_type> (X->nrow),
6057  static_cast<octave_idx_type> (X->ncol),
6058  static_cast<octave_idx_type> (X->nzmax));
6059  for (octave_idx_type j = 0;
6060  j <= static_cast<octave_idx_type> (X->ncol); j++)
6061  retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6062  for (octave_idx_type j = 0;
6063  j < static_cast<octave_idx_type> (X->nzmax); j++)
6064  {
6065  retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6066  retval.xdata (j) = static_cast<double *>(X->x)[j];
6067  }
6068 
6069  CHOLMOD_NAME(free_sparse) (&X, cm);
6070  CHOLMOD_NAME(free_factor) (&L, cm);
6071  CHOLMOD_NAME(finish) (cm);
6072  static char blank_name[] = " ";
6073  CHOLMOD_NAME(print_common) (blank_name, cm);
6074  }
6075 #else
6076  (*current_liboctave_warning_with_id_handler)
6077  ("Octave:missing-dependency",
6078  "support for CHOLMOD was unavailable or disabled "
6079  "when liboctave was built");
6080 
6081  mattype.mark_as_unsymmetric ();
6082  typ = MatrixType::Full;
6083 #endif
6084  }
6085 
6086  if (typ == MatrixType::Full)
6087  {
6088 #if defined (HAVE_UMFPACK)
6089  Matrix Control, Info;
6090  void *Numeric = factorize (err, rcond, Control, Info,
6091  sing_handler, calc_cond);
6092 
6093  if (err == 0)
6094  {
6095  // one iterative refinement instead of the default two in UMFPACK
6096  Control (UMFPACK_IRSTEP) = 1;
6097  octave_idx_type b_nr = b.rows ();
6098  octave_idx_type b_nc = b.cols ();
6099  int status = 0;
6100  double *control = Control.fortran_vec ();
6101  double *info = Info.fortran_vec ();
6102  const octave_idx_type *Ap = cidx ();
6103  const octave_idx_type *Ai = ridx ();
6104  const double *Ax = data ();
6105 
6106  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6107  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6108 
6109  // Take a first guess that the number of nonzero terms
6110  // will be as many as in b
6111  octave_idx_type x_nz = b.nnz ();
6112  octave_idx_type ii = 0;
6113  retval = SparseMatrix (b_nr, b_nc, x_nz);
6114 
6115  retval.xcidx (0) = 0;
6116  for (octave_idx_type j = 0; j < b_nc; j++)
6117  {
6118 
6119  for (octave_idx_type i = 0; i < b_nr; i++)
6120  Bx[i] = b.elem (i, j);
6121 
6122  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6125  Ax, Xx, Bx, Numeric,
6126  control, info);
6127  if (status < 0)
6128  {
6129  UMFPACK_DNAME (report_status) (control, status);
6130 
6131  // FIXME: Should this be a warning?
6132  (*current_liboctave_error_handler)
6133  ("SparseMatrix::solve solve failed");
6134 
6135  err = -1;
6136  break;
6137  }
6138 
6139  for (octave_idx_type i = 0; i < b_nr; i++)
6140  {
6141  double tmp = Xx[i];
6142  if (tmp != 0.0)
6143  {
6144  if (ii == x_nz)
6145  {
6146  // Resize the sparse matrix
6147  octave_idx_type sz;
6148  sz = (static_cast<double> (b_nc) - j) / b_nc
6149  * x_nz;
6150  sz = x_nz + (sz > 100 ? sz : 100);
6151  retval.change_capacity (sz);
6152  x_nz = sz;
6153  }
6154  retval.xdata (ii) = tmp;
6155  retval.xridx (ii++) = i;
6156  }
6157  }
6158  retval.xcidx (j+1) = ii;
6159  }
6160 
6161  retval.maybe_compress ();
6162 
6163  UMFPACK_DNAME (report_info) (control, info);
6164 
6165  UMFPACK_DNAME (free_numeric) (&Numeric);
6166  }
6167  else
6168  mattype.mark_as_rectangular ();
6169 
6170 #else
6171  octave_unused_parameter (rcond);
6172  octave_unused_parameter (sing_handler);
6173  octave_unused_parameter (calc_cond);
6174 
6175  (*current_liboctave_error_handler)
6176  ("support for UMFPACK was unavailable or disabled "
6177  "when liboctave was built");
6178 #endif
6179  }
6180  else if (typ != MatrixType::Hermitian)
6181  (*current_liboctave_error_handler) ("incorrect matrix type");
6182  }
6183 
6184  return retval;
6185 }
6186 
6189  octave_idx_type& err, double& rcond,
6190  solve_singularity_handler sing_handler,
6191  bool calc_cond) const
6192 {
6193  ComplexMatrix retval;
6194 
6195  octave_idx_type nr = rows ();
6196  octave_idx_type nc = cols ();
6197  err = 0;
6198 
6199  if (nr != nc || nr != b.rows ())
6200  (*current_liboctave_error_handler)
6201  ("matrix dimension mismatch solution of linear equations");
6202 
6203  if (nr == 0 || b.cols () == 0)
6204  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
6205  else
6206  {
6207  // Print spparms("spumoni") info if requested
6208  volatile int typ = mattype.type ();
6209  mattype.info ();
6210 
6211  if (typ == MatrixType::Hermitian)
6212  {
6213 #if defined (HAVE_CHOLMOD)
6214  cholmod_common Common;
6215  cholmod_common *cm = &Common;
6216 
6217  // Setup initial parameters
6218  CHOLMOD_NAME(start) (cm);
6219  cm->prefer_zomplex = false;
6220 
6221  double spu = octave::sparse_params::get_key ("spumoni");
6222  if (spu == 0.)
6223  {
6224  cm->print = -1;
6225  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6226  }
6227  else
6228  {
6229  cm->print = static_cast<int> (spu) + 2;
6230  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6231  }
6232 
6233  cm->error_handler = &SparseCholError;
6234  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6235  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6236 
6237  cm->final_ll = true;
6238 
6239  cholmod_sparse Astore;
6240  cholmod_sparse *A = &Astore;
6241  A->nrow = nr;
6242  A->ncol = nc;
6243 
6244  A->p = cidx ();
6245  A->i = ridx ();
6246  A->nzmax = nnz ();
6247  A->packed = true;
6248  A->sorted = true;
6249  A->nz = nullptr;
6250 #if defined (OCTAVE_ENABLE_64)
6251  A->itype = CHOLMOD_LONG;
6252 #else
6253  A->itype = CHOLMOD_INT;
6254 #endif
6255  A->dtype = CHOLMOD_DOUBLE;
6256  A->stype = 1;
6257  A->xtype = CHOLMOD_REAL;
6258 
6259  A->x = data ();
6260 
6261  cholmod_dense Bstore;
6262  cholmod_dense *B = &Bstore;
6263  B->nrow = b.rows ();
6264  B->ncol = b.cols ();
6265  B->d = B->nrow;
6266  B->nzmax = B->nrow * B->ncol;
6267  B->dtype = CHOLMOD_DOUBLE;
6268  B->xtype = CHOLMOD_COMPLEX;
6269 
6270  B->x = const_cast<Complex *> (b.data ());
6271 
6272  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6273  CHOLMOD_NAME(factorize) (A, L, cm);
6274  if (calc_cond)
6275  rcond = CHOLMOD_NAME(rcond)(L, cm);
6276  else
6277  rcond = 1.0;
6278 
6279  if (rcond == 0.0)
6280  {
6281  // Either its indefinite or singular. Try UMFPACK
6282  mattype.mark_as_unsymmetric ();
6283  typ = MatrixType::Full;
6284  }
6285  else
6286  {
6287  volatile double rcond_plus_one = rcond + 1.0;
6288 
6289  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
6290  {
6291  err = -2;
6292 
6293  if (sing_handler)
6294  {
6295  sing_handler (rcond);
6296  mattype.mark_as_rectangular ();
6297  }
6298  else
6300 
6301  return retval;
6302  }
6303 
6304  cholmod_dense *X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
6305 
6306  retval.resize (b.rows (), b.cols ());
6307  for (octave_idx_type j = 0; j < b.cols (); j++)
6308  {
6309  octave_idx_type jr = j * b.rows ();
6310  for (octave_idx_type i = 0; i < b.rows (); i++)
6311  retval.xelem (i, j) = static_cast<Complex *>(X->x)[jr + i];
6312  }
6313 
6314  CHOLMOD_NAME(free_dense) (&X, cm);
6315  CHOLMOD_NAME(free_factor) (&L, cm);
6316  CHOLMOD_NAME(finish) (cm);
6317  static char blank_name[] = " ";
6318  CHOLMOD_NAME(print_common) (blank_name, cm);
6319  }
6320 #else
6321  (*current_liboctave_warning_with_id_handler)
6322  ("Octave:missing-dependency",
6323  "support for CHOLMOD was unavailable or disabled "
6324  "when liboctave was built");
6325 
6326  mattype.mark_as_unsymmetric ();
6327  typ = MatrixType::Full;
6328 #endif
6329  }
6330 
6331  if (typ == MatrixType::Full)
6332  {
6333 #if defined (HAVE_UMFPACK)
6334  Matrix Control, Info;
6335  void *Numeric = factorize (err, rcond, Control, Info,
6336  sing_handler, calc_cond);
6337 
6338  if (err == 0)
6339  {
6340  // one iterative refinement instead of the default two in UMFPACK
6341  Control (UMFPACK_IRSTEP) = 1;
6342  octave_idx_type b_nr = b.rows ();
6343  octave_idx_type b_nc = b.cols ();
6344  int status = 0;
6345  double *control = Control.fortran_vec ();
6346  double *info = Info.fortran_vec ();
6347  const octave_idx_type *Ap = cidx ();
6348  const octave_idx_type *Ai = ridx ();
6349  const double *Ax = data ();
6350 
6351  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6352  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6353 
6354  retval.resize (b_nr, b_nc);
6355 
6356  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6357  OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6358 
6359  for (octave_idx_type j = 0; j < b_nc; j++)
6360  {
6361  for (octave_idx_type i = 0; i < b_nr; i++)
6362  {
6363  Complex c = b(i, j);
6364  Bx[i] = c.real ();
6365  Bz[i] = c.imag ();
6366  }
6367 
6368  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6371  Ax, Xx, Bx, Numeric,
6372  control, info);
6373  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6376  Ax, Xz, Bz,
6377  Numeric, control, info);
6378 
6379  if (status < 0 || status2 < 0)
6380  {
6381  UMFPACK_DNAME (report_status) (control, status);
6382 
6383  // FIXME: Should this be a warning?
6384  (*current_liboctave_error_handler)
6385  ("SparseMatrix::solve solve failed");
6386 
6387  err = -1;
6388  break;
6389  }
6390 
6391  for (octave_idx_type i = 0; i < b_nr; i++)
6392  retval(i, j) = Complex (Xx[i], Xz[i]);
6393  }
6394 
6395  UMFPACK_DNAME (report_info) (control, info);
6396 
6397  UMFPACK_DNAME (free_numeric) (&Numeric);
6398  }
6399  else
6400  mattype.mark_as_rectangular ();
6401 
6402 #else
6403  octave_unused_parameter (rcond);
6404  octave_unused_parameter (sing_handler);
6405  octave_unused_parameter (calc_cond);
6406 
6407  (*current_liboctave_error_handler)
6408  ("support for UMFPACK was unavailable or disabled "
6409  "when liboctave was built");
6410 #endif
6411  }
6412  else if (typ != MatrixType::Hermitian)
6413  (*current_liboctave_error_handler) ("incorrect matrix type");
6414  }
6415 
6416  return retval;
6417 }
6418 
6421  octave_idx_type& err, double& rcond,
6422  solve_singularity_handler sing_handler,
6423  bool calc_cond) const
6424 {
6425  SparseComplexMatrix retval;
6426 
6427  octave_idx_type nr = rows ();
6428  octave_idx_type nc = cols ();
6429  err = 0;
6430 
6431  if (nr != nc || nr != b.rows ())
6432  (*current_liboctave_error_handler)
6433  ("matrix dimension mismatch solution of linear equations");
6434 
6435  if (nr == 0 || b.cols () == 0)
6436  retval = SparseComplexMatrix (nc, b.cols ());
6437  else
6438  {
6439  // Print spparms("spumoni") info if requested
6440  volatile int typ = mattype.type ();
6441  mattype.info ();
6442 
6443  if (typ == MatrixType::Hermitian)
6444  {
6445 #if defined (HAVE_CHOLMOD)
6446  cholmod_common Common;
6447  cholmod_common *cm = &Common;
6448 
6449  // Setup initial parameters
6450  CHOLMOD_NAME(start) (cm);
6451  cm->prefer_zomplex = false;
6452 
6453  double spu = octave::sparse_params::get_key ("spumoni");
6454  if (spu == 0.)
6455  {
6456  cm->print = -1;
6457  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6458  }
6459  else
6460  {
6461  cm->print = static_cast<int> (spu) + 2;
6462  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6463  }
6464 
6465  cm->error_handler = &SparseCholError;
6466  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6467  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6468 
6469  cm->final_ll = true;
6470 
6471  cholmod_sparse Astore;
6472  cholmod_sparse *A = &Astore;
6473  A->nrow = nr;
6474  A->ncol = nc;
6475 
6476  A->p = cidx ();
6477  A->i = ridx ();
6478  A->nzmax = nnz ();
6479  A->packed = true;
6480  A->sorted = true;
6481  A->nz = nullptr;
6482 #if defined (OCTAVE_ENABLE_64)
6483  A->itype = CHOLMOD_LONG;
6484 #else
6485  A->itype = CHOLMOD_INT;
6486 #endif
6487  A->dtype = CHOLMOD_DOUBLE;
6488  A->stype = 1;
6489  A->xtype = CHOLMOD_REAL;
6490 
6491  A->x = data ();
6492 
6493  cholmod_sparse Bstore;
6494  cholmod_sparse *B = &Bstore;
6495  B->nrow = b.rows ();
6496  B->ncol = b.cols ();
6497  B->p = b.cidx ();
6498  B->i = b.ridx ();
6499  B->nzmax = b.nnz ();
6500  B->packed = true;
6501  B->sorted = true;
6502  B->nz = nullptr;
6503 #if defined (OCTAVE_ENABLE_64)
6504  B->itype = CHOLMOD_LONG;
6505 #else
6506  B->itype = CHOLMOD_INT;
6507 #endif
6508  B->dtype = CHOLMOD_DOUBLE;
6509  B->stype = 0;
6510  B->xtype = CHOLMOD_COMPLEX;
6511 
6512  B->x = b.data ();
6513 
6514  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6515  CHOLMOD_NAME(factorize) (A, L, cm);
6516  if (calc_cond)
6517  rcond = CHOLMOD_NAME(rcond)(L, cm);
6518  else
6519  rcond = 1.0;
6520 
6521  if (rcond == 0.0)
6522  {
6523  // Either its indefinite or singular. Try UMFPACK
6524  mattype.mark_as_unsymmetric ();
6525  typ = MatrixType::Full;
6526  }
6527  else
6528  {
6529  volatile double rcond_plus_one = rcond + 1.0;
6530 
6531  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
6532  {
6533  err = -2;
6534 
6535  if (sing_handler)
6536  {
6537  sing_handler (rcond);
6538  mattype.mark_as_rectangular ();
6539  }
6540  else
6542 
6543  return retval;
6544  }
6545 
6546  cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6547 
6548  retval = SparseComplexMatrix
6549  (static_cast<octave_idx_type> (X->nrow),
6550  static_cast<octave_idx_type> (X->ncol),
6551  static_cast<octave_idx_type> (X->nzmax));
6552  for (octave_idx_type j = 0;
6553  j <= static_cast<octave_idx_type> (X->ncol); j++)
6554  retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6555  for (octave_idx_type j = 0;
6556  j < static_cast<octave_idx_type> (X->nzmax); j++)
6557  {
6558  retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6559  retval.xdata (j) = static_cast<Complex *>(X->x)[j];
6560  }
6561 
6562  CHOLMOD_NAME(free_sparse) (&X, cm);
6563  CHOLMOD_NAME(free_factor) (&L, cm);
6564  CHOLMOD_NAME(finish) (cm);
6565  static char blank_name[] = " ";
6566  CHOLMOD_NAME(print_common) (blank_name, cm);
6567  }
6568 #else
6569  (*current_liboctave_warning_with_id_handler)
6570  ("Octave:missing-dependency",
6571  "support for CHOLMOD was unavailable or disabled "
6572  "when liboctave was built");
6573 
6574  mattype.mark_as_unsymmetric ();
6575  typ = MatrixType::Full;
6576 #endif
6577  }
6578 
6579  if (typ == MatrixType::Full)
6580  {
6581 #if defined (HAVE_UMFPACK)
6582  Matrix Control, Info;
6583  void *Numeric = factorize (err, rcond, Control, Info,
6584  sing_handler, calc_cond);
6585 
6586  if (err == 0)
6587  {
6588  // one iterative refinement instead of the default two in UMFPACK
6589  Control (UMFPACK_IRSTEP) = 1;
6590  octave_idx_type b_nr = b.rows ();
6591  octave_idx_type b_nc = b.cols ();
6592  int status = 0;
6593  double *control = Control.fortran_vec ();
6594  double *info = Info.fortran_vec ();
6595  const octave_idx_type *Ap = cidx ();
6596  const octave_idx_type *Ai = ridx ();
6597  const double *Ax = data ();
6598 
6599  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6600  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6601 
6602  // Take a first guess that the number of nonzero terms
6603  // will be as many as in b
6604  octave_idx_type x_nz = b.nnz ();
6605  octave_idx_type ii = 0;
6606  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6607 
6608  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6609  OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6610 
6611  retval.xcidx (0) = 0;
6612  for (octave_idx_type j = 0; j < b_nc; j++)
6613  {
6614  for (octave_idx_type i = 0; i < b_nr; i++)
6615  {
6616  Complex c = b(i, j);
6617  Bx[i] = c.real ();
6618  Bz[i] = c.imag ();
6619  }
6620 
6621  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6624  Ax, Xx, Bx, Numeric,
6625  control, info);
6626  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6629  Ax, Xz, Bz,
6630  Numeric, control, info);
6631 
6632  if (status < 0 || status2 < 0)
6633  {
6634  UMFPACK_DNAME (report_status) (control, status);
6635 
6636  // FIXME: Should this be a warning?
6637  (*current_liboctave_error_handler)
6638  ("SparseMatrix::solve solve failed");
6639 
6640  err = -1;
6641  break;
6642  }
6643 
6644  for (octave_idx_type i = 0; i < b_nr; i++)
6645  {
6646  Complex tmp = Complex (Xx[i], Xz[i]);
6647  if (tmp != 0.0)
6648  {
6649  if (ii == x_nz)
6650  {
6651  // Resize the sparse matrix
6652  octave_idx_type sz;
6653  sz = (static_cast<double> (b_nc) - j) / b_nc
6654  * x_nz;
6655  sz = x_nz + (sz > 100 ? sz : 100);
6656  retval.change_capacity (sz);
6657  x_nz = sz;
6658  }
6659  retval.xdata (ii) = tmp;
6660  retval.xridx (ii++) = i;
6661  }
6662  }
6663  retval.xcidx (j+1) = ii;
6664  }
6665 
6666  retval.maybe_compress ();
6667 
6668  UMFPACK_DNAME (report_info) (control, info);
6669 
6670  UMFPACK_DNAME (free_numeric) (&Numeric);
6671  }
6672  else
6673  mattype.mark_as_rectangular ();
6674 #else
6675  octave_unused_parameter (rcond);
6676  octave_unused_parameter (sing_handler);
6677  octave_unused_parameter (calc_cond);
6678 
6679  (*current_liboctave_error_handler)
6680  ("support for UMFPACK was unavailable or disabled "
6681  "when liboctave was built");
6682 #endif
6683  }
6684  else if (typ != MatrixType::Hermitian)
6685  (*current_liboctave_error_handler) ("incorrect matrix type");
6686  }
6687 
6688  return retval;
6689 }
6690 
6691 Matrix
6692 SparseMatrix::solve (MatrixType& mattype, const Matrix& b) const
6693 {
6694  octave_idx_type info;
6695  double rcond;
6696  return solve (mattype, b, info, rcond, nullptr);
6697 }
6698 
6699 Matrix
6701  octave_idx_type& info) const
6702 {
6703  double rcond;
6704  return solve (mattype, b, info, rcond, nullptr);
6705 }
6706 
6707 Matrix
6709  octave_idx_type& info, double& rcond) const
6710 {
6711  return solve (mattype, b, info, rcond, nullptr);
6712 }
6713 
6714 Matrix
6716  double& rcond, solve_singularity_handler sing_handler,
6717  bool singular_fallback) const
6718 {
6719  Matrix retval;
6720  int typ = mattype.type (false);
6721 
6722  if (typ == MatrixType::Unknown)
6723  typ = mattype.type (*this);
6724 
6725  // Only calculate the condition number for CHOLMOD/UMFPACK
6727  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6728  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6729  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6730  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6731  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6732  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6733  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6734  else if (typ == MatrixType::Tridiagonal
6736  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6737  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6738  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6739  else if (typ != MatrixType::Rectangular)
6740  (*current_liboctave_error_handler) ("unknown matrix type");
6741 
6742  // Rectangular or one of the above solvers flags a singular matrix
6743  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6744  {
6745  rcond = 1.;
6746 #if defined (USE_QRSOLVE)
6747  retval = qrsolve (*this, b, err);
6748 #else
6749  retval = dmsolve<Matrix, SparseMatrix, Matrix> (*this, b, err);
6750 #endif
6751  }
6752 
6753  return retval;
6754 }
6755 
6757 SparseMatrix::solve (MatrixType& mattype, const SparseMatrix& b) const
6758 {
6759  octave_idx_type info;
6760  double rcond;
6761  return solve (mattype, b, info, rcond, nullptr);
6762 }
6763 
6766  octave_idx_type& info) const
6767 {
6768  double rcond;
6769  return solve (mattype, b, info, rcond, nullptr);
6770 }
6771 
6774  octave_idx_type& info, double& rcond) const
6775 {
6776  return solve (mattype, b, info, rcond, nullptr);
6777 }
6778 
6781  octave_idx_type& err, double& rcond,
6782  solve_singularity_handler sing_handler,
6783  bool singular_fallback) const
6784 {
6785  SparseMatrix retval;
6786  int typ = mattype.type (false);
6787 
6788  if (typ == MatrixType::Unknown)
6789  typ = mattype.type (*this);
6790 
6792  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6793  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6794  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6795  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6796  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6797  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6798  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6799  else if (typ == MatrixType::Tridiagonal
6801  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6802  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6803  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6804  else if (typ != MatrixType::Rectangular)
6805  (*current_liboctave_error_handler) ("unknown matrix type");
6806 
6807  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6808  {
6809  rcond = 1.;
6810 #if defined (USE_QRSOLVE)
6811  retval = qrsolve (*this, b, err);
6812 #else
6814  (*this, b, err);
6815 #endif
6816  }
6817 
6818  return retval;
6819 }
6820 
6823 {
6824  octave_idx_type info;
6825  double rcond;
6826  return solve (mattype, b, info, rcond, nullptr);
6827 }
6828 
6831  octave_idx_type& info) const
6832 {
6833  double rcond;
6834  return solve (mattype, b, info, rcond, nullptr);
6835 }
6836 
6839  octave_idx_type& info, double& rcond) const
6840 {
6841  return solve (mattype, b, info, rcond, nullptr);
6842 }
6843 
6846  octave_idx_type& err, double& rcond,
6847  solve_singularity_handler sing_handler,
6848  bool singular_fallback) const
6849 {
6850  ComplexMatrix retval;
6851  int typ = mattype.type (false);
6852 
6853  if (typ == MatrixType::Unknown)
6854  typ = mattype.type (*this);
6855 
6857  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6858  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6859  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6860  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6861  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6862  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6863  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6864  else if (typ == MatrixType::Tridiagonal
6866  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6867  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6868  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6869  else if (typ != MatrixType::Rectangular)
6870  (*current_liboctave_error_handler) ("unknown matrix type");
6871 
6872  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6873  {
6874  rcond = 1.;
6875 #if defined (USE_QRSOLVE)
6876  retval = qrsolve (*this, b, err);
6877 #else
6879  (*this, b, err);
6880 #endif
6881  }
6882 
6883  return retval;
6884 }
6885 
6888 {
6889  octave_idx_type info;
6890  double rcond;
6891  return solve (mattype, b, info, rcond, nullptr);
6892 }
6893 
6896  octave_idx_type& info) const
6897 {
6898  double rcond;
6899  return solve (mattype, b, info, rcond, nullptr);
6900 }
6901 
6904  octave_idx_type& info, double& rcond) const
6905 {
6906  return solve (mattype, b, info, rcond, nullptr);
6907 }
6908 
6911  octave_idx_type& err, double& rcond,
6912  solve_singularity_handler sing_handler,
6913  bool singular_fallback) const
6914 {
6915  SparseComplexMatrix retval;
6916  int typ = mattype.type (false);
6917 
6918  if (typ == MatrixType::Unknown)
6919  typ = mattype.type (*this);
6920 
6922  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6923  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6924  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6925  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6926  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6927  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6928  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6929  else if (typ == MatrixType::Tridiagonal
6931  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6932  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6933  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6934  else if (typ != MatrixType::Rectangular)
6935  (*current_liboctave_error_handler) ("unknown matrix type");
6936 
6937  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6938  {
6939  rcond = 1.;
6940 #if defined (USE_QRSOLVE)
6941  retval = qrsolve (*this, b, err);
6942 #else
6944  (*this, b, err);
6945 #endif
6946  }
6947 
6948  return retval;
6949 }
6950 
6952 SparseMatrix::solve (MatrixType& mattype, const ColumnVector& b) const
6953 {
6954  octave_idx_type info; double rcond;
6955  return solve (mattype, b, info, rcond);
6956 }
6957 
6960  octave_idx_type& info) const
6961 {
6962  double rcond;
6963  return solve (mattype, b, info, rcond);
6964 }
6965 
6968  octave_idx_type& info, double& rcond) const
6969 {
6970  return solve (mattype, b, info, rcond, nullptr);
6971 }
6972 
6975  octave_idx_type& info, double& rcond,
6976  solve_singularity_handler sing_handler) const
6977 {
6978  Matrix tmp (b);
6979  return solve (mattype, tmp, info, rcond,
6980  sing_handler).column (static_cast<octave_idx_type> (0));
6981 }
6982 
6985 {
6986  octave_idx_type info;
6987  double rcond;
6988  return solve (mattype, b, info, rcond, nullptr);
6989 }
6990 
6993  octave_idx_type& info) const
6994 {
6995  double rcond;
6996  return solve (mattype, b, info, rcond, nullptr);
6997 }
6998 
7001  octave_idx_type& info,
7002  double& rcond) const
7003 {
7004  return solve (mattype, b, info, rcond, nullptr);
7005 }
7006 
7009  octave_idx_type& info, double& rcond,
7010  solve_singularity_handler sing_handler) const
7011 {
7012  ComplexMatrix tmp (b);
7013  return solve (mattype, tmp, info, rcond,
7014  sing_handler).column (static_cast<octave_idx_type> (0));
7015 }
7016 
7017 Matrix
7018 SparseMatrix::solve (const Matrix& b) const
7019 {
7020  octave_idx_type info;
7021  double rcond;
7022  return solve (b, info, rcond, nullptr);
7023 }
7024 
7025 Matrix
7027 {
7028  double rcond;
7029  return solve (b, info, rcond, nullptr);
7030 }
7031 
7032 Matrix
7034  double& rcond) const
7035 {
7036  return solve (b, info, rcond, nullptr);
7037 }
7038 
7039 Matrix
7040 SparseMatrix::solve (const Matrix& b, octave_idx_type& err, double& rcond,
7041  solve_singularity_handler sing_handler) const
7042 {
7043  MatrixType mattype (*this);
7044  return solve (mattype, b, err, rcond, sing_handler);
7045 }
7046 
7049 {
7050  octave_idx_type info;
7051  double rcond;
7052  return solve (b, info, rcond, nullptr);
7053 }
7054 
7057  octave_idx_type& info) const
7058 {
7059  double rcond;
7060  return solve (b, info, rcond, nullptr);
7061 }
7062 
7065  octave_idx_type& info, double& rcond) const
7066 {
7067  return solve (b, info, rcond, nullptr);
7068 }
7069 
7071 SparseMatrix::solve (const SparseMatrix& b, octave_idx_type& err, double& rcond,
7072  solve_singularity_handler sing_handler) const
7073 {
7074  MatrixType mattype (*this);
7075  return solve (mattype, b, err, rcond, sing_handler);
7076 }
7077 
7080 {
7081  double rcond;
7082  return solve (b, info, rcond, nullptr);
7083 }
7084 
7087  double& rcond) const
7088 {
7089  return solve (b, info, rcond, nullptr);
7090 }
7091 
7094  double& rcond,
7095  solve_singularity_handler sing_handler) const
7096 {
7097  MatrixType mattype (*this);
7098  return solve (mattype, b, err, rcond, sing_handler);
7099 }
7100 
7103 {
7104  octave_idx_type info;
7105  double rcond;
7106  return solve (b, info, rcond, nullptr);
7107 }
7108 
7111 {
7112  double rcond;
7113  return solve (b, info, rcond, nullptr);
7114 }
7115 
7118  double& rcond) const
7119 {
7120  return solve (b, info, rcond, nullptr);
7121 }
7122 
7125  octave_idx_type& err, double& rcond,
7126  solve_singularity_handler sing_handler) const
7127 {
7128  MatrixType mattype (*this);
7129  return solve (mattype, b, err, rcond, sing_handler);
7130 }
7131 
7134 {
7135  octave_idx_type info; double rcond;
7136  return solve (b, info, rcond);
7137 }
7138 
7141 {
7142  double rcond;
7143  return solve (b, info, rcond);
7144 }
7145 
7148  double& rcond) const
7149 {
7150  return solve (b, info, rcond, nullptr);
7151 }
7152 
7155  double& rcond,
7156  solve_singularity_handler sing_handler) const
7157 {
7158  Matrix tmp (b);
7159  return solve (tmp, info, rcond,
7160  sing_handler).column (static_cast<octave_idx_type> (0));
7161 }
7162 
7165 {
7166  octave_idx_type info;
7167  double rcond;
7168  return solve (b, info, rcond, nullptr);
7169 }
7170 
7173 {
7174  double rcond;
7175  return solve (b, info, rcond, nullptr);
7176 }
7177 
7180  double& rcond) const
7181 {
7182  return solve (b, info, rcond, nullptr);
7183 }
7184 
7187  double& rcond,
7188  solve_singularity_handler sing_handler) const
7189 {
7190  ComplexMatrix tmp (b);
7191  return solve (tmp, info, rcond,
7192  sing_handler).column (static_cast<octave_idx_type> (0));
7193 }
7194 
7195 // other operations.
7196 
7197 bool
7199 {
7200  octave_idx_type nel = nnz ();
7201 
7202  if (neg_zero)
7203  {
7204  for (octave_idx_type i = 0; i < nel; i++)
7205  if (lo_ieee_signbit (data (i)))
7206  return true;
7207  }
7208  else
7209  {
7210  for (octave_idx_type i = 0; i < nel; i++)
7211  if (data (i) < 0)
7212  return true;
7213  }
7214 
7215  return false;
7216 }
7217 
7218 bool
7220 {
7221  octave_idx_type nel = nnz ();
7222 
7223  for (octave_idx_type i = 0; i < nel; i++)
7224  {
7225  double val = data (i);
7226  if (octave::math::isnan (val))
7227  return true;
7228  }
7229 
7230  return false;
7231 }
7232 
7233 bool
7235 {
7236  octave_idx_type nel = nnz ();
7237 
7238  for (octave_idx_type i = 0; i < nel; i++)
7239  {
7240  double val = data (i);
7241  if (octave::math::isinf (val) || octave::math::isnan (val))
7242  return true;
7243  }
7244 
7245  return false;
7246 }
7247 
7248 bool
7250 {
7251  octave_idx_type nel = nnz ();
7252 
7253  for (octave_idx_type i = 0; i < nel; i++)
7254  {
7255  double val = data (i);
7256  if (val != 0.0 && val != 1.0)
7257  return true;
7258  }
7259 
7260  return false;
7261 }
7262 
7263 bool
7265 {
7266  octave_idx_type nel = nnz ();
7267 
7268  for (octave_idx_type i = 0; i < nel; i++)
7269  if (data (i) != 0)
7270  return false;
7271 
7272  return true;
7273 }
7274 
7275 bool
7277 {
7278  octave_idx_type nel = nnz ();
7279 
7280  for (octave_idx_type i = 0; i < nel; i++)
7281  {
7282  double val = data (i);
7283  if (octave::math::isnan (val) || octave::math::x_nint (val) == val)
7284  continue;
7285  else
7286  return false;
7287  }
7288 
7289  return true;
7290 }
7291 
7292 // Return nonzero if any element of M is not an integer. Also extract
7293 // the largest and smallest values and return them in MAX_VAL and MIN_VAL.
7294 
7295 bool
7296 SparseMatrix::all_integers (double& max_val, double& min_val) const
7297 {
7298  octave_idx_type nel = nnz ();
7299 
7300  if (nel == 0)
7301  return false;
7302 
7303  max_val = data (0);
7304  min_val = data (0);
7305 
7306  for (octave_idx_type i = 0; i < nel; i++)
7307  {
7308  double val = data (i);
7309 
7310  if (val > max_val)
7311  max_val = val;
7312 
7313  if (val < min_val)
7314  min_val = val;
7315 
7316  if (octave::math::x_nint (val) != val)
7317  return false;
7318  }
7319 
7320  return true;
7321 }
7322 
7323 bool
7325 {
7327 }
7328 
7331 {
7332  if (any_element_is_nan ())
7334 
7335  octave_idx_type nr = rows ();
7336  octave_idx_type nc = cols ();
7337  octave_idx_type nz1 = nnz ();
7338  octave_idx_type nz2 = nr*nc - nz1;
7339 
7340  SparseBoolMatrix r (nr, nc, nz2);
7341 
7342  octave_idx_type ii = 0;
7343  octave_idx_type jj = 0;
7344  r.cidx (0) = 0;
7345  for (octave_idx_type i = 0; i < nc; i++)
7346  {
7347  for (octave_idx_type j = 0; j < nr; j++)
7348  {
7349  if (jj < cidx (i+1) && ridx (jj) == j)
7350  jj++;
7351  else
7352  {
7353  r.data (ii) = true;
7354  r.ridx (ii++) = j;
7355  }
7356  }
7357  r.cidx (i+1) = ii;
7358  }
7359 
7360  return r;
7361 }
7362 
7363 // FIXME: Do these really belong here? Maybe they should be in a base class?
7364 
7366 SparseMatrix::all (int dim) const
7367 {
7368  SPARSE_ALL_OP (dim);
7369 }
7370 
7372 SparseMatrix::any (int dim) const
7373 {
7374  SPARSE_ANY_OP (dim);
7375 }
7376 
7378 SparseMatrix::cumprod (int dim) const
7379 {
7380  SPARSE_CUMPROD (SparseMatrix, double, cumprod);
7381 }
7382 
7384 SparseMatrix::cumsum (int dim) const
7385 {
7386  SPARSE_CUMSUM (SparseMatrix, double, cumsum);
7387 }
7388 
7390 SparseMatrix::prod (int dim) const
7391 {
7392  if ((rows () == 1 && dim == -1) || dim == 1)
7393  return transpose ().prod (0).transpose ();
7394  else
7395  {
7396  SPARSE_REDUCTION_OP (SparseMatrix, double, *=,
7397  (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
7398  }
7399 }
7400 
7402 SparseMatrix::sum (int dim) const
7403 {
7404  SPARSE_REDUCTION_OP (SparseMatrix, double, +=, 0.0, 0.0);
7405 }
7406 
7408 SparseMatrix::sumsq (int dim) const
7409 {
7410 #define ROW_EXPR \
7411  double d = data (i); \
7412  tmp[ridx (i)] += d * d
7413 
7414 #define COL_EXPR \
7415  double d = data (i); \
7416  tmp[j] += d * d
7417 
7419  0.0, 0.0);
7420 
7421 #undef ROW_EXPR
7422 #undef COL_EXPR
7423 }
7424 
7426 SparseMatrix::abs (void) const
7427 {
7428  octave_idx_type nz = nnz ();
7429 
7430  SparseMatrix retval (*this);
7431 
7432  for (octave_idx_type i = 0; i < nz; i++)
7433  retval.data (i) = fabs (retval.data (i));
7434 
7435  return retval;
7436 }
7437 
7440 {
7441  return MSparse<double>::diag (k);
7442 }
7443 
7444 Matrix
7446 {
7447  return Sparse<double>::array_value ();
7448 }
7449 
7450 std::ostream&
7451 operator << (std::ostream& os, const SparseMatrix& a)
7452 {
7453  octave_idx_type nc = a.cols ();
7454 
7455  // add one to the printed indices to go from
7456  // zero-based to one-based arrays
7457  for (octave_idx_type j = 0; j < nc; j++)
7458  {
7459  octave_quit ();
7460  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7461  {
7462  os << a.ridx (i) + 1 << ' ' << j + 1 << ' ';
7463  octave::write_value<double> (os, a.data (i));
7464  os << "\n";
7465  }
7466  }
7467 
7468  return os;
7469 }
7470 
7471 std::istream&
7472 operator >> (std::istream& is, SparseMatrix& a)
7473 {
7474  typedef SparseMatrix::element_type elt_type;
7475 
7476  return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7477 }
7478 
7481 {
7482  return MSparse<double>::squeeze ();
7483 }
7484 
7486 SparseMatrix::reshape (const dim_vector& new_dims) const
7487 {
7488  return MSparse<double>::reshape (new_dims);
7489 }
7490 
7492 SparseMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const
7493 {
7494  return MSparse<double>::permute (vec, inv);
7495 }
7496 
7499 {
7500  return MSparse<double>::ipermute (vec);
7501 }
7502 
7503 // matrix by matrix -> matrix operations
7504 
7507 {
7508  SPARSE_SPARSE_MUL (SparseMatrix, double, double);
7509 }
7510 
7511 Matrix
7513 {
7514  FULL_SPARSE_MUL (Matrix, double);
7515 }
7516 
7517 Matrix
7518 mul_trans (const Matrix& m, const SparseMatrix& a)
7519 {
7520  FULL_SPARSE_MUL_TRANS (Matrix, double, );
7521 }
7522 
7523 Matrix
7525 {
7526  SPARSE_FULL_MUL (Matrix, double);
7527 }
7528 
7529 Matrix
7530 trans_mul (const SparseMatrix& m, const Matrix& a)
7531 {
7532  SPARSE_FULL_TRANS_MUL (Matrix, double, );
7533 }
7534 
7535 // diag * sparse and sparse * diag
7536 
7539 {
7540  return do_mul_dm_sm<SparseMatrix> (d, a);
7541 }
7542 
7545 {
7546  return do_mul_sm_dm<SparseMatrix> (a, d);
7547 }
7548 
7551 {
7552  return do_add_dm_sm<SparseMatrix> (d, a);
7553 }
7554 
7557 {
7558  return do_sub_dm_sm<SparseMatrix> (d, a);
7559 }
7560 
7563 {
7564  return do_add_sm_dm<SparseMatrix> (a, d);
7565 }
7566 
7569 {
7570  return do_sub_sm_dm<SparseMatrix> (a, d);
7571 }
7572 
7573 // perm * sparse and sparse * perm
7574 
7577 {
7578  return octinternal_do_mul_pm_sm (p, a);
7579 }
7580 
7583 {
7584  return octinternal_do_mul_sm_pm (a, p);
7585 }
7586 
7587 // FIXME: it would be nice to share code among the min/max functions below.
7588 
7589 #define EMPTY_RETURN_CHECK(T) \
7590  if (nr == 0 || nc == 0) \
7591  return T (nr, nc);
7592 
7594 min (double d, const SparseMatrix& m)
7595 {
7596  SparseMatrix result;
7597 
7598  octave_idx_type nr = m.rows ();
7599  octave_idx_type nc = m.columns ();
7600 
7602 
7603  // Count the number of nonzero elements
7604  if (d < 0.)
7605  {
7606  result = SparseMatrix (nr, nc, d);
7607  for (octave_idx_type j = 0; j < nc; j++)
7608  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7609  {
7610  double tmp = octave::math::min (d, m.data (i));
7611  if (tmp != 0.)
7612  {
7613  octave_idx_type idx = m.ridx (i) + j * nr;
7614  result.xdata (idx) = tmp;
7615  result.xridx (idx) = m.ridx (i);
7616  }
7617  }
7618  }
7619  else
7620  {
7621  octave_idx_type nel = 0;
7622  for (octave_idx_type j = 0; j < nc; j++)
7623  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7624  if (octave::math::min (d, m.data (i)) != 0.)
7625  nel++;
7626 
7627  result = SparseMatrix (nr, nc, nel);
7628 
7629  octave_idx_type ii = 0;
7630  result.xcidx (0) = 0;
7631  for (octave_idx_type j = 0; j < nc; j++)
7632  {
7633  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7634  {
7635  double tmp = octave::math::min (d, m.data (i));
7636 
7637  if (tmp != 0.)
7638  {
7639  result.xdata (ii) = tmp;
7640  result.xridx (ii++) = m.ridx (i);
7641  }
7642  }
7643  result.xcidx (j+1) = ii;
7644  }
7645  }
7646 
7647  return result;
7648 }
7649 
7651 min (const SparseMatrix& m, double d)
7652 {
7653  return min (d, m);
7654 }
7655 
7657 min (const SparseMatrix& a, const SparseMatrix& b)
7658 {
7659  SparseMatrix r;
7660 
7661  octave_idx_type a_nr = a.rows ();
7662  octave_idx_type a_nc = a.cols ();
7663  octave_idx_type b_nr = b.rows ();
7664  octave_idx_type b_nc = b.cols ();
7665 
7666  if (a_nr == b_nr && a_nc == b_nc)
7667  {
7668  r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7669 
7670  octave_idx_type jx = 0;
7671  r.cidx (0) = 0;
7672  for (octave_idx_type i = 0 ; i < a_nc ; i++)
7673  {
7674  octave_idx_type ja = a.cidx (i);
7675  octave_idx_type ja_max = a.cidx (i+1);
7676  bool ja_lt_max = ja < ja_max;
7677 
7678  octave_idx_type jb = b.cidx (i);
7679  octave_idx_type jb_max = b.cidx (i+1);
7680  bool jb_lt_max = jb < jb_max;
7681 
7682  while (ja_lt_max || jb_lt_max)
7683  {
7684  octave_quit ();
7685  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7686  {
7687  double tmp = octave::math::min (a.data (ja), 0.);
7688  if (tmp != 0.)
7689  {
7690  r.ridx (jx) = a.ridx (ja);
7691  r.data (jx) = tmp;
7692  jx++;
7693  }
7694  ja++;
7695  ja_lt_max= ja < ja_max;
7696  }
7697  else if ((! ja_lt_max)
7698  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
7699  {
7700  double tmp = octave::math::min (0., b.data (jb));
7701  if (tmp != 0.)
7702  {
7703  r.ridx (jx) = b.ridx (jb);
7704  r.data (jx) = tmp;
7705  jx++;
7706  }
7707  jb++;
7708  jb_lt_max= jb < jb_max;
7709  }
7710  else
7711  {
7712  double tmp = octave::math::min (a.data (ja), b.data (jb));
7713  if (tmp != 0.)
7714  {
7715  r.data (jx) = tmp;
7716  r.ridx (jx) = a.ridx (ja);
7717  jx++;
7718  }
7719  ja++;
7720  ja_lt_max= ja < ja_max;
7721  jb++;
7722  jb_lt_max= jb < jb_max;
7723  }
7724  }
7725  r.cidx (i+1) = jx;
7726  }
7727 
7728  r.maybe_compress ();
7729  }
7730  else
7731  {
7732  if (a_nr == 0 || a_nc == 0)
7733  r.resize (a_nr, a_nc);
7734  else if (b_nr == 0 || b_nc == 0)
7735  r.resize (b_nr, b_nc);
7736  else
7737  octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
7738  }
7739 
7740  return r;
7741 }
7742 
7744 max (double d, const SparseMatrix& m)
7745 {
7746  SparseMatrix result;
7747 
7748  octave_idx_type nr = m.rows ();
7749  octave_idx_type nc = m.columns ();
7750 
7752 
7753  // Count the number of nonzero elements
7754  if (d > 0.)
7755  {
7756  result = SparseMatrix (nr, nc, d);
7757  for (octave_idx_type j = 0; j < nc; j++)
7758  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7759  {
7760  double tmp = octave::math::max (d, m.data (i));
7761 
7762  if (tmp != 0.)
7763  {
7764  octave_idx_type idx = m.ridx (i) + j * nr;
7765  result.xdata (idx) = tmp;
7766  result.xridx (idx) = m.ridx (i);
7767  }
7768  }
7769  }
7770  else
7771  {
7772  octave_idx_type nel = 0;
7773  for (octave_idx_type j = 0; j < nc; j++)
7774  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7775  if (octave::math::max (d, m.data (i)) != 0.)
7776  nel++;
7777 
7778  result = SparseMatrix (nr, nc, nel);
7779 
7780  octave_idx_type ii = 0;
7781  result.xcidx (0) = 0;
7782  for (octave_idx_type j = 0; j < nc; j++)
7783  {
7784  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7785  {
7786  double tmp = octave::math::max (d, m.data (i));
7787  if (tmp != 0.)
7788  {
7789  result.xdata (ii) = tmp;
7790  result.xridx (ii++) = m.ridx (i);
7791  }
7792  }
7793  result.xcidx (j+1) = ii;
7794  }
7795  }
7796 
7797  return result;
7798 }
7799 
7801 max (const SparseMatrix& m, double d)
7802 {
7803  return max (d, m);
7804 }
7805 
7807 max (const SparseMatrix& a, const SparseMatrix& b)
7808 {
7809  SparseMatrix r;
7810 
7811  octave_idx_type a_nr = a.rows ();
7812  octave_idx_type a_nc = a.cols ();
7813  octave_idx_type b_nr = b.rows ();
7814  octave_idx_type b_nc = b.cols ();
7815 
7816  if (a_nr == b_nr && a_nc == b_nc)
7817  {
7818  r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7819 
7820  octave_idx_type jx = 0;
7821  r.cidx (0) = 0;
7822  for (octave_idx_type i = 0 ; i < a_nc ; i++)
7823  {
7824  octave_idx_type ja = a.cidx (i);
7825  octave_idx_type ja_max = a.cidx (i+1);
7826  bool ja_lt_max = ja < ja_max;
7827 
7828  octave_idx_type jb = b.cidx (i);
7829  octave_idx_type jb_max = b.cidx (i+1);
7830  bool jb_lt_max = jb < jb_max;
7831 
7832  while (ja_lt_max || jb_lt_max)
7833  {
7834  octave_quit ();
7835  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7836  {
7837  double tmp = octave::math::max (a.data (ja), 0.);
7838  if (tmp != 0.)
7839  {
7840  r.ridx (jx) = a.ridx (ja);
7841  r.data (jx) = tmp;
7842  jx++;
7843  }
7844  ja++;
7845  ja_lt_max= ja < ja_max;
7846  }
7847  else if ((! ja_lt_max)
7848  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
7849  {
7850  double tmp = octave::math::max (0., b.data (jb));
7851  if (tmp != 0.)
7852  {
7853  r.ridx (jx) = b.ridx (jb);
7854  r.data (jx) = tmp;
7855  jx++;
7856  }
7857  jb++;
7858  jb_lt_max= jb < jb_max;
7859  }
7860  else
7861  {
7862  double tmp = octave::math::max (a.data (ja), b.data (jb));
7863  if (tmp != 0.)
7864  {
7865  r.data (jx) = tmp;
7866  r.ridx (jx) = a.ridx (ja);
7867  jx++;
7868  }
7869  ja++;
7870  ja_lt_max= ja < ja_max;
7871  jb++;
7872  jb_lt_max= jb < jb_max;
7873  }
7874  }
7875  r.cidx (i+1) = jx;
7876  }
7877 
7878  r.maybe_compress ();
7879  }
7880  else
7881  {
7882  if (a_nr == 0 || a_nc == 0)
7883  r.resize (a_nr, a_nc);
7884  else if (b_nr == 0 || b_nc == 0)
7885  r.resize (b_nr, b_nc);
7886  else
7887  octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
7888  }
7889 
7890  return r;
7891 }
7892 
7895 
7898 
base_det< double > DET
Definition: DET.h:93
#define Inf
Definition: Faddeeva.cc:260
#define NaN
Definition: Faddeeva.cc:261
#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)
OCTARRAY_OVERRIDABLE_FUNC_API const T * data(void) const
Size of the specified dimension.
Definition: Array.h:663
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array-base.cc:1032
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type rows(void) const
Definition: Array.h:459
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array-base.cc:1766
OCTARRAY_OVERRIDABLE_FUNC_API T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
OCTARRAY_OVERRIDABLE_FUNC_API octave_idx_type cols(void) const
Definition: Array.h:469
OCTARRAY_OVERRIDABLE_FUNC_API T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:524
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:193
octave_idx_type length(void) const
Definition: DiagArray2.h:95
octave_idx_type cols(void) const
Definition: DiagArray2.h:90
MSparse< T > permute(const Array< octave_idx_type > &vec, bool inv=false) const
Definition: MSparse.h:105
MSparse< T > squeeze(void) const
Definition: MSparse.h:100
MSparse< T > ipermute(const Array< octave_idx_type > &vec) const
Definition: MSparse.h:108
MSparse< T > reshape(const dim_vector &new_dims) const
Definition: MSparse.h:102
MSparse< T > & insert(const Sparse< T > &a, octave_idx_type r, octave_idx_type c)
Definition: MSparse.h:86
MSparse< T > diag(octave_idx_type k=0) const
Definition: MSparse.h:111
bool ishermitian(void) const
Definition: MatrixType.h:122
OCTAVE_API MatrixType transpose(void) const
Definition: MatrixType.cc:973
void mark_as_rectangular(void)
Definition: MatrixType.h:155
bool is_dense(void) const
Definition: MatrixType.h:105
@ Tridiagonal_Hermitian
Definition: MatrixType.h:53
@ Permuted_Lower
Definition: MatrixType.h:48
@ Banded_Hermitian
Definition: MatrixType.h:51
@ Permuted_Diagonal
Definition: MatrixType.h:44
@ Permuted_Upper
Definition: MatrixType.h:47
OCTAVE_API void mark_as_unsymmetric(void)
Definition: MatrixType.cc:920
octave_idx_type * triangular_perm(void) const
Definition: MatrixType.h:136
int nupper(void) const
Definition: MatrixType.h:101
int nlower(void) const
Definition: MatrixType.h:103
OCTAVE_API void info(void) const
Definition: MatrixType.cc:846
OCTAVE_API int type(bool quiet=true)
Definition: MatrixType.cc:656
Definition: dMatrix.h:42
OCTAVE_API Matrix sum(int dim=-1) const
Definition: dMatrix.cc:2389
OCTAVE_API RowVector row(octave_idx_type i) const
Definition: dMatrix.cc:416
OCTAVE_API Matrix abs(void) const
Definition: dMatrix.cc:2401
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
OCTAVE_API ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:422
OCTAVE_API double max(void) const
Definition: dRowVector.cc:223
OCTAVE_API SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
Definition: CSparse.cc:556
OCTAVE_API bool issymmetric(void) const
Definition: dSparse.cc:131
OCTAVE_API SparseMatrix dinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
Definition: dSparse.cc:627
OCTAVE_API SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
Definition: dSparse.cc:171
SparseMatrix(void)
Definition: dSparse.h:53
OCTAVE_API bool any_element_is_inf_or_nan(void) const
Definition: dSparse.cc:7234
OCTAVE_API 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
OCTAVE_API SparseMatrix reshape(const dim_vector &new_dims) const
Definition: dSparse.cc:7486
OCTAVE_API SparseMatrix sum(int dim=-1) const
Definition: dSparse.cc:7402
OCTAVE_API SparseBoolMatrix any(int dim=-1) const
Definition: dSparse.cc:7372
OCTAVE_API 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
OCTAVE_API 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
OCTAVE_API bool all_integers(double &max_val, double &min_val) const
Definition: dSparse.cc:7296
OCTAVE_API bool operator==(const SparseMatrix &a) const
Definition: dSparse.cc:101
OCTAVE_API SparseMatrix abs(void) const
Definition: dSparse.cc:7426
OCTAVE_API bool operator!=(const SparseMatrix &a) const
Definition: dSparse.cc:125
OCTAVE_API RowVector row(octave_idx_type i) const
Definition: dSparse.cc:503
OCTAVE_API SparseMatrix cumprod(int dim=-1) const
Definition: dSparse.cc:7378
OCTAVE_API SparseMatrix diag(octave_idx_type k=0) const
Definition: dSparse.cc:7439
OCTAVE_API SparseMatrix prod(int dim=-1) const
Definition: dSparse.cc:7390
OCTAVE_API SparseBoolMatrix all(int dim=-1) const
Definition: dSparse.cc:7366
OCTAVE_API Matrix matrix_value(void) const
Definition: dSparse.cc:7445
OCTAVE_API SparseMatrix min(int dim=-1) const
Definition: dSparse.cc:337
OCTAVE_API bool any_element_is_nan(void) const
Definition: dSparse.cc:7219
OCTAVE_API SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
Definition: dSparse.cc:7492
OCTAVE_API SparseMatrix cumsum(int dim=-1) const
Definition: dSparse.cc:7384
OCTAVE_API 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
OCTAVE_API ColumnVector column(octave_idx_type i) const
Definition: dSparse.cc:522
SparseMatrix transpose(void) const
Definition: dSparse.h:124
OCTAVE_API SparseMatrix inverse(void) const
Definition: dSparse.cc:603
OCTAVE_API SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
Definition: dSparse.cc:534
OCTAVE_API bool too_large_for_float(void) const
Definition: dSparse.cc:7324
OCTAVE_API SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
Definition: dSparse.cc:7498
OCTAVE_API bool any_element_not_one_or_zero(void) const
Definition: dSparse.cc:7249
OCTAVE_API DET determinant(void) const
Definition: dSparse.cc:1012
OCTAVE_API bool all_elements_are_zero(void) const
Definition: dSparse.cc:7264
OCTAVE_API SparseMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcond, const bool force=false, const bool calccond=true) const
Definition: dSparse.cc:677
OCTAVE_API SparseMatrix sumsq(int dim=-1) const
Definition: dSparse.cc:7408
OCTAVE_API SparseMatrix squeeze(void) const
Definition: dSparse.cc:7480
OCTAVE_API Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6692
OCTAVE_API 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
OCTAVE_API bool all_elements_are_int_or_inf_or_nan(void) const
Definition: dSparse.cc:7276
OCTAVE_API 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
OCTAVE_API SparseBoolMatrix operator!(void) const
Definition: dSparse.cc:7330
OCTAVE_API SparseMatrix max(int dim=-1) const
Definition: dSparse.cc:186
OCTAVE_API bool any_element_is_negative(bool=false) const
Definition: dSparse.cc:7198
OCTAVE_API 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
octave_idx_type rows(void) const
Definition: Sparse.h:351
octave_idx_type * cidx(void)
Definition: Sparse.h:596
octave_idx_type columns(void) const
Definition: Sparse.h:353
T * xdata(void)
Definition: Sparse.h:576
octave_idx_type * ridx(void)
Definition: Sparse.h:583
T * data(void)
Definition: Sparse.h:574
octave_idx_type nnz(void) const
Actual number of nonzero terms.
Definition: Sparse.h:339
OCTAVE_API Array< T > array_value(void) const
Definition: Sparse.cc:2766
dim_vector dims(void) const
Definition: Sparse.h:371
bool test_any(F fcn) const
Definition: Sparse.h:663
T element_type
Definition: Sparse.h:52
Sparse< T, Alloc > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:531
T & elem(octave_idx_type n)
Definition: Sparse.h:456
octave_idx_type cols(void) const
Definition: Sparse.h:352
void change_capacity(octave_idx_type nz)
Definition: Sparse.h:556
octave_idx_type * xcidx(void)
Definition: Sparse.h:602
octave_idx_type * xridx(void)
Definition: Sparse.h:589
Definition: DET.h:39
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
int first_non_singleton(int def=0) const
Definition: dim-vector.h:444
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:257
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
Definition: dSparse.cc:7530
SparseMatrix real(const SparseComplexMatrix &a)
Definition: dSparse.cc:554
#define EMPTY_RETURN_CHECK(T)
Definition: dSparse.cc:7589
SparseMatrix min(double d, const SparseMatrix &m)
Definition: dSparse.cc:7594
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
Definition: dSparse.cc:7518
#define COL_EXPR
#define ROW_EXPR
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
Definition: dSparse.cc:7556
SparseMatrix max(double d, const SparseMatrix &m)
Definition: dSparse.cc:7744
std::istream & operator>>(std::istream &is, SparseMatrix &a)
Definition: dSparse.cc:7472
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
Definition: dSparse.cc:7506
SparseMatrix imag(const SparseComplexMatrix &a)
Definition: dSparse.cc:575
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
Definition: dSparse.cc:7550
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)
Definition: dSparse.cc:7451
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:316
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:45
octave_f77_int_type F77_INT
Definition: f77-fcn.h:306
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void err_nan_to_logical_conversion(void)
void warn_singular_matrix(double rcond)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:41
#define lo_ieee_signbit(x)
Definition: lo-ieee.h:120
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
bool isinf(double x)
Definition: lo-mappers.h:203
bool isnan(bool)
Definition: lo-mappers.h:178
T x_nint(T x)
Definition: lo-mappers.h:269
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
bool too_large_for_float(double x)
Definition: lo-utils.cc:55
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
class OCTAVE_API SparseMatrix
Definition: mx-fwd.h:55
class OCTAVE_API SparseComplexMatrix
Definition: mx-fwd.h:56
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
std::complex< double > Complex
Definition: oct-cmplx.h:33
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition: oct-sparse.cc:51
#define CHOLMOD_NAME(name)
Definition: oct-sparse.h:140
#define UMFPACK_DNAME(name)
Definition: oct-sparse.h:169
const octave_base_value const Array< octave_idx_type > & ra_idx
template class OCTAVE_API sparse_chol< SparseMatrix >
Definition: sparse-chol.cc:561
template OCTAVE_API SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
template OCTAVE_API ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
template OCTAVE_API Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
template OCTAVE_API SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
template class OCTAVE_API sparse_lu< SparseMatrix >
Definition: sparse-lu.cc:984
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition: sparse-qr.cc:3288
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