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)