GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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