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