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