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