GNU Octave  9.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-2024 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 
627 SparseMatrix::dinverse (MatrixType& mattype, octave_idx_type& info,
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 
677 SparseMatrix::tinverse (MatrixType& mattype, octave_idx_type& info,
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);
982  octave::math::sparse_lu<SparseMatrix> fact (*this,
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
1159 SparseMatrix::dsolve (MatrixType& mattype, const Matrix& b,
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 
1219 SparseMatrix::dsolve (MatrixType& mattype, const SparseMatrix& b,
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 
1308 SparseMatrix::dsolve (MatrixType& mattype, const ComplexMatrix& b,
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 
1367 SparseMatrix::dsolve (MatrixType& mattype, const SparseComplexMatrix& b,
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
1456 SparseMatrix::utsolve (MatrixType& mattype, const Matrix& b,
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 
1685 SparseMatrix::utsolve (MatrixType& mattype, const SparseMatrix& b,
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 
1966 SparseMatrix::utsolve (MatrixType& mattype, const ComplexMatrix& b,
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 
2198 SparseMatrix::utsolve (MatrixType& mattype, const SparseComplexMatrix& b,
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
2482 SparseMatrix::ltsolve (MatrixType& mattype, const Matrix& b,
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 
2735 SparseMatrix::ltsolve (MatrixType& mattype, const SparseMatrix& b,
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 
3035 SparseMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
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 
3289 SparseMatrix::ltsolve (MatrixType& mattype, const SparseComplexMatrix& b,
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
3591 SparseMatrix::trisolve (MatrixType& mattype, const Matrix& b,
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 
3757 SparseMatrix::trisolve (MatrixType& mattype, const SparseMatrix& b,
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 
3911 SparseMatrix::trisolve (MatrixType& mattype, const ComplexMatrix& b,
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 
4081 SparseMatrix::trisolve (MatrixType& mattype, const SparseComplexMatrix& b,
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
4271 SparseMatrix::bsolve (MatrixType& mattype, const Matrix& b,
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 
4540 SparseMatrix::bsolve (MatrixType& mattype, const SparseMatrix& b,
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 
4880 SparseMatrix::bsolve (MatrixType& mattype, const ComplexMatrix& b,
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 
5204 SparseMatrix::bsolve (MatrixType& mattype, const SparseComplexMatrix& b,
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
5717 SparseMatrix::fsolve (MatrixType& mattype, const Matrix& b,
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 
5928 SparseMatrix::fsolve (MatrixType& mattype, const SparseMatrix& b,
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
6057  (octave::from_size_t (X->nrow),
6058  octave::from_size_t (X->ncol),
6060  (X, cm)));
6061  for (octave_idx_type j = 0;
6062  j <= static_cast<octave_idx_type> (X->ncol); j++)
6063  retval.xcidx (j) = static_cast<octave_idx_type *> (X->p)[j];
6064  for (octave_idx_type j = 0;
6066  j++)
6067  {
6068  retval.xridx (j) = static_cast<octave_idx_type *> (X->i)[j];
6069  retval.xdata (j) = static_cast<double *> (X->x)[j];
6070  }
6071 
6072  CHOLMOD_NAME(free_sparse) (&X, cm);
6073  CHOLMOD_NAME(free_factor) (&L, cm);
6074  CHOLMOD_NAME(finish) (cm);
6075  static char blank_name[] = " ";
6076  CHOLMOD_NAME(print_common) (blank_name, cm);
6077  }
6078 #else
6079  (*current_liboctave_warning_with_id_handler)
6080  ("Octave:missing-dependency",
6081  "support for CHOLMOD was unavailable or disabled "
6082  "when liboctave was built");
6083 
6084  mattype.mark_as_unsymmetric ();
6085  typ = MatrixType::Full;
6086 #endif
6087  }
6088 
6089  if (typ == MatrixType::Full)
6090  {
6091 #if defined (HAVE_UMFPACK)
6092  Matrix Control, Info;
6093  void *Numeric = factorize (err, rcond, Control, Info,
6094  sing_handler, calc_cond);
6095 
6096  if (err == 0)
6097  {
6098  // one iterative refinement instead of the default two in UMFPACK
6099  Control (UMFPACK_IRSTEP) = 1;
6100  octave_idx_type b_nr = b.rows ();
6101  octave_idx_type b_nc = b.cols ();
6102  int status = 0;
6103  double *control = Control.fortran_vec ();
6104  double *info = Info.fortran_vec ();
6105  const octave_idx_type *Ap = cidx ();
6106  const octave_idx_type *Ai = ridx ();
6107  const double *Ax = data ();
6108 
6109  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6110  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6111 
6112  // Take a first guess that the number of nonzero terms
6113  // will be as many as in b
6114  octave_idx_type x_nz = b.nnz ();
6115  octave_idx_type ii = 0;
6116  retval = SparseMatrix (b_nr, b_nc, x_nz);
6117 
6118  retval.xcidx (0) = 0;
6119  for (octave_idx_type j = 0; j < b_nc; j++)
6120  {
6121 
6122  for (octave_idx_type i = 0; i < b_nr; i++)
6123  Bx[i] = b.elem (i, j);
6124 
6125  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6128  Ax, Xx, Bx, Numeric,
6129  control, info);
6130  if (status < 0)
6131  {
6132  UMFPACK_DNAME (report_status) (control, status);
6133 
6134  // FIXME: Should this be a warning?
6135  (*current_liboctave_error_handler)
6136  ("SparseMatrix::solve solve failed");
6137 
6138  err = -1;
6139  break;
6140  }
6141 
6142  for (octave_idx_type i = 0; i < b_nr; i++)
6143  {
6144  double tmp = Xx[i];
6145  if (tmp != 0.0)
6146  {
6147  if (ii == x_nz)
6148  {
6149  // Resize the sparse matrix
6150  octave_idx_type sz;
6151  sz = (static_cast<double> (b_nc) - j) / b_nc
6152  * x_nz;
6153  sz = x_nz + (sz > 100 ? sz : 100);
6154  retval.change_capacity (sz);
6155  x_nz = sz;
6156  }
6157  retval.xdata (ii) = tmp;
6158  retval.xridx (ii++) = i;
6159  }
6160  }
6161  retval.xcidx (j+1) = ii;
6162  }
6163 
6164  retval.maybe_compress ();
6165 
6166  UMFPACK_DNAME (report_info) (control, info);
6167 
6168  UMFPACK_DNAME (free_numeric) (&Numeric);
6169  }
6170  else
6171  mattype.mark_as_rectangular ();
6172 
6173 #else
6174  octave_unused_parameter (rcond);
6175  octave_unused_parameter (sing_handler);
6176  octave_unused_parameter (calc_cond);
6177 
6178  (*current_liboctave_error_handler)
6179  ("support for UMFPACK was unavailable or disabled "
6180  "when liboctave was built");
6181 #endif
6182  }
6183  else if (typ != MatrixType::Hermitian)
6184  (*current_liboctave_error_handler) ("incorrect matrix type");
6185  }
6186 
6187  return retval;
6188 }
6189 
6191 SparseMatrix::fsolve (MatrixType& mattype, const ComplexMatrix& b,
6192  octave_idx_type& err, double& rcond,
6193  solve_singularity_handler sing_handler,
6194  bool calc_cond) const
6195 {
6196  ComplexMatrix retval;
6197 
6198  octave_idx_type nr = rows ();
6199  octave_idx_type nc = cols ();
6200  err = 0;
6201 
6202  if (nr != nc || nr != b.rows ())
6203  (*current_liboctave_error_handler)
6204  ("matrix dimension mismatch solution of linear equations");
6205 
6206  if (nr == 0 || b.cols () == 0)
6207  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
6208  else
6209  {
6210  // Print spparms("spumoni") info if requested
6211  volatile int typ = mattype.type ();
6212  mattype.info ();
6213 
6214  if (typ == MatrixType::Hermitian)
6215  {
6216 #if defined (HAVE_CHOLMOD)
6217  cholmod_common Common;
6218  cholmod_common *cm = &Common;
6219 
6220  // Setup initial parameters
6221  CHOLMOD_NAME(start) (cm);
6222  cm->prefer_zomplex = false;
6223 
6224  double spu = octave::sparse_params::get_key ("spumoni");
6225  if (spu == 0.)
6226  {
6227  cm->print = -1;
6228  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6229  }
6230  else
6231  {
6232  cm->print = static_cast<int> (spu) + 2;
6233  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6234  }
6235 
6236  cm->error_handler = &SparseCholError;
6237  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6238  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6239 
6240  cm->final_ll = true;
6241 
6242  cholmod_sparse Astore;
6243  cholmod_sparse *A = &Astore;
6244  A->nrow = nr;
6245  A->ncol = nc;
6246 
6247  A->p = cidx ();
6248  A->i = ridx ();
6249  A->nzmax = nnz ();
6250  A->packed = true;
6251  A->sorted = true;
6252  A->nz = nullptr;
6253 #if defined (OCTAVE_ENABLE_64)
6254  A->itype = CHOLMOD_LONG;
6255 #else
6256  A->itype = CHOLMOD_INT;
6257 #endif
6258  A->dtype = CHOLMOD_DOUBLE;
6259  A->stype = 1;
6260  A->xtype = CHOLMOD_REAL;
6261 
6262  A->x = data ();
6263 
6264  cholmod_dense Bstore;
6265  cholmod_dense *B = &Bstore;
6266  B->nrow = b.rows ();
6267  B->ncol = b.cols ();
6268  B->d = B->nrow;
6269  B->nzmax = B->nrow * B->ncol;
6270  B->dtype = CHOLMOD_DOUBLE;
6271  B->xtype = CHOLMOD_COMPLEX;
6272 
6273  B->x = const_cast<Complex *> (b.data ());
6274 
6275  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6276  CHOLMOD_NAME(factorize) (A, L, cm);
6277  if (calc_cond)
6278  rcond = CHOLMOD_NAME(rcond)(L, cm);
6279  else
6280  rcond = 1.0;
6281 
6282  if (rcond == 0.0)
6283  {
6284  // Either its indefinite or singular. Try UMFPACK
6285  mattype.mark_as_unsymmetric ();
6286  typ = MatrixType::Full;
6287  }
6288  else
6289  {
6290  volatile double rcond_plus_one = rcond + 1.0;
6291 
6292  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
6293  {
6294  err = -2;
6295 
6296  if (sing_handler)
6297  {
6298  sing_handler (rcond);
6299  mattype.mark_as_rectangular ();
6300  }
6301  else
6303 
6304  return retval;
6305  }
6306 
6307  cholmod_dense *X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
6308 
6309  retval.resize (b.rows (), b.cols ());
6310  for (octave_idx_type j = 0; j < b.cols (); j++)
6311  {
6312  octave_idx_type jr = j * b.rows ();
6313  for (octave_idx_type i = 0; i < b.rows (); i++)
6314  retval.xelem (i, j) = static_cast<Complex *> (X->x)[jr + i];
6315  }
6316 
6317  CHOLMOD_NAME(free_dense) (&X, cm);
6318  CHOLMOD_NAME(free_factor) (&L, cm);
6319  CHOLMOD_NAME(finish) (cm);
6320  static char blank_name[] = " ";
6321  CHOLMOD_NAME(print_common) (blank_name, cm);
6322  }
6323 #else
6324  (*current_liboctave_warning_with_id_handler)
6325  ("Octave:missing-dependency",
6326  "support for CHOLMOD was unavailable or disabled "
6327  "when liboctave was built");
6328 
6329  mattype.mark_as_unsymmetric ();
6330  typ = MatrixType::Full;
6331 #endif
6332  }
6333 
6334  if (typ == MatrixType::Full)
6335  {
6336 #if defined (HAVE_UMFPACK)
6337  Matrix Control, Info;
6338  void *Numeric = factorize (err, rcond, Control, Info,
6339  sing_handler, calc_cond);
6340 
6341  if (err == 0)
6342  {
6343  // one iterative refinement instead of the default two in UMFPACK
6344  Control (UMFPACK_IRSTEP) = 1;
6345  octave_idx_type b_nr = b.rows ();
6346  octave_idx_type b_nc = b.cols ();
6347  int status = 0;
6348  double *control = Control.fortran_vec ();
6349  double *info = Info.fortran_vec ();
6350  const octave_idx_type *Ap = cidx ();
6351  const octave_idx_type *Ai = ridx ();
6352  const double *Ax = data ();
6353 
6354  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6355  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6356 
6357  retval.resize (b_nr, b_nc);
6358 
6359  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6360  OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6361 
6362  for (octave_idx_type j = 0; j < b_nc; j++)
6363  {
6364  for (octave_idx_type i = 0; i < b_nr; i++)
6365  {
6366  Complex c = b(i, j);
6367  Bx[i] = c.real ();
6368  Bz[i] = c.imag ();
6369  }
6370 
6371  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6374  Ax, Xx, Bx, Numeric,
6375  control, info);
6376  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6379  Ax, Xz, Bz,
6380  Numeric, control, info);
6381 
6382  if (status < 0 || status2 < 0)
6383  {
6384  UMFPACK_DNAME (report_status) (control, status);
6385 
6386  // FIXME: Should this be a warning?
6387  (*current_liboctave_error_handler)
6388  ("SparseMatrix::solve solve failed");
6389 
6390  err = -1;
6391  break;
6392  }
6393 
6394  for (octave_idx_type i = 0; i < b_nr; i++)
6395  retval(i, j) = Complex (Xx[i], Xz[i]);
6396  }
6397 
6398  UMFPACK_DNAME (report_info) (control, info);
6399 
6400  UMFPACK_DNAME (free_numeric) (&Numeric);
6401  }
6402  else
6403  mattype.mark_as_rectangular ();
6404 
6405 #else
6406  octave_unused_parameter (rcond);
6407  octave_unused_parameter (sing_handler);
6408  octave_unused_parameter (calc_cond);
6409 
6410  (*current_liboctave_error_handler)
6411  ("support for UMFPACK was unavailable or disabled "
6412  "when liboctave was built");
6413 #endif
6414  }
6415  else if (typ != MatrixType::Hermitian)
6416  (*current_liboctave_error_handler) ("incorrect matrix type");
6417  }
6418 
6419  return retval;
6420 }
6421 
6423 SparseMatrix::fsolve (MatrixType& mattype, const SparseComplexMatrix& b,
6424  octave_idx_type& err, double& rcond,
6425  solve_singularity_handler sing_handler,
6426  bool calc_cond) const
6427 {
6428  SparseComplexMatrix retval;
6429 
6430  octave_idx_type nr = rows ();
6431  octave_idx_type nc = cols ();
6432  err = 0;
6433 
6434  if (nr != nc || nr != b.rows ())
6435  (*current_liboctave_error_handler)
6436  ("matrix dimension mismatch solution of linear equations");
6437 
6438  if (nr == 0 || b.cols () == 0)
6439  retval = SparseComplexMatrix (nc, b.cols ());
6440  else
6441  {
6442  // Print spparms("spumoni") info if requested
6443  volatile int typ = mattype.type ();
6444  mattype.info ();
6445 
6446  if (typ == MatrixType::Hermitian)
6447  {
6448 #if defined (HAVE_CHOLMOD)
6449  cholmod_common Common;
6450  cholmod_common *cm = &Common;
6451 
6452  // Setup initial parameters
6453  CHOLMOD_NAME(start) (cm);
6454  cm->prefer_zomplex = false;
6455 
6456  double spu = octave::sparse_params::get_key ("spumoni");
6457  if (spu == 0.)
6458  {
6459  cm->print = -1;
6460  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr);
6461  }
6462  else
6463  {
6464  cm->print = static_cast<int> (spu) + 2;
6465  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
6466  }
6467 
6468  cm->error_handler = &SparseCholError;
6469  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
6470  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
6471 
6472  cm->final_ll = true;
6473 
6474  cholmod_sparse Astore;
6475  cholmod_sparse *A = &Astore;
6476  A->nrow = nr;
6477  A->ncol = nc;
6478 
6479  A->p = cidx ();
6480  A->i = ridx ();
6481  A->nzmax = nnz ();
6482  A->packed = true;
6483  A->sorted = true;
6484  A->nz = nullptr;
6485 #if defined (OCTAVE_ENABLE_64)
6486  A->itype = CHOLMOD_LONG;
6487 #else
6488  A->itype = CHOLMOD_INT;
6489 #endif
6490  A->dtype = CHOLMOD_DOUBLE;
6491  A->stype = 1;
6492  A->xtype = CHOLMOD_REAL;
6493 
6494  A->x = data ();
6495 
6496  cholmod_sparse Bstore;
6497  cholmod_sparse *B = &Bstore;
6498  B->nrow = b.rows ();
6499  B->ncol = b.cols ();
6500  B->p = b.cidx ();
6501  B->i = b.ridx ();
6502  B->nzmax = b.nnz ();
6503  B->packed = true;
6504  B->sorted = true;
6505  B->nz = nullptr;
6506 #if defined (OCTAVE_ENABLE_64)
6507  B->itype = CHOLMOD_LONG;
6508 #else
6509  B->itype = CHOLMOD_INT;
6510 #endif
6511  B->dtype = CHOLMOD_DOUBLE;
6512  B->stype = 0;
6513  B->xtype = CHOLMOD_COMPLEX;
6514 
6515  B->x = b.data ();
6516 
6517  cholmod_factor *L = CHOLMOD_NAME(analyze) (A, cm);
6518  CHOLMOD_NAME(factorize) (A, L, cm);
6519  if (calc_cond)
6520  rcond = CHOLMOD_NAME(rcond)(L, cm);
6521  else
6522  rcond = 1.0;
6523 
6524  if (rcond == 0.0)
6525  {
6526  // Either its indefinite or singular. Try UMFPACK
6527  mattype.mark_as_unsymmetric ();
6528  typ = MatrixType::Full;
6529  }
6530  else
6531  {
6532  volatile double rcond_plus_one = rcond + 1.0;
6533 
6534  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
6535  {
6536  err = -2;
6537 
6538  if (sing_handler)
6539  {
6540  sing_handler (rcond);
6541  mattype.mark_as_rectangular ();
6542  }
6543  else
6545 
6546  return retval;
6547  }
6548 
6549  cholmod_sparse *X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6550 
6551  retval = SparseComplexMatrix
6552  (octave::from_size_t (X->nrow),
6553  octave::from_size_t (X->ncol),
6555  (X, cm)));
6556  for (octave_idx_type j = 0;
6557  j <= static_cast<octave_idx_type> (X->ncol); j++)
6558  retval.xcidx (j) = static_cast<octave_idx_type *> (X->p)[j];
6559  for (octave_idx_type j = 0;
6561  j++)
6562  {
6563  retval.xridx (j) = static_cast<octave_idx_type *> (X->i)[j];
6564  retval.xdata (j) = static_cast<Complex *> (X->x)[j];
6565  }
6566 
6567  CHOLMOD_NAME(free_sparse) (&X, cm);
6568  CHOLMOD_NAME(free_factor) (&L, cm);
6569  CHOLMOD_NAME(finish) (cm);
6570  static char blank_name[] = " ";
6571  CHOLMOD_NAME(print_common) (blank_name, cm);
6572  }
6573 #else
6574  (*current_liboctave_warning_with_id_handler)
6575  ("Octave:missing-dependency",
6576  "support for CHOLMOD was unavailable or disabled "
6577  "when liboctave was built");
6578 
6579  mattype.mark_as_unsymmetric ();
6580  typ = MatrixType::Full;
6581 #endif
6582  }
6583 
6584  if (typ == MatrixType::Full)
6585  {
6586 #if defined (HAVE_UMFPACK)
6587  Matrix Control, Info;
6588  void *Numeric = factorize (err, rcond, Control, Info,
6589  sing_handler, calc_cond);
6590 
6591  if (err == 0)
6592  {
6593  // one iterative refinement instead of the default two in UMFPACK
6594  Control (UMFPACK_IRSTEP) = 1;
6595  octave_idx_type b_nr = b.rows ();
6596  octave_idx_type b_nc = b.cols ();
6597  int status = 0;
6598  double *control = Control.fortran_vec ();
6599  double *info = Info.fortran_vec ();
6600  const octave_idx_type *Ap = cidx ();
6601  const octave_idx_type *Ai = ridx ();
6602  const double *Ax = data ();
6603 
6604  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6605  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6606 
6607  // Take a first guess that the number of nonzero terms
6608  // will be as many as in b
6609  octave_idx_type x_nz = b.nnz ();
6610  octave_idx_type ii = 0;
6611  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6612 
6613  OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
6614  OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
6615 
6616  retval.xcidx (0) = 0;
6617  for (octave_idx_type j = 0; j < b_nc; j++)
6618  {
6619  for (octave_idx_type i = 0; i < b_nr; i++)
6620  {
6621  Complex c = b(i, j);
6622  Bx[i] = c.real ();
6623  Bz[i] = c.imag ();
6624  }
6625 
6626  status = UMFPACK_DNAME (solve) (UMFPACK_A,
6629  Ax, Xx, Bx, Numeric,
6630  control, info);
6631  int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
6634  Ax, Xz, Bz,
6635  Numeric, control, info);
6636 
6637  if (status < 0 || status2 < 0)
6638  {
6639  UMFPACK_DNAME (report_status) (control, status);
6640 
6641  // FIXME: Should this be a warning?
6642  (*current_liboctave_error_handler)
6643  ("SparseMatrix::solve solve failed");
6644 
6645  err = -1;
6646  break;
6647  }
6648 
6649  for (octave_idx_type i = 0; i < b_nr; i++)
6650  {
6651  Complex tmp = Complex (Xx[i], Xz[i]);
6652  if (tmp != 0.0)
6653  {
6654  if (ii == x_nz)
6655  {
6656  // Resize the sparse matrix
6657  octave_idx_type sz;
6658  sz = (static_cast<double> (b_nc) - j) / b_nc
6659  * x_nz;
6660  sz = x_nz + (sz > 100 ? sz : 100);
6661  retval.change_capacity (sz);
6662  x_nz = sz;
6663  }
6664  retval.xdata (ii) = tmp;
6665  retval.xridx (ii++) = i;
6666  }
6667  }
6668  retval.xcidx (j+1) = ii;
6669  }
6670 
6671  retval.maybe_compress ();
6672 
6673  UMFPACK_DNAME (report_info) (control, info);
6674 
6675  UMFPACK_DNAME (free_numeric) (&Numeric);
6676  }
6677  else
6678  mattype.mark_as_rectangular ();
6679 #else
6680  octave_unused_parameter (rcond);
6681  octave_unused_parameter (sing_handler);
6682  octave_unused_parameter (calc_cond);
6683 
6684  (*current_liboctave_error_handler)
6685  ("support for UMFPACK was unavailable or disabled "
6686  "when liboctave was built");
6687 #endif
6688  }
6689  else if (typ != MatrixType::Hermitian)
6690  (*current_liboctave_error_handler) ("incorrect matrix type");
6691  }
6692 
6693  return retval;
6694 }
6695 
6696 Matrix
6697 SparseMatrix::solve (MatrixType& mattype, const Matrix& b) const
6698 {
6699  octave_idx_type info;
6700  double rcond;
6701  return solve (mattype, b, info, rcond, nullptr);
6702 }
6703 
6704 Matrix
6706  octave_idx_type& info) const
6707 {
6708  double rcond;
6709  return solve (mattype, b, info, rcond, nullptr);
6710 }
6711 
6712 Matrix
6714  octave_idx_type& info, double& rcond) const
6715 {
6716  return solve (mattype, b, info, rcond, nullptr);
6717 }
6718 
6719 Matrix
6721  double& rcond, solve_singularity_handler sing_handler,
6722  bool singular_fallback) const
6723 {
6724  Matrix retval;
6725  int typ = mattype.type (false);
6726 
6727  if (typ == MatrixType::Unknown)
6728  typ = mattype.type (*this);
6729 
6730  // Only calculate the condition number for CHOLMOD/UMFPACK
6732  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6733  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6734  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6735  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6736  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6737  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6738  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6739  else if (typ == MatrixType::Tridiagonal
6741  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6742  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6743  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6744  else if (typ != MatrixType::Rectangular)
6745  (*current_liboctave_error_handler) ("unknown matrix type");
6746 
6747  // Rectangular or one of the above solvers flags a singular matrix
6748  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6749  {
6750  rcond = 1.;
6751 #if defined (USE_QRSOLVE)
6752  retval = qrsolve (*this, b, err);
6753 #else
6754  retval = dmsolve<Matrix, SparseMatrix, Matrix> (*this, b, err);
6755 #endif
6756  }
6757 
6758  return retval;
6759 }
6760 
6762 SparseMatrix::solve (MatrixType& mattype, const SparseMatrix& b) const
6763 {
6764  octave_idx_type info;
6765  double rcond;
6766  return solve (mattype, b, info, rcond, nullptr);
6767 }
6768 
6771  octave_idx_type& info) const
6772 {
6773  double rcond;
6774  return solve (mattype, b, info, rcond, nullptr);
6775 }
6776 
6779  octave_idx_type& info, double& rcond) const
6780 {
6781  return solve (mattype, b, info, rcond, nullptr);
6782 }
6783 
6786  octave_idx_type& err, double& rcond,
6787  solve_singularity_handler sing_handler,
6788  bool singular_fallback) const
6789 {
6790  SparseMatrix retval;
6791  int typ = mattype.type (false);
6792 
6793  if (typ == MatrixType::Unknown)
6794  typ = mattype.type (*this);
6795 
6797  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6798  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6799  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6800  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6801  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6802  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6803  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6804  else if (typ == MatrixType::Tridiagonal
6806  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6807  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6808  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6809  else if (typ != MatrixType::Rectangular)
6810  (*current_liboctave_error_handler) ("unknown matrix type");
6811 
6812  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6813  {
6814  rcond = 1.;
6815 #if defined (USE_QRSOLVE)
6816  retval = qrsolve (*this, b, err);
6817 #else
6819  (*this, b, err);
6820 #endif
6821  }
6822 
6823  return retval;
6824 }
6825 
6828 {
6829  octave_idx_type info;
6830  double rcond;
6831  return solve (mattype, b, info, rcond, nullptr);
6832 }
6833 
6836  octave_idx_type& info) const
6837 {
6838  double rcond;
6839  return solve (mattype, b, info, rcond, nullptr);
6840 }
6841 
6844  octave_idx_type& info, double& rcond) const
6845 {
6846  return solve (mattype, b, info, rcond, nullptr);
6847 }
6848 
6851  octave_idx_type& err, double& rcond,
6852  solve_singularity_handler sing_handler,
6853  bool singular_fallback) const
6854 {
6855  ComplexMatrix retval;
6856  int typ = mattype.type (false);
6857 
6858  if (typ == MatrixType::Unknown)
6859  typ = mattype.type (*this);
6860 
6862  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6863  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6864  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6865  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6866  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6867  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6868  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6869  else if (typ == MatrixType::Tridiagonal
6871  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6872  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6873  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6874  else if (typ != MatrixType::Rectangular)
6875  (*current_liboctave_error_handler) ("unknown matrix type");
6876 
6877  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6878  {
6879  rcond = 1.;
6880 #if defined (USE_QRSOLVE)
6881  retval = qrsolve (*this, b, err);
6882 #else
6884  (*this, b, err);
6885 #endif
6886  }
6887 
6888  return retval;
6889 }
6890 
6893 {
6894  octave_idx_type info;
6895  double rcond;
6896  return solve (mattype, b, info, rcond, nullptr);
6897 }
6898 
6901  octave_idx_type& info) const
6902 {
6903  double rcond;
6904  return solve (mattype, b, info, rcond, nullptr);
6905 }
6906 
6909  octave_idx_type& info, double& rcond) const
6910 {
6911  return solve (mattype, b, info, rcond, nullptr);
6912 }
6913 
6916  octave_idx_type& err, double& rcond,
6917  solve_singularity_handler sing_handler,
6918  bool singular_fallback) const
6919 {
6920  SparseComplexMatrix retval;
6921  int typ = mattype.type (false);
6922 
6923  if (typ == MatrixType::Unknown)
6924  typ = mattype.type (*this);
6925 
6927  retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6928  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6929  retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6930  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6931  retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6932  else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6933  retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6934  else if (typ == MatrixType::Tridiagonal
6936  retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6937  else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6938  retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6939  else if (typ != MatrixType::Rectangular)
6940  (*current_liboctave_error_handler) ("unknown matrix type");
6941 
6942  if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6943  {
6944  rcond = 1.;
6945 #if defined (USE_QRSOLVE)
6946  retval = qrsolve (*this, b, err);
6947 #else
6949  (*this, b, err);
6950 #endif
6951  }
6952 
6953  return retval;
6954 }
6955 
6957 SparseMatrix::solve (MatrixType& mattype, const ColumnVector& b) const
6958 {
6959  octave_idx_type info; double rcond;
6960  return solve (mattype, b, info, rcond);
6961 }
6962 
6965  octave_idx_type& info) const
6966 {
6967  double rcond;
6968  return solve (mattype, b, info, rcond);
6969 }
6970 
6973  octave_idx_type& info, double& rcond) const
6974 {
6975  return solve (mattype, b, info, rcond, nullptr);
6976 }
6977 
6980  octave_idx_type& info, double& rcond,
6981  solve_singularity_handler sing_handler) const
6982 {
6983  Matrix tmp (b);
6984  return solve (mattype, tmp, info, rcond,
6985  sing_handler).column (static_cast<octave_idx_type> (0));
6986 }
6987 
6990 {
6991  octave_idx_type info;
6992  double rcond;
6993  return solve (mattype, b, info, rcond, nullptr);
6994 }
6995 
6998  octave_idx_type& info) const
6999 {
7000  double rcond;
7001  return solve (mattype, b, info, rcond, nullptr);
7002 }
7003 
7006  octave_idx_type& info,
7007  double& rcond) const
7008 {
7009  return solve (mattype, b, info, rcond, nullptr);
7010 }
7011 
7014  octave_idx_type& info, double& rcond,
7015  solve_singularity_handler sing_handler) const
7016 {
7017  ComplexMatrix tmp (b);
7018  return solve (mattype, tmp, info, rcond,
7019  sing_handler).column (static_cast<octave_idx_type> (0));
7020 }
7021 
7022 Matrix
7023 SparseMatrix::solve (const Matrix& b) const
7024 {
7025  octave_idx_type info;
7026  double rcond;
7027  return solve (b, info, rcond, nullptr);
7028 }
7029 
7030 Matrix
7032 {
7033  double rcond;
7034  return solve (b, info, rcond, nullptr);
7035 }
7036 
7037 Matrix
7039  double& rcond) const
7040 {
7041  return solve (b, info, rcond, nullptr);
7042 }
7043 
7044 Matrix
7045 SparseMatrix::solve (const Matrix& b, octave_idx_type& err, double& rcond,
7046  solve_singularity_handler sing_handler) const
7047 {
7048  MatrixType mattype (*this);
7049  return solve (mattype, b, err, rcond, sing_handler);
7050 }
7051 
7054 {
7055  octave_idx_type info;
7056  double rcond;
7057  return solve (b, info, rcond, nullptr);
7058 }
7059 
7062  octave_idx_type& info) const
7063 {
7064  double rcond;
7065  return solve (b, info, rcond, nullptr);
7066 }
7067 
7070  octave_idx_type& info, double& rcond) const
7071 {
7072  return solve (b, info, rcond, nullptr);
7073 }
7074 
7076 SparseMatrix::solve (const SparseMatrix& b, octave_idx_type& err, double& rcond,
7077  solve_singularity_handler sing_handler) const
7078 {
7079  MatrixType mattype (*this);
7080  return solve (mattype, b, err, rcond, sing_handler);
7081 }
7082 
7085 {
7086  double rcond;
7087  return solve (b, info, rcond, nullptr);
7088 }
7089 
7092  double& rcond) const
7093 {
7094  return solve (b, info, rcond, nullptr);
7095 }
7096 
7099  double& rcond,
7100  solve_singularity_handler sing_handler) const
7101 {
7102  MatrixType mattype (*this);
7103  return solve (mattype, b, err, rcond, sing_handler);
7104 }
7105 
7108 {
7109  octave_idx_type info;
7110  double rcond;
7111  return solve (b, info, rcond, nullptr);
7112 }
7113 
7116 {
7117  double rcond;
7118  return solve (b, info, rcond, nullptr);
7119 }
7120 
7123  double& rcond) const
7124 {
7125  return solve (b, info, rcond, nullptr);
7126 }
7127 
7130  octave_idx_type& err, double& rcond,
7131  solve_singularity_handler sing_handler) const
7132 {
7133  MatrixType mattype (*this);
7134  return solve (mattype, b, err, rcond, sing_handler);
7135 }
7136 
7139 {
7140  octave_idx_type info; double rcond;
7141  return solve (b, info, rcond);
7142 }
7143 
7146 {
7147  double rcond;
7148  return solve (b, info, rcond);
7149 }
7150 
7153  double& rcond) const
7154 {
7155  return solve (b, info, rcond, nullptr);
7156 }
7157 
7160  double& rcond,
7161  solve_singularity_handler sing_handler) const
7162 {
7163  Matrix tmp (b);
7164  return solve (tmp, info, rcond,
7165  sing_handler).column (static_cast<octave_idx_type> (0));
7166 }
7167 
7170 {
7171  octave_idx_type info;
7172  double rcond;
7173  return solve (b, info, rcond, nullptr);
7174 }
7175 
7178 {
7179  double rcond;
7180  return solve (b, info, rcond, nullptr);
7181 }
7182 
7185  double& rcond) const
7186 {
7187  return solve (b, info, rcond, nullptr);
7188 }
7189 
7192  double& rcond,
7193  solve_singularity_handler sing_handler) const
7194 {
7195  ComplexMatrix tmp (b);
7196  return solve (tmp, info, rcond,
7197  sing_handler).column (static_cast<octave_idx_type> (0));
7198 }
7199 
7200 // other operations.
7201 
7202 bool
7204 {
7205  octave_idx_type nel = nnz ();
7206 
7207  if (neg_zero)
7208  {
7209  for (octave_idx_type i = 0; i < nel; i++)
7210  if (lo_ieee_signbit (data (i)))
7211  return true;
7212  }
7213  else
7214  {
7215  for (octave_idx_type i = 0; i < nel; i++)
7216  if (data (i) < 0)
7217  return true;
7218  }
7219 
7220  return false;
7221 }
7222 
7223 bool
7225 {
7226  octave_idx_type nel = nnz ();
7227 
7228  for (octave_idx_type i = 0; i < nel; i++)
7229  {
7230  double val = data (i);
7231  if (octave::math::isnan (val))
7232  return true;
7233  }
7234 
7235  return false;
7236 }
7237 
7238 bool
7240 {
7241  octave_idx_type nel = nnz ();
7242 
7243  for (octave_idx_type i = 0; i < nel; i++)
7244  {
7245  double val = data (i);
7246  if (octave::math::isinf (val) || octave::math::isnan (val))
7247  return true;
7248  }
7249 
7250  return false;
7251 }
7252 
7253 bool
7255 {
7256  octave_idx_type nel = nnz ();
7257 
7258  for (octave_idx_type i = 0; i < nel; i++)
7259  {
7260  double val = data (i);
7261  if (val != 0.0 && val != 1.0)
7262  return true;
7263  }
7264 
7265  return false;
7266 }
7267 
7268 bool
7270 {
7271  octave_idx_type nel = nnz ();
7272 
7273  for (octave_idx_type i = 0; i < nel; i++)
7274  if (data (i) != 0)
7275  return false;
7276 
7277  return true;
7278 }
7279 
7280 bool
7282 {
7283  octave_idx_type nel = nnz ();
7284 
7285  for (octave_idx_type i = 0; i < nel; i++)
7286  {
7287  double val = data (i);
7288  if (octave::math::isnan (val) || octave::math::x_nint (val) == val)
7289  continue;
7290  else
7291  return false;
7292  }
7293 
7294  return true;
7295 }
7296 
7297 // Return nonzero if any element of M is not an integer. Also extract
7298 // the largest and smallest values and return them in MAX_VAL and MIN_VAL.
7299 
7300 bool
7301 SparseMatrix::all_integers (double& max_val, double& min_val) const
7302 {
7303  octave_idx_type nel = nnz ();
7304 
7305  if (nel == 0)
7306  return false;
7307 
7308  max_val = data (0);
7309  min_val = data (0);
7310 
7311  for (octave_idx_type i = 0; i < nel; i++)
7312  {
7313  double val = data (i);
7314 
7315  if (val > max_val)
7316  max_val = val;
7317 
7318  if (val < min_val)
7319  min_val = val;
7320 
7321  if (octave::math::x_nint (val) != val)
7322  return false;
7323  }
7324 
7325  return true;
7326 }
7327 
7328 bool
7330 {
7332 }
7333 
7336 {
7337  if (any_element_is_nan ())
7339 
7340  octave_idx_type nr = rows ();
7341  octave_idx_type nc = cols ();
7342  octave_idx_type nz1 = nnz ();
7343  octave_idx_type nz2 = nr*nc - nz1;
7344 
7345  SparseBoolMatrix r (nr, nc, nz2);
7346 
7347  octave_idx_type ii = 0;
7348  octave_idx_type jj = 0;
7349  r.cidx (0) = 0;
7350  for (octave_idx_type i = 0; i < nc; i++)
7351  {
7352  for (octave_idx_type j = 0; j < nr; j++)
7353  {
7354  if (jj < cidx (i+1) && ridx (jj) == j)
7355  jj++;
7356  else
7357  {
7358  r.data (ii) = true;
7359  r.ridx (ii++) = j;
7360  }
7361  }
7362  r.cidx (i+1) = ii;
7363  }
7364 
7365  return r;
7366 }
7367 
7368 // FIXME: Do these really belong here? Maybe they should be in a base class?
7369 
7371 SparseMatrix::all (int dim) const
7372 {
7373  SPARSE_ALL_OP (dim);
7374 }
7375 
7377 SparseMatrix::any (int dim) const
7378 {
7379  SPARSE_ANY_OP (dim);
7380 }
7381 
7383 SparseMatrix::cumprod (int dim) const
7384 {
7385  SPARSE_CUMPROD (SparseMatrix, double, cumprod);
7386 }
7387 
7389 SparseMatrix::cumsum (int dim) const
7390 {
7391  SPARSE_CUMSUM (SparseMatrix, double, cumsum);
7392 }
7393 
7395 SparseMatrix::prod (int dim) const
7396 {
7397  if ((rows () == 1 && dim == -1) || dim == 1)
7398  return transpose ().prod (0).transpose ();
7399  else
7400  {
7401  SPARSE_REDUCTION_OP (SparseMatrix, double, *=,
7402  (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
7403  }
7404 }
7405 
7407 SparseMatrix::sum (int dim) const
7408 {
7409  SPARSE_REDUCTION_OP (SparseMatrix, double, +=, 0.0, 0.0);
7410 }
7411 
7413 SparseMatrix::sumsq (int dim) const
7414 {
7415 #define ROW_EXPR \
7416  double d = data (i); \
7417  tmp[ridx (i)] += d * d
7418 
7419 #define COL_EXPR \
7420  double d = data (i); \
7421  tmp[j] += d * d
7422 
7424  0.0, 0.0);
7425 
7426 #undef ROW_EXPR
7427 #undef COL_EXPR
7428 }
7429 
7432 {
7433  octave_idx_type nz = nnz ();
7434 
7435  SparseMatrix retval (*this);
7436 
7437  for (octave_idx_type i = 0; i < nz; i++)
7438  retval.data (i) = fabs (retval.data (i));
7439 
7440  return retval;
7441 }
7442 
7445 {
7446  return MSparse<double>::diag (k);
7447 }
7448 
7449 Matrix
7451 {
7452  return Sparse<double>::array_value ();
7453 }
7454 
7455 std::ostream&
7456 operator << (std::ostream& os, const SparseMatrix& a)
7457 {
7458  octave_idx_type nc = a.cols ();
7459 
7460  // add one to the printed indices to go from
7461  // zero-based to one-based arrays
7462  for (octave_idx_type j = 0; j < nc; j++)
7463  {
7464  octave_quit ();
7465  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7466  {
7467  os << a.ridx (i) + 1 << ' ' << j + 1 << ' ';
7468  octave::write_value<double> (os, a.data (i));
7469  os << "\n";
7470  }
7471  }
7472 
7473  return os;
7474 }
7475 
7476 std::istream&
7477 operator >> (std::istream& is, SparseMatrix& a)
7478 {
7479  typedef SparseMatrix::element_type elt_type;
7480 
7481  return read_sparse_matrix<elt_type> (is, a, octave::read_value<double>);
7482 }
7483 
7486 {
7487  return MSparse<double>::squeeze ();
7488 }
7489 
7491 SparseMatrix::reshape (const dim_vector& new_dims) const
7492 {
7493  return MSparse<double>::reshape (new_dims);
7494 }
7495 
7497 SparseMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const
7498 {
7499  return MSparse<double>::permute (vec, inv);
7500 }
7501 
7504 {
7505  return MSparse<double>::ipermute (vec);
7506 }
7507 
7508 // matrix by matrix -> matrix operations
7509 
7512 {
7513  SPARSE_SPARSE_MUL (SparseMatrix, double, double);
7514 }
7515 
7516 Matrix
7518 {
7519  FULL_SPARSE_MUL (Matrix, double);
7520 }
7521 
7522 Matrix
7523 mul_trans (const Matrix& m, const SparseMatrix& a)
7524 {
7525  FULL_SPARSE_MUL_TRANS (Matrix, double, );
7526 }
7527 
7528 Matrix
7530 {
7531  SPARSE_FULL_MUL (Matrix, double);
7532 }
7533 
7534 Matrix
7535 trans_mul (const SparseMatrix& m, const Matrix& a)
7536 {
7537  SPARSE_FULL_TRANS_MUL (Matrix, double, );
7538 }
7539 
7540 // diag * sparse and sparse * diag
7541 
7544 {
7545  return do_mul_dm_sm<SparseMatrix> (d, a);
7546 }
7547 
7550 {
7551  return do_mul_sm_dm<SparseMatrix> (a, d);
7552 }
7553 
7556 {
7557  return do_add_dm_sm<SparseMatrix> (d, a);
7558 }
7559 
7562 {
7563  return do_sub_dm_sm<SparseMatrix> (d, a);
7564 }
7565 
7568 {
7569  return do_add_sm_dm<SparseMatrix> (a, d);
7570 }
7571 
7574 {
7575  return do_sub_sm_dm<SparseMatrix> (a, d);
7576 }
7577 
7578 // perm * sparse and sparse * perm
7579 
7582 {
7583  return octinternal_do_mul_pm_sm (p, a);
7584 }
7585 
7588 {
7589  return octinternal_do_mul_sm_pm (a, p);
7590 }
7591 
7592 // FIXME: it would be nice to share code among the min/max functions below.
7593 
7594 #define EMPTY_RETURN_CHECK(T) \
7595  if (nr == 0 || nc == 0) \
7596  return T (nr, nc);
7597 
7599 min (double d, const SparseMatrix& m)
7600 {
7601  SparseMatrix result;
7602 
7603  octave_idx_type nr = m.rows ();
7604  octave_idx_type nc = m.columns ();
7605 
7607 
7608  // Count the number of nonzero elements
7609  if (d < 0.)
7610  {
7611  result = SparseMatrix (nr, nc, d);
7612  for (octave_idx_type j = 0; j < nc; j++)
7613  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7614  {
7615  double tmp = octave::math::min (d, m.data (i));
7616  if (tmp != 0.)
7617  {
7618  octave_idx_type idx = m.ridx (i) + j * nr;
7619  result.xdata (idx) = tmp;
7620  result.xridx (idx) = m.ridx (i);
7621  }
7622  }
7623  }
7624  else
7625  {
7626  octave_idx_type nel = 0;
7627  for (octave_idx_type j = 0; j < nc; j++)
7628  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7629  if (octave::math::min (d, m.data (i)) != 0.)
7630  nel++;
7631 
7632  result = SparseMatrix (nr, nc, nel);
7633 
7634  octave_idx_type ii = 0;
7635  result.xcidx (0) = 0;
7636  for (octave_idx_type j = 0; j < nc; j++)
7637  {
7638  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7639  {
7640  double tmp = octave::math::min (d, m.data (i));
7641 
7642  if (tmp != 0.)
7643  {
7644  result.xdata (ii) = tmp;
7645  result.xridx (ii++) = m.ridx (i);
7646  }
7647  }
7648  result.xcidx (j+1) = ii;
7649  }
7650  }
7651 
7652  return result;
7653 }
7654 
7656 min (const SparseMatrix& m, double d)
7657 {
7658  return min (d, m);
7659 }
7660 
7662 min (const SparseMatrix& a, const SparseMatrix& b)
7663 {
7664  SparseMatrix r;
7665 
7666  octave_idx_type a_nr = a.rows ();
7667  octave_idx_type a_nc = a.cols ();
7668  octave_idx_type b_nr = b.rows ();
7669  octave_idx_type b_nc = b.cols ();
7670 
7671  if (a_nr == b_nr && a_nc == b_nc)
7672  {
7673  r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7674 
7675  octave_idx_type jx = 0;
7676  r.cidx (0) = 0;
7677  for (octave_idx_type i = 0 ; i < a_nc ; i++)
7678  {
7679  octave_idx_type ja = a.cidx (i);
7680  octave_idx_type ja_max = a.cidx (i+1);
7681  bool ja_lt_max = ja < ja_max;
7682 
7683  octave_idx_type jb = b.cidx (i);
7684  octave_idx_type jb_max = b.cidx (i+1);
7685  bool jb_lt_max = jb < jb_max;
7686 
7687  while (ja_lt_max || jb_lt_max)
7688  {
7689  octave_quit ();
7690  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7691  {
7692  double tmp = octave::math::min (a.data (ja), 0.);
7693  if (tmp != 0.)
7694  {
7695  r.ridx (jx) = a.ridx (ja);
7696  r.data (jx) = tmp;
7697  jx++;
7698  }
7699  ja++;
7700  ja_lt_max= ja < ja_max;
7701  }
7702  else if ((! ja_lt_max)
7703  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
7704  {
7705  double tmp = octave::math::min (0., b.data (jb));
7706  if (tmp != 0.)
7707  {
7708  r.ridx (jx) = b.ridx (jb);
7709  r.data (jx) = tmp;
7710  jx++;
7711  }
7712  jb++;
7713  jb_lt_max= jb < jb_max;
7714  }
7715  else
7716  {
7717  double tmp = octave::math::min (a.data (ja), b.data (jb));
7718  if (tmp != 0.)
7719  {
7720  r.data (jx) = tmp;
7721  r.ridx (jx) = a.ridx (ja);
7722  jx++;
7723  }
7724  ja++;
7725  ja_lt_max= ja < ja_max;
7726  jb++;
7727  jb_lt_max= jb < jb_max;
7728  }
7729  }
7730  r.cidx (i+1) = jx;
7731  }
7732 
7733  r.maybe_compress ();
7734  }
7735  else
7736  {
7737  if (a_nr == 0 || a_nc == 0)
7738  r.resize (a_nr, a_nc);
7739  else if (b_nr == 0 || b_nc == 0)
7740  r.resize (b_nr, b_nc);
7741  else
7742  octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
7743  }
7744 
7745  return r;
7746 }
7747 
7749 max (double d, const SparseMatrix& m)
7750 {
7751  SparseMatrix result;
7752 
7753  octave_idx_type nr = m.rows ();
7754  octave_idx_type nc = m.columns ();
7755 
7757 
7758  // Count the number of nonzero elements
7759  if (d > 0.)
7760  {
7761  result = SparseMatrix (nr, nc, d);
7762  for (octave_idx_type j = 0; j < nc; j++)
7763  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7764  {
7765  double tmp = octave::math::max (d, m.data (i));
7766 
7767  if (tmp != 0.)
7768  {
7769  octave_idx_type idx = m.ridx (i) + j * nr;
7770  result.xdata (idx) = tmp;
7771  result.xridx (idx) = m.ridx (i);
7772  }
7773  }
7774  }
7775  else
7776  {
7777  octave_idx_type nel = 0;
7778  for (octave_idx_type j = 0; j < nc; j++)
7779  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7780  if (octave::math::max (d, m.data (i)) != 0.)
7781  nel++;
7782 
7783  result = SparseMatrix (nr, nc, nel);
7784 
7785  octave_idx_type ii = 0;
7786  result.xcidx (0) = 0;
7787  for (octave_idx_type j = 0; j < nc; j++)
7788  {
7789  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7790  {
7791  double tmp = octave::math::max (d, m.data (i));
7792  if (tmp != 0.)
7793  {
7794  result.xdata (ii) = tmp;
7795  result.xridx (ii++) = m.ridx (i);
7796  }
7797  }
7798  result.xcidx (j+1) = ii;
7799  }
7800  }
7801 
7802  return result;
7803 }
7804 
7806 max (const SparseMatrix& m, double d)
7807 {
7808  return max (d, m);
7809 }
7810 
7812 max (const SparseMatrix& a, const SparseMatrix& b)
7813 {
7814  SparseMatrix r;
7815 
7816  octave_idx_type a_nr = a.rows ();
7817  octave_idx_type a_nc = a.cols ();
7818  octave_idx_type b_nr = b.rows ();
7819  octave_idx_type b_nc = b.cols ();
7820 
7821  if (a_nr == b_nr && a_nc == b_nc)
7822  {
7823  r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7824 
7825  octave_idx_type jx = 0;
7826  r.cidx (0) = 0;
7827  for (octave_idx_type i = 0 ; i < a_nc ; i++)
7828  {
7829  octave_idx_type ja = a.cidx (i);
7830  octave_idx_type ja_max = a.cidx (i+1);
7831  bool ja_lt_max = ja < ja_max;
7832 
7833  octave_idx_type jb = b.cidx (i);
7834  octave_idx_type jb_max = b.cidx (i+1);
7835  bool jb_lt_max = jb < jb_max;
7836 
7837  while (ja_lt_max || jb_lt_max)
7838  {
7839  octave_quit ();
7840  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7841  {
7842  double tmp = octave::math::max (a.data (ja), 0.);
7843  if (tmp != 0.)
7844  {
7845  r.ridx (jx) = a.ridx (ja);
7846  r.data (jx) = tmp;
7847  jx++;
7848  }
7849  ja++;
7850  ja_lt_max= ja < ja_max;
7851  }
7852  else if ((! ja_lt_max)
7853  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
7854  {
7855  double tmp = octave::math::max (0., b.data (jb));
7856  if (tmp != 0.)
7857  {
7858  r.ridx (jx) = b.ridx (jb);
7859  r.data (jx) = tmp;
7860  jx++;
7861  }
7862  jb++;
7863  jb_lt_max= jb < jb_max;
7864  }
7865  else
7866  {
7867  double tmp = octave::math::max (a.data (ja), b.data (jb));
7868  if (tmp != 0.)
7869  {
7870  r.data (jx) = tmp;
7871  r.ridx (jx) = a.ridx (ja);
7872  jx++;
7873  }
7874  ja++;
7875  ja_lt_max= ja < ja_max;
7876  jb++;
7877  jb_lt_max= jb < jb_max;
7878  }
7879  }
7880  r.cidx (i+1) = jx;
7881  }
7882 
7883  r.maybe_compress ();
7884  }
7885  else
7886  {
7887  if (a_nr == 0 || a_nc == 0)
7888  r.resize (a_nr, a_nc);
7889  else if (b_nr == 0 || b_nc == 0)
7890  r.resize (b_nr, b_nc);
7891  else
7892  octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc);
7893  }
7894 
7895  return r;
7896 }
7897 
7900 
7903 
base_det< double > DET
Definition: DET.h:86
#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)
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
T * fortran_vec()
Size of the specified dimension.
Definition: Array-base.cc:1764
octave_idx_type rows() const
Definition: Array.h:459
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array-base.cc:1023
const T * data() const
Size of the specified dimension.
Definition: Array.h:663
octave_idx_type cols() const
Definition: Array.h:469
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() const
Definition: DiagArray2.h:95
octave_idx_type cols() 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() 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() const
Definition: MatrixType.h:122
void mark_as_unsymmetric()
Definition: MatrixType.cc:920
@ 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_idx_type * triangular_perm() const
Definition: MatrixType.h:136
int nupper() const
Definition: MatrixType.h:101
int nlower() const
Definition: MatrixType.h:103
void info() const
Definition: MatrixType.cc:846
MatrixType transpose() const
Definition: MatrixType.cc:973
int type(bool quiet=true)
Definition: MatrixType.cc:656
void mark_as_rectangular()
Definition: MatrixType.h:155
bool is_dense() const
Definition: MatrixType.h:105
Definition: dMatrix.h:42
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:422
SparseComplexMatrix & insert(const SparseComplexMatrix &a, octave_idx_type r, octave_idx_type c)
Definition: CSparse.cc:556
SparseMatrix & insert(const SparseMatrix &a, octave_idx_type r, octave_idx_type c)
Definition: dSparse.cc:171
SparseMatrix squeeze() const
Definition: dSparse.cc:7485
bool issymmetric() const
Definition: dSparse.cc:131
bool any_element_not_one_or_zero() const
Definition: dSparse.cc:7254
SparseMatrix reshape(const dim_vector &new_dims) const
Definition: dSparse.cc:7491
SparseMatrix sum(int dim=-1) const
Definition: dSparse.cc:7407
SparseBoolMatrix any(int dim=-1) const
Definition: dSparse.cc:7377
bool any_element_is_inf_or_nan() const
Definition: dSparse.cc:7239
bool all_integers(double &max_val, double &min_val) const
Definition: dSparse.cc:7301
bool any_element_is_nan() const
Definition: dSparse.cc:7224
bool too_large_for_float() const
Definition: dSparse.cc:7329
bool operator==(const SparseMatrix &a) const
Definition: dSparse.cc:101
SparseMatrix transpose() const
Definition: dSparse.h:126
DET determinant() const
Definition: dSparse.cc:1012
SparseMatrix abs() const
Definition: dSparse.cc:7431
SparseMatrix inverse() const
Definition: dSparse.cc:603
bool operator!=(const SparseMatrix &a) const
Definition: dSparse.cc:125
RowVector row(octave_idx_type i) const
Definition: dSparse.cc:503
SparseMatrix cumprod(int dim=-1) const
Definition: dSparse.cc:7383
SparseMatrix diag(octave_idx_type k=0) const
Definition: dSparse.cc:7444
SparseMatrix prod(int dim=-1) const
Definition: dSparse.cc:7395
bool all_elements_are_zero() const
Definition: dSparse.cc:7269
SparseBoolMatrix all(int dim=-1) const
Definition: dSparse.cc:7371
SparseMatrix min(int dim=-1) const
Definition: dSparse.cc:337
SparseMatrix permute(const Array< octave_idx_type > &vec, bool inv=false) const
Definition: dSparse.cc:7497
SparseMatrix cumsum(int dim=-1) const
Definition: dSparse.cc:7389
ColumnVector column(octave_idx_type i) const
Definition: dSparse.cc:522
SparseMatrix concat(const SparseMatrix &rb, const Array< octave_idx_type > &ra_idx)
Definition: dSparse.cc:534
SparseMatrix ipermute(const Array< octave_idx_type > &vec) const
Definition: dSparse.cc:7503
SparseBoolMatrix operator!() const
Definition: dSparse.cc:7335
SparseMatrix sumsq(int dim=-1) const
Definition: dSparse.cc:7413
SparseMatrix()
Definition: dSparse.h:53
bool all_elements_are_int_or_inf_or_nan() const
Definition: dSparse.cc:7281
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6697
SparseMatrix max(int dim=-1) const
Definition: dSparse.cc:186
bool any_element_is_negative(bool=false) const
Definition: dSparse.cc:7203
Matrix matrix_value() const
Definition: dSparse.cc:7450
octave_idx_type cols() const
Definition: Sparse.h:352
octave_idx_type * xcidx()
Definition: Sparse.h:602
T * xdata()
Definition: Sparse.h:576
Array< T > array_value() const
Definition: Sparse.cc:2767
octave_idx_type columns() const
Definition: Sparse.h:353
octave_idx_type * cidx()
Definition: Sparse.h:596
T * data()
Definition: Sparse.h:574
bool test_any(F fcn) const
Definition: Sparse.h:663
octave_idx_type * ridx()
Definition: Sparse.h:583
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 nnz() const
Actual number of nonzero terms.
Definition: Sparse.h:339
octave_idx_type rows() const
Definition: Sparse.h:351
void change_capacity(octave_idx_type nz)
Definition: Sparse.h:556
dim_vector dims() const
Definition: Sparse.h:371
octave_idx_type * xridx()
Definition: Sparse.h:589
Definition: DET.h:39
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
octave_idx_type ndims() const
Number of dimensions.
Definition: dim-vector.h:257
int first_non_singleton(int def=0) const
Definition: dim-vector.h:444
Matrix trans_mul(const SparseMatrix &m, const Matrix &a)
Definition: dSparse.cc:7535
SparseMatrix real(const SparseComplexMatrix &a)
Definition: dSparse.cc:554
#define EMPTY_RETURN_CHECK(T)
Definition: dSparse.cc:7594
SparseMatrix min(double d, const SparseMatrix &m)
Definition: dSparse.cc:7599
Matrix mul_trans(const Matrix &m, const SparseMatrix &a)
Definition: dSparse.cc:7523
#define COL_EXPR
#define ROW_EXPR
SparseMatrix operator-(const DiagMatrix &d, const SparseMatrix &a)
Definition: dSparse.cc:7561
SparseMatrix max(double d, const SparseMatrix &m)
Definition: dSparse.cc:7749
std::istream & operator>>(std::istream &is, SparseMatrix &a)
Definition: dSparse.cc:7477
SparseMatrix operator*(const SparseMatrix &m, const SparseMatrix &a)
Definition: dSparse.cc:7511
SparseMatrix imag(const SparseComplexMatrix &a)
Definition: dSparse.cc:575
SparseMatrix operator+(const DiagMatrix &d, const SparseMatrix &a)
Definition: dSparse.cc:7555
std::ostream & operator<<(std::ostream &os, const SparseMatrix &a)
Definition: dSparse.cc:7456
#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 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:132
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
F77_RET_T const F77_INT F77_CMPLX * A
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:56
T octave_idx_type m
Definition: mx-inlines.cc:781
octave_idx_type n
Definition: mx-inlines.cc:761
T * r
Definition: mx-inlines.cc:781
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
octave_idx_type from_suitesparse_long(SuiteSparse_long x)
Definition: oct-sparse.h:200
#define UMFPACK_DNAME(name)
Definition: oct-sparse.h:169
octave_idx_type from_size_t(std::size_t x)
Definition: oct-sparse.h:211
const octave_base_value const Array< octave_idx_type > & ra_idx
template SparseMatrix dmsolve< SparseMatrix, SparseMatrix, SparseMatrix >(const SparseMatrix &, const SparseMatrix &, octave_idx_type &)
template Matrix dmsolve< Matrix, SparseMatrix, Matrix >(const SparseMatrix &, const Matrix &, octave_idx_type &)
template ComplexMatrix dmsolve< ComplexMatrix, SparseMatrix, ComplexMatrix >(const SparseMatrix &, const ComplexMatrix &, octave_idx_type &)
template SparseComplexMatrix dmsolve< SparseComplexMatrix, SparseMatrix, SparseComplexMatrix >(const SparseMatrix &, const SparseComplexMatrix &, octave_idx_type &)
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition: sparse-qr.cc:3270
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