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;
6