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
CSparse.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 
34 #include "quit.h"
35 #include "lo-ieee.h"
36 #include "lo-mappers.h"
37 #include "f77-fcn.h"
38 #include "dRowVector.h"
39 #include "mx-m-cs.h"
40 #include "mx-cs-m.h"
41 #include "mx-cm-s.h"
42 #include "mx-fcm-fs.h"
43 #include "mx-s-cm.h"
44 #include "mx-fs-fcm.h"
45 #include "oct-locbuf.h"
46 
47 #include "dDiagMatrix.h"
48 #include "CDiagMatrix.h"
49 #include "CSparse.h"
50 #include "boolSparse.h"
51 #include "dSparse.h"
52 #include "functor.h"
53 #include "oct-spparms.h"
54 #include "SparseCmplxLU.h"
55 #include "oct-sparse.h"
56 #include "sparse-util.h"
57 #include "SparseCmplxCHOL.h"
58 #include "SparseCmplxQR.h"
59 
60 #include "Sparse-op-defs.h"
61 
62 #include "Sparse-diag-op-defs.h"
63 
64 #include "Sparse-perm-op-defs.h"
65 
66 // Define whether to use a basic QR solver or one that uses a Dulmange
67 // Mendelsohn factorization to seperate the problem into under-determined,
68 // well-determined and over-determined parts and solves them seperately
69 #ifndef USE_QRSOLVE
70 #include "sparse-dmsolve.cc"
71 #endif
72 
73 // Fortran functions we call.
74 extern "C"
75 {
76  F77_RET_T
77  F77_FUNC (zgbtrf, ZGBTRF) (const octave_idx_type&, const octave_idx_type&,
78  const octave_idx_type&, const octave_idx_type&,
79  Complex*, const octave_idx_type&,
80  octave_idx_type*, octave_idx_type&);
81 
82  F77_RET_T
83  F77_FUNC (zgbtrs, ZGBTRS) (F77_CONST_CHAR_ARG_DECL,
84  const octave_idx_type&, const octave_idx_type&,
85  const octave_idx_type&, const octave_idx_type&,
86  const Complex*, const octave_idx_type&,
87  const octave_idx_type*, Complex*,
88  const octave_idx_type&, octave_idx_type&
90 
91  F77_RET_T
92  F77_FUNC (zgbcon, ZGBCON) (F77_CONST_CHAR_ARG_DECL,
93  const octave_idx_type&, const octave_idx_type&,
94  const octave_idx_type&, Complex*,
95  const octave_idx_type&, const octave_idx_type*,
96  const double&, double&, Complex*, double*,
97  octave_idx_type&
99 
100  F77_RET_T
101  F77_FUNC (zpbtrf, ZPBTRF) (F77_CONST_CHAR_ARG_DECL,
102  const octave_idx_type&, const octave_idx_type&,
103  Complex*, const octave_idx_type&, octave_idx_type&
105 
106  F77_RET_T
107  F77_FUNC (zpbtrs, ZPBTRS) (F77_CONST_CHAR_ARG_DECL,
108  const octave_idx_type&, const octave_idx_type&,
109  const octave_idx_type&, Complex*,
110  const octave_idx_type&, Complex*,
111  const octave_idx_type&, octave_idx_type&
113 
114  F77_RET_T
115  F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL,
116  const octave_idx_type&, const octave_idx_type&,
117  Complex*, const octave_idx_type&, const double&,
118  double&, Complex*, double*, octave_idx_type&
120 
121  F77_RET_T
122  F77_FUNC (zgttrf, ZGTTRF) (const octave_idx_type&, Complex*, Complex*,
123  Complex*, Complex*, octave_idx_type*,
124  octave_idx_type&);
125 
126  F77_RET_T
127  F77_FUNC (zgttrs, ZGTTRS) (F77_CONST_CHAR_ARG_DECL,
128  const octave_idx_type&, const octave_idx_type&,
129  const Complex*, const Complex*, const Complex*,
130  const Complex*, const octave_idx_type*,
131  Complex *, const octave_idx_type&,
132  octave_idx_type&
134 
135  F77_RET_T
136  F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&,
137  double*, Complex*, Complex*,
138  const octave_idx_type&, octave_idx_type&);
139 
140  F77_RET_T
141  F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&,
142  Complex*, Complex*, Complex*, Complex*,
143  const octave_idx_type&, octave_idx_type&);
144 }
145 
147  : MSparse<Complex> (a)
148 {
149 }
150 
152  : MSparse<Complex> (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) = Complex (a.data (i));
163  ridx (i) = a.ridx (i);
164  }
165 }
166 
168  : MSparse<Complex> (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 bool
187 {
188  octave_idx_type nr = rows ();
189  octave_idx_type nc = cols ();
190  octave_idx_type nz = nnz ();
191  octave_idx_type nr_a = a.rows ();
192  octave_idx_type nc_a = a.cols ();
193  octave_idx_type nz_a = a.nnz ();
194 
195  if (nr != nr_a || nc != nc_a || nz != nz_a)
196  return false;
197 
198  for (octave_idx_type i = 0; i < nc + 1; i++)
199  if (cidx (i) != a.cidx (i))
200  return false;
201 
202  for (octave_idx_type i = 0; i < nz; i++)
203  if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
204  return false;
205 
206  return true;
207 }
208 
209 bool
211 {
212  return !(*this == a);
213 }
214 
215 bool
217 {
218  octave_idx_type nr = rows ();
219  octave_idx_type nc = cols ();
220 
221  if (nr == nc && nr > 0)
222  {
223  for (octave_idx_type j = 0; j < nc; j++)
224  {
225  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
226  {
227  octave_idx_type ri = ridx (i);
228 
229  if (ri != j)
230  {
231  bool found = false;
232 
233  for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
234  {
235  if (ridx (k) == j)
236  {
237  if (data (i) == conj (data (k)))
238  found = true;
239  break;
240  }
241  }
242 
243  if (! found)
244  return false;
245  }
246  }
247  }
248 
249  return true;
250  }
251 
252  return false;
253 }
254 
256 
259 {
260  Array<octave_idx_type> dummy_idx;
261  return max (dummy_idx, dim);
262 }
263 
266 {
267  SparseComplexMatrix result;
268  dim_vector dv = dims ();
269  octave_idx_type nr = dv(0);
270  octave_idx_type nc = dv(1);
271 
272 
273  if (dim >= dv.length ())
274  {
275  idx_arg.resize (dim_vector (nr, nc), 0);
276  return *this;
277  }
278 
279  if (dim < 0)
280  dim = dv.first_non_singleton ();
281 
282  if (dim == 0)
283  {
284  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
285 
286  if (nr == 0 || nc == 0 || dim >= dv.length ())
287  return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
288 
289  octave_idx_type nel = 0;
290  for (octave_idx_type j = 0; j < nc; j++)
291  {
292  Complex tmp_max;
293  double abs_max = octave_NaN;
294  octave_idx_type idx_j = 0;
295  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
296  {
297  if (ridx (i) != idx_j)
298  break;
299  else
300  idx_j++;
301  }
302 
303  if (idx_j != nr)
304  {
305  tmp_max = 0.;
306  abs_max = 0.;
307  }
308 
309  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
310  {
311  Complex tmp = data (i);
312 
313  if (xisnan (tmp))
314  continue;
315 
316  double abs_tmp = std::abs (tmp);
317 
318  if (xisnan (abs_max) || abs_tmp > abs_max)
319  {
320  idx_j = ridx (i);
321  tmp_max = tmp;
322  abs_max = abs_tmp;
323  }
324  }
325 
326  idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_j;
327  if (abs_max != 0.)
328  nel++;
329  }
330 
331  result = SparseComplexMatrix (1, nc, nel);
332 
333  octave_idx_type ii = 0;
334  result.xcidx (0) = 0;
335  for (octave_idx_type j = 0; j < nc; j++)
336  {
337  Complex tmp = elem (idx_arg(j), j);
338  if (tmp != 0.)
339  {
340  result.xdata (ii) = tmp;
341  result.xridx (ii++) = 0;
342  }
343  result.xcidx (j+1) = ii;
344  }
345  }
346  else
347  {
348  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
349 
350  if (nr == 0 || nc == 0 || dim >= dv.length ())
351  return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
352 
354 
355  for (octave_idx_type i = 0; i < nr; i++)
356  found[i] = 0;
357 
358  for (octave_idx_type j = 0; j < nc; j++)
359  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
360  if (found[ridx (i)] == -j)
361  found[ridx (i)] = -j - 1;
362 
363  for (octave_idx_type i = 0; i < nr; i++)
364  if (found[i] > -nc && found[i] < 0)
365  idx_arg.elem (i) = -found[i];
366 
367  for (octave_idx_type j = 0; j < nc; j++)
368  {
369  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
370  {
371  octave_idx_type ir = ridx (i);
372  octave_idx_type ix = idx_arg.elem (ir);
373  Complex tmp = data (i);
374 
375  if (xisnan (tmp))
376  continue;
377  else if (ix == -1 || std::abs (tmp) > std::abs (elem (ir, ix)))
378  idx_arg.elem (ir) = j;
379  }
380  }
381 
382  octave_idx_type nel = 0;
383  for (octave_idx_type j = 0; j < nr; j++)
384  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
385  nel++;
386 
387  result = SparseComplexMatrix (nr, 1, nel);
388 
389  octave_idx_type ii = 0;
390  result.xcidx (0) = 0;
391  result.xcidx (1) = nel;
392  for (octave_idx_type j = 0; j < nr; j++)
393  {
394  if (idx_arg(j) == -1)
395  {
396  idx_arg(j) = 0;
397  result.xdata (ii) = Complex_NaN_result;
398  result.xridx (ii++) = j;
399  }
400  else
401  {
402  Complex tmp = elem (j, idx_arg(j));
403  if (tmp != 0.)
404  {
405  result.xdata (ii) = tmp;
406  result.xridx (ii++) = j;
407  }
408  }
409  }
410  }
411 
412  return result;
413 }
414 
417 {
418  Array<octave_idx_type> dummy_idx;
419  return min (dummy_idx, dim);
420 }
421 
424 {
425  SparseComplexMatrix result;
426  dim_vector dv = dims ();
427  octave_idx_type nr = dv(0);
428  octave_idx_type nc = dv(1);
429 
430  if (dim >= dv.length ())
431  {
432  idx_arg.resize (dim_vector (nr, nc), 0);
433  return *this;
434  }
435 
436  if (dim < 0)
437  dim = dv.first_non_singleton ();
438 
439  if (dim == 0)
440  {
441  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
442 
443  if (nr == 0 || nc == 0 || dim >= dv.length ())
444  return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
445 
446  octave_idx_type nel = 0;
447  for (octave_idx_type j = 0; j < nc; j++)
448  {
449  Complex tmp_min;
450  double abs_min = octave_NaN;
451  octave_idx_type idx_j = 0;
452  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
453  {
454  if (ridx (i) != idx_j)
455  break;
456  else
457  idx_j++;
458  }
459 
460  if (idx_j != nr)
461  {
462  tmp_min = 0.;
463  abs_min = 0.;
464  }
465 
466  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
467  {
468  Complex tmp = data (i);
469 
470  if (xisnan (tmp))
471  continue;
472 
473  double abs_tmp = std::abs (tmp);
474 
475  if (xisnan (abs_min) || abs_tmp < abs_min)
476  {
477  idx_j = ridx (i);
478  tmp_min = tmp;
479  abs_min = abs_tmp;
480  }
481  }
482 
483  idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_j;
484  if (abs_min != 0.)
485  nel++;
486  }
487 
488  result = SparseComplexMatrix (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  Complex 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  else
504  {
505  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
506 
507  if (nr == 0 || nc == 0 || dim >= dv.length ())
508  return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
509 
511 
512  for (octave_idx_type i = 0; i < nr; i++)
513  found[i] = 0;
514 
515  for (octave_idx_type j = 0; j < nc; j++)
516  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
517  if (found[ridx (i)] == -j)
518  found[ridx (i)] = -j - 1;
519 
520  for (octave_idx_type i = 0; i < nr; i++)
521  if (found[i] > -nc && found[i] < 0)
522  idx_arg.elem (i) = -found[i];
523 
524  for (octave_idx_type j = 0; j < nc; j++)
525  {
526  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
527  {
528  octave_idx_type ir = ridx (i);
529  octave_idx_type ix = idx_arg.elem (ir);
530  Complex tmp = data (i);
531 
532  if (xisnan (tmp))
533  continue;
534  else if (ix == -1 || std::abs (tmp) < std::abs (elem (ir, ix)))
535  idx_arg.elem (ir) = j;
536  }
537  }
538 
539  octave_idx_type nel = 0;
540  for (octave_idx_type j = 0; j < nr; j++)
541  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
542  nel++;
543 
544  result = SparseComplexMatrix (nr, 1, nel);
545 
546  octave_idx_type ii = 0;
547  result.xcidx (0) = 0;
548  result.xcidx (1) = nel;
549  for (octave_idx_type j = 0; j < nr; j++)
550  {
551  if (idx_arg(j) == -1)
552  {
553  idx_arg(j) = 0;
554  result.xdata (ii) = Complex_NaN_result;
555  result.xridx (ii++) = j;
556  }
557  else
558  {
559  Complex tmp = elem (j, idx_arg(j));
560  if (tmp != 0.)
561  {
562  result.xdata (ii) = tmp;
563  result.xridx (ii++) = j;
564  }
565  }
566  }
567  }
568 
569  return result;
570 }
571 
572 /*
573 
574 %!assert (max (max (speye (65536) * 1i)), sparse (1i))
575 %!assert (min (min (speye (65536) * 1i)), sparse (0))
576 %!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
577 %!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
578 %!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
579 %!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
580 %!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
581 %!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
582 %!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
583 %!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
584 
585 */
586 
589 {
590  octave_idx_type nc = columns ();
591  ComplexRowVector retval (nc, 0);
592 
593  for (octave_idx_type j = 0; j < nc; j++)
594  for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
595  {
596  if (ridx (k) == i)
597  {
598  retval(j) = data (k);
599  break;
600  }
601  }
602 
603  return retval;
604 }
605 
608 {
609  octave_idx_type nr = rows ();
610  ComplexColumnVector retval (nr, 0);
611 
612  for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
613  retval(ridx (k)) = data (k);
614 
615  return retval;
616 }
617 
618 // destructive insert/delete/reorder operations
619 
623 {
624  SparseComplexMatrix tmp (a);
625  return insert (tmp /*a*/, r, c);
626 }
627 
631 {
632  MSparse<Complex>::insert (a, r, c);
633  return *this;
634 }
635 
638  const Array<octave_idx_type>& indx)
639 {
640  SparseComplexMatrix tmp (a);
641  return insert (tmp /*a*/, indx);
642 }
643 
646  const Array<octave_idx_type>& indx)
647 {
648  MSparse<Complex>::insert (a, indx);
649  return *this;
650 }
651 
655 {
656  // Don't use numel to avoid all possiblity of an overflow
657  if (rb.rows () > 0 && rb.cols () > 0)
658  insert (rb, ra_idx(0), ra_idx(1));
659  return *this;
660 }
661 
665 {
666  SparseComplexMatrix tmp (rb);
667  if (rb.rows () > 0 && rb.cols () > 0)
668  insert (tmp, ra_idx(0), ra_idx(1));
669  return *this;
670 }
671 
674 {
676 }
677 
680 {
681  octave_idx_type nr = rows ();
682  octave_idx_type nc = cols ();
683  octave_idx_type nz = nnz ();
684  SparseComplexMatrix retval (nc, nr, nz);
685 
686  for (octave_idx_type i = 0; i < nz; i++)
687  retval.xcidx (ridx (i) + 1)++;
688  // retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
689  nz = 0;
690  for (octave_idx_type i = 1; i <= nr; i++)
691  {
692  const octave_idx_type tmp = retval.xcidx (i);
693  retval.xcidx (i) = nz;
694  nz += tmp;
695  }
696  // retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
697 
698  for (octave_idx_type j = 0; j < nc; j++)
699  for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
700  {
701  octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
702  retval.xridx (q) = j;
703  retval.xdata (q) = conj (data (k));
704  }
705  assert (nnz () == retval.xcidx (nr));
706  // retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
707  // and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
708 
709  return retval;
710 }
711 
714 {
715  octave_idx_type nr = a.rows ();
716  octave_idx_type nc = a.cols ();
717  octave_idx_type nz = a.nnz ();
718  SparseComplexMatrix retval (nc, nr, nz);
719 
720  for (octave_idx_type i = 0; i < nc + 1; i++)
721  retval.cidx (i) = a.cidx (i);
722 
723  for (octave_idx_type i = 0; i < nz; i++)
724  {
725  retval.data (i) = conj (a.data (i));
726  retval.ridx (i) = a.ridx (i);
727  }
728 
729  return retval;
730 }
731 
734 {
735  octave_idx_type info;
736  double rcond;
737  MatrixType mattype (*this);
738  return inverse (mattype, info, rcond, 0, 0);
739 }
740 
743 {
744  octave_idx_type info;
745  double rcond;
746  return inverse (mattype, info, rcond, 0, 0);
747 }
748 
751 {
752  double rcond;
753  return inverse (mattype, info, rcond, 0, 0);
754 }
755 
758  double& rcond, const bool,
759  const bool calccond) const
760 {
761  SparseComplexMatrix retval;
762 
763  octave_idx_type nr = rows ();
764  octave_idx_type nc = cols ();
765  info = 0;
766 
767  if (nr == 0 || nc == 0 || nr != nc)
768  (*current_liboctave_error_handler) ("inverse requires square matrix");
769  else
770  {
771  // Print spparms("spumoni") info if requested
772  int typ = mattyp.type ();
773  mattyp.info ();
774 
776  {
778  retval = transpose ();
779  else
780  retval = *this;
781 
782  // Force make_unique to be called
783  Complex *v = retval.data ();
784 
785  if (calccond)
786  {
787  double dmax = 0.;
788  double dmin = octave_Inf;
789  for (octave_idx_type i = 0; i < nr; i++)
790  {
791  double tmp = std::abs (v[i]);
792  if (tmp > dmax)
793  dmax = tmp;
794  if (tmp < dmin)
795  dmin = tmp;
796  }
797  rcond = dmin / dmax;
798  }
799 
800  for (octave_idx_type i = 0; i < nr; i++)
801  v[i] = 1.0 / v[i];
802  }
803  else
804  (*current_liboctave_error_handler) ("incorrect matrix type");
805  }
806 
807  return retval;
808 }
809 
812  double& rcond, const bool,
813  const bool calccond) const
814 {
815  SparseComplexMatrix retval;
816 
817  octave_idx_type nr = rows ();
818  octave_idx_type nc = cols ();
819  info = 0;
820 
821  if (nr == 0 || nc == 0 || nr != nc)
822  (*current_liboctave_error_handler) ("inverse requires square matrix");
823  else
824  {
825  // Print spparms("spumoni") info if requested
826  int typ = mattyp.type ();
827  mattyp.info ();
828 
829  if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper
830  || typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
831  {
832  double anorm = 0.;
833  double ainvnorm = 0.;
834 
835  if (calccond)
836  {
837  // Calculate the 1-norm of matrix for rcond calculation
838  for (octave_idx_type j = 0; j < nr; j++)
839  {
840  double atmp = 0.;
841  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
842  atmp += std::abs (data (i));
843  if (atmp > anorm)
844  anorm = atmp;
845  }
846  }
847 
848  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
849  {
850  octave_idx_type nz = nnz ();
851  octave_idx_type cx = 0;
852  octave_idx_type nz2 = nz;
853  retval = SparseComplexMatrix (nr, nc, nz2);
854 
855  for (octave_idx_type i = 0; i < nr; i++)
856  {
857  octave_quit ();
858  // place the 1 in the identity position
859  octave_idx_type cx_colstart = cx;
860 
861  if (cx == nz2)
862  {
863  nz2 *= 2;
864  retval.change_capacity (nz2);
865  }
866 
867  retval.xcidx (i) = cx;
868  retval.xridx (cx) = i;
869  retval.xdata (cx) = 1.0;
870  cx++;
871 
872  // iterate accross columns of input matrix
873  for (octave_idx_type j = i+1; j < nr; j++)
874  {
875  Complex v = 0.;
876  // iterate to calculate sum
877  octave_idx_type colXp = retval.xcidx (i);
878  octave_idx_type colUp = cidx (j);
879  octave_idx_type rpX, rpU;
880 
881  if (cidx (j) == cidx (j+1))
882  {
883  (*current_liboctave_error_handler)
884  ("division by zero");
885  goto inverse_singular;
886  }
887 
888  do
889  {
890  octave_quit ();
891  rpX = retval.xridx (colXp);
892  rpU = ridx (colUp);
893 
894  if (rpX < rpU)
895  colXp++;
896  else if (rpX > rpU)
897  colUp++;
898  else
899  {
900  v -= retval.xdata (colXp) * data (colUp);
901  colXp++;
902  colUp++;
903  }
904  }
905  while (rpX < j && rpU < j && colXp < cx && colUp < nz);
906 
907 
908  // get A(m,m)
909  if (typ == MatrixType::Upper)
910  colUp = cidx (j+1) - 1;
911  else
912  colUp = cidx (j);
913  Complex pivot = data (colUp);
914  if (pivot == 0. || ridx (colUp) != j)
915  {
916  (*current_liboctave_error_handler)
917  ("division by zero");
918  goto inverse_singular;
919  }
920 
921  if (v != 0.)
922  {
923  if (cx == nz2)
924  {
925  nz2 *= 2;
926  retval.change_capacity (nz2);
927  }
928 
929  retval.xridx (cx) = j;
930  retval.xdata (cx) = v / pivot;
931  cx++;
932  }
933  }
934 
935  // get A(m,m)
936  octave_idx_type colUp;
937  if (typ == MatrixType::Upper)
938  colUp = cidx (i+1) - 1;
939  else
940  colUp = cidx (i);
941  Complex pivot = data (colUp);
942  if (pivot == 0. || ridx (colUp) != i)
943  {
944  (*current_liboctave_error_handler) ("division by zero");
945  goto inverse_singular;
946  }
947 
948  if (pivot != 1.0)
949  for (octave_idx_type j = cx_colstart; j < cx; j++)
950  retval.xdata (j) /= pivot;
951  }
952  retval.xcidx (nr) = cx;
953  retval.maybe_compress ();
954  }
955  else
956  {
957  octave_idx_type nz = nnz ();
958  octave_idx_type cx = 0;
959  octave_idx_type nz2 = nz;
960  retval = SparseComplexMatrix (nr, nc, nz2);
961 
962  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
964 
965  octave_idx_type *perm = mattyp.triangular_perm ();
966  if (typ == MatrixType::Permuted_Upper)
967  {
968  for (octave_idx_type i = 0; i < nr; i++)
969  rperm[perm[i]] = i;
970  }
971  else
972  {
973  for (octave_idx_type i = 0; i < nr; i++)
974  rperm[i] = perm[i];
975  for (octave_idx_type i = 0; i < nr; i++)
976  perm[rperm[i]] = i;
977  }
978 
979  for (octave_idx_type i = 0; i < nr; i++)
980  {
981  octave_quit ();
982  octave_idx_type iidx = rperm[i];
983 
984  for (octave_idx_type j = 0; j < nr; j++)
985  work[j] = 0.;
986 
987  // place the 1 in the identity position
988  work[iidx] = 1.0;
989 
990  // iterate accross columns of input matrix
991  for (octave_idx_type j = iidx+1; j < nr; j++)
992  {
993  Complex v = 0.;
994  octave_idx_type jidx = perm[j];
995  // iterate to calculate sum
996  for (octave_idx_type k = cidx (jidx);
997  k < cidx (jidx+1); k++)
998  {
999  octave_quit ();
1000  v -= work[ridx (k)] * data (k);
1001  }
1002 
1003  // get A(m,m)
1004  Complex pivot;
1005  if (typ == MatrixType::Permuted_Upper)
1006  pivot = data (cidx (jidx+1) - 1);
1007  else
1008  pivot = data (cidx (jidx));
1009  if (pivot == 0.)
1010  {
1011  (*current_liboctave_error_handler)
1012  ("division by zero");
1013  goto inverse_singular;
1014  }
1015 
1016  work[j] = v / pivot;
1017  }
1018 
1019  // get A(m,m)
1020  octave_idx_type colUp;
1021  if (typ == MatrixType::Permuted_Upper)
1022  colUp = cidx (perm[iidx]+1) - 1;
1023  else
1024  colUp = cidx (perm[iidx]);
1025 
1026  Complex pivot = data (colUp);
1027  if (pivot == 0.)
1028  {
1029  (*current_liboctave_error_handler)
1030  ("division by zero");
1031  goto inverse_singular;
1032  }
1033 
1034  octave_idx_type new_cx = cx;
1035  for (octave_idx_type j = iidx; j < nr; j++)
1036  if (work[j] != 0.0)
1037  {
1038  new_cx++;
1039  if (pivot != 1.0)
1040  work[j] /= pivot;
1041  }
1042 
1043  if (cx < new_cx)
1044  {
1045  nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1046  retval.change_capacity (nz2);
1047  }
1048 
1049  retval.xcidx (i) = cx;
1050  for (octave_idx_type j = iidx; j < nr; j++)
1051  if (work[j] != 0.)
1052  {
1053  retval.xridx (cx) = j;
1054  retval.xdata (cx++) = work[j];
1055  }
1056  }
1057 
1058  retval.xcidx (nr) = cx;
1059  retval.maybe_compress ();
1060  }
1061 
1062  if (calccond)
1063  {
1064  // Calculate the 1-norm of inverse matrix for rcond calculation
1065  for (octave_idx_type j = 0; j < nr; j++)
1066  {
1067  double atmp = 0.;
1068  for (octave_idx_type i = retval.cidx (j);
1069  i < retval.cidx (j+1); i++)
1070  atmp += std::abs (retval.data (i));
1071  if (atmp > ainvnorm)
1072  ainvnorm = atmp;
1073  }
1074 
1075  rcond = 1. / ainvnorm / anorm;
1076  }
1077  }
1078  else
1079  (*current_liboctave_error_handler) ("incorrect matrix type");
1080  }
1081 
1082  return retval;
1083 
1084 inverse_singular:
1085  return SparseComplexMatrix ();
1086 }
1087 
1090  double& rcond, int, int calc_cond) const
1091 {
1092  int typ = mattype.type (false);
1093  SparseComplexMatrix ret;
1094 
1095  if (typ == MatrixType::Unknown)
1096  typ = mattype.type (*this);
1097 
1099  ret = dinverse (mattype, info, rcond, true, calc_cond);
1100  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1101  ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
1102  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1103  {
1104  MatrixType newtype = mattype.transpose ();
1105  ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
1106  }
1107  else
1108  {
1109  if (mattype.is_hermitian ())
1110  {
1111  MatrixType tmp_typ (MatrixType::Upper);
1112  SparseComplexCHOL fact (*this, info, false);
1113  rcond = fact.rcond ();
1114  if (info == 0)
1115  {
1116  double rcond2;
1117  SparseMatrix Q = fact.Q ();
1118  SparseComplexMatrix InvL = fact.L ().transpose ().
1119  tinverse (tmp_typ, info, rcond2,
1120  true, false);
1121  ret = Q * InvL.hermitian () * InvL * Q.transpose ();
1122  }
1123  else
1124  {
1125  // Matrix is either singular or not positive definite
1126  mattype.mark_as_unsymmetric ();
1127  }
1128  }
1129 
1130  if (!mattype.is_hermitian ())
1131  {
1132  octave_idx_type n = rows ();
1133  ColumnVector Qinit(n);
1134  for (octave_idx_type i = 0; i < n; i++)
1135  Qinit(i) = i;
1136 
1137  MatrixType tmp_typ (MatrixType::Upper);
1138  SparseComplexLU fact (*this, Qinit, Matrix (), false, false);
1139  rcond = fact.rcond ();
1140  double rcond2;
1141  SparseComplexMatrix InvL = fact.L ().transpose ().
1142  tinverse (tmp_typ, info, rcond2,
1143  true, false);
1144  SparseComplexMatrix InvU = fact.U ().
1145  tinverse (tmp_typ, info, rcond2,
1146  true, false).transpose ();
1147  ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1148  }
1149  }
1150 
1151  return ret;
1152 }
1153 
1154 ComplexDET
1156 {
1157  octave_idx_type info;
1158  double rcond;
1159  return determinant (info, rcond, 0);
1160 }
1161 
1162 ComplexDET
1164 {
1165  double rcond;
1166  return determinant (info, rcond, 0);
1167 }
1168 
1169 ComplexDET
1171  int) const
1172 {
1173  ComplexDET retval;
1174 #ifdef HAVE_UMFPACK
1175 
1176  octave_idx_type nr = rows ();
1177  octave_idx_type nc = cols ();
1178 
1179  if (nr == 0 || nc == 0 || nr != nc)
1180  {
1181  retval = ComplexDET (1.0);
1182  }
1183  else
1184  {
1185  err = 0;
1186 
1187  // Setup the control parameters
1188  Matrix Control (UMFPACK_CONTROL, 1);
1189  double *control = Control.fortran_vec ();
1190  UMFPACK_ZNAME (defaults) (control);
1191 
1192  double tmp = octave_sparse_params::get_key ("spumoni");
1193  if (!xisnan (tmp))
1194  Control (UMFPACK_PRL) = tmp;
1195 
1196  tmp = octave_sparse_params::get_key ("piv_tol");
1197  if (!xisnan (tmp))
1198  {
1199  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1200  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1201  }
1202 
1203  // Set whether we are allowed to modify Q or not
1204  tmp = octave_sparse_params::get_key ("autoamd");
1205  if (!xisnan (tmp))
1206  Control (UMFPACK_FIXQ) = tmp;
1207 
1208  // Turn-off UMFPACK scaling for LU
1209  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1210 
1211  UMFPACK_ZNAME (report_control) (control);
1212 
1213  const octave_idx_type *Ap = cidx ();
1214  const octave_idx_type *Ai = ridx ();
1215  const Complex *Ax = data ();
1216 
1217  UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
1218  reinterpret_cast<const double *> (Ax),
1219  0, 1, control);
1220 
1221  void *Symbolic;
1222  Matrix Info (1, UMFPACK_INFO);
1223  double *info = Info.fortran_vec ();
1224  int status = UMFPACK_ZNAME (qsymbolic)
1225  (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0,
1226  0, &Symbolic, control, info);
1227 
1228  if (status < 0)
1229  {
1230  (*current_liboctave_error_handler)
1231  ("SparseComplexMatrix::determinant symbolic factorization failed");
1232 
1233  UMFPACK_ZNAME (report_status) (control, status);
1234  UMFPACK_ZNAME (report_info) (control, info);
1235 
1236  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1237  }
1238  else
1239  {
1240  UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
1241 
1242  void *Numeric;
1243  status
1244  = UMFPACK_ZNAME (numeric) (Ap, Ai,
1245  reinterpret_cast<const double *> (Ax),
1246  0, Symbolic, &Numeric, control, info);
1247  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1248 
1249  rcond = Info (UMFPACK_RCOND);
1250 
1251  if (status < 0)
1252  {
1253  (*current_liboctave_error_handler)
1254  ("SparseComplexMatrix::determinant numeric factorization failed");
1255 
1256  UMFPACK_ZNAME (report_status) (control, status);
1257  UMFPACK_ZNAME (report_info) (control, info);
1258 
1259  UMFPACK_ZNAME (free_numeric) (&Numeric);
1260  }
1261  else
1262  {
1263  UMFPACK_ZNAME (report_numeric) (Numeric, control);
1264 
1265  double c10[2], e10;
1266 
1267  status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10,
1268  Numeric, info);
1269 
1270  if (status < 0)
1271  {
1272  (*current_liboctave_error_handler)
1273  ("SparseComplexMatrix::determinant error calculating determinant");
1274 
1275  UMFPACK_ZNAME (report_status) (control, status);
1276  UMFPACK_ZNAME (report_info) (control, info);
1277  }
1278  else
1279  retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10);
1280 
1281  UMFPACK_ZNAME (free_numeric) (&Numeric);
1282  }
1283  }
1284  }
1285 #else
1286  (*current_liboctave_error_handler) ("UMFPACK not installed");
1287 #endif
1288 
1289  return retval;
1290 }
1291 
1294  octave_idx_type& err, double& rcond,
1295  solve_singularity_handler, bool calc_cond) const
1296 {
1297  ComplexMatrix retval;
1298 
1299  octave_idx_type nr = rows ();
1300  octave_idx_type nc = cols ();
1301  octave_idx_type nm = (nc < nr ? nc : nr);
1302  err = 0;
1303 
1304  if (nr != b.rows ())
1306  ("matrix dimension mismatch solution of linear equations");
1307  else if (nr == 0 || nc == 0 || b.cols () == 0)
1308  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1309  else
1310  {
1311  // Print spparms("spumoni") info if requested
1312  int typ = mattype.type ();
1313  mattype.info ();
1314 
1316  {
1317  retval.resize (nc, b.cols (), Complex (0.,0.));
1318  if (typ == MatrixType::Diagonal)
1319  for (octave_idx_type j = 0; j < b.cols (); j++)
1320  for (octave_idx_type i = 0; i < nm; i++)
1321  retval(i,j) = b(i,j) / data (i);
1322  else
1323  for (octave_idx_type j = 0; j < b.cols (); j++)
1324  for (octave_idx_type k = 0; k < nc; k++)
1325  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1326  retval(k,j) = b(ridx (i),j) / data (i);
1327 
1328  if (calc_cond)
1329  {
1330  double dmax = 0.;
1331  double dmin = octave_Inf;
1332  for (octave_idx_type i = 0; i < nm; i++)
1333  {
1334  double tmp = std::abs (data (i));
1335  if (tmp > dmax)
1336  dmax = tmp;
1337  if (tmp < dmin)
1338  dmin = tmp;
1339  }
1340  rcond = dmin / dmax;
1341  }
1342  else
1343  rcond = 1.0;
1344  }
1345  else
1346  (*current_liboctave_error_handler) ("incorrect matrix type");
1347  }
1348 
1349  return retval;
1350 }
1351 
1354  octave_idx_type& err, double& rcond,
1355  solve_singularity_handler,
1356  bool calc_cond) const
1357 {
1358  SparseComplexMatrix retval;
1359 
1360  octave_idx_type nr = rows ();
1361  octave_idx_type nc = cols ();
1362  octave_idx_type nm = (nc < nr ? nc : nr);
1363  err = 0;
1364 
1365  if (nr != b.rows ())
1367  ("matrix dimension mismatch solution of linear equations");
1368  else if (nr == 0 || nc == 0 || b.cols () == 0)
1369  retval = SparseComplexMatrix (nc, b.cols ());
1370  else
1371  {
1372  // Print spparms("spumoni") info if requested
1373  int typ = mattype.type ();
1374  mattype.info ();
1375 
1377  {
1378  octave_idx_type b_nc = b.cols ();
1379  octave_idx_type b_nz = b.nnz ();
1380  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1381 
1382  retval.xcidx (0) = 0;
1383  octave_idx_type ii = 0;
1384  if (typ == MatrixType::Diagonal)
1385  for (octave_idx_type j = 0; j < b.cols (); j++)
1386  {
1387  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1388  {
1389  if (b.ridx (i) >= nm)
1390  break;
1391  retval.xridx (ii) = b.ridx (i);
1392  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1393  }
1394  retval.xcidx (j+1) = ii;
1395  }
1396  else
1397  for (octave_idx_type j = 0; j < b.cols (); j++)
1398  {
1399  for (octave_idx_type l = 0; l < nc; l++)
1400  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1401  {
1402  bool found = false;
1403  octave_idx_type k;
1404  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1405  if (ridx (i) == b.ridx (k))
1406  {
1407  found = true;
1408  break;
1409  }
1410  if (found)
1411  {
1412  retval.xridx (ii) = l;
1413  retval.xdata (ii++) = b.data (k) / data (i);
1414  }
1415  }
1416  retval.xcidx (j+1) = ii;
1417  }
1418 
1419  if (calc_cond)
1420  {
1421  double dmax = 0.;
1422  double dmin = octave_Inf;
1423  for (octave_idx_type i = 0; i < nm; i++)
1424  {
1425  double tmp = std::abs (data (i));
1426  if (tmp > dmax)
1427  dmax = tmp;
1428  if (tmp < dmin)
1429  dmin = tmp;
1430  }
1431  rcond = dmin / dmax;
1432  }
1433  else
1434  rcond = 1.0;
1435  }
1436  else
1437  (*current_liboctave_error_handler) ("incorrect matrix type");
1438  }
1439 
1440  return retval;
1441 }
1442 
1445  octave_idx_type& err, double& rcond,
1446  solve_singularity_handler,
1447  bool calc_cond) const
1448 {
1449  ComplexMatrix retval;
1450 
1451  octave_idx_type nr = rows ();
1452  octave_idx_type nc = cols ();
1453  octave_idx_type nm = (nc < nr ? nc : nr);
1454  err = 0;
1455 
1456  if (nr != b.rows ())
1458  ("matrix dimension mismatch solution of linear equations");
1459  else if (nr == 0 || nc == 0 || b.cols () == 0)
1460  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1461  else
1462  {
1463  // Print spparms("spumoni") info if requested
1464  int typ = mattype.type ();
1465  mattype.info ();
1466 
1468  {
1469  retval.resize (nc, b.cols (), Complex (0.,0.));
1470  if (typ == MatrixType::Diagonal)
1471  for (octave_idx_type j = 0; j < b.cols (); j++)
1472  for (octave_idx_type i = 0; i < nm; i++)
1473  retval(i,j) = b(i,j) / data (i);
1474  else
1475  for (octave_idx_type j = 0; j < b.cols (); j++)
1476  for (octave_idx_type k = 0; k < nc; k++)
1477  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1478  retval(k,j) = b(ridx (i),j) / data (i);
1479 
1480  if (calc_cond)
1481  {
1482  double dmax = 0.;
1483  double dmin = octave_Inf;
1484  for (octave_idx_type i = 0; i < nr; i++)
1485  {
1486  double tmp = std::abs (data (i));
1487  if (tmp > dmax)
1488  dmax = tmp;
1489  if (tmp < dmin)
1490  dmin = tmp;
1491  }
1492  rcond = dmin / dmax;
1493  }
1494  else
1495  rcond = 1.0;
1496  }
1497  else
1498  (*current_liboctave_error_handler) ("incorrect matrix type");
1499  }
1500 
1501  return retval;
1502 }
1503 
1506  octave_idx_type& err, double& rcond,
1507  solve_singularity_handler,
1508  bool calc_cond) const
1509 {
1510  SparseComplexMatrix retval;
1511 
1512  octave_idx_type nr = rows ();
1513  octave_idx_type nc = cols ();
1514  octave_idx_type nm = (nc < nr ? nc : nr);
1515  err = 0;
1516 
1517  if (nr != b.rows ())
1519  ("matrix dimension mismatch solution of linear equations");
1520  else if (nr == 0 || nc == 0 || b.cols () == 0)
1521  retval = SparseComplexMatrix (nc, b.cols ());
1522  else
1523  {
1524  // Print spparms("spumoni") info if requested
1525  int typ = mattype.type ();
1526  mattype.info ();
1527 
1529  {
1530  octave_idx_type b_nc = b.cols ();
1531  octave_idx_type b_nz = b.nnz ();
1532  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1533 
1534  retval.xcidx (0) = 0;
1535  octave_idx_type ii = 0;
1536  if (typ == MatrixType::Diagonal)
1537  for (octave_idx_type j = 0; j < b.cols (); j++)
1538  {
1539  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1540  {
1541  if (b.ridx (i) >= nm)
1542  break;
1543  retval.xridx (ii) = b.ridx (i);
1544  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1545  }
1546  retval.xcidx (j+1) = ii;
1547  }
1548  else
1549  for (octave_idx_type j = 0; j < b.cols (); j++)
1550  {
1551  for (octave_idx_type l = 0; l < nc; l++)
1552  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1553  {
1554  bool found = false;
1555  octave_idx_type k;
1556  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1557  if (ridx (i) == b.ridx (k))
1558  {
1559  found = true;
1560  break;
1561  }
1562  if (found)
1563  {
1564  retval.xridx (ii) = l;
1565  retval.xdata (ii++) = b.data (k) / data (i);
1566  }
1567  }
1568  retval.xcidx (j+1) = ii;
1569  }
1570 
1571  if (calc_cond)
1572  {
1573  double dmax = 0.;
1574  double dmin = octave_Inf;
1575  for (octave_idx_type i = 0; i < nm; i++)
1576  {
1577  double tmp = std::abs (data (i));
1578  if (tmp > dmax)
1579  dmax = tmp;
1580  if (tmp < dmin)
1581  dmin = tmp;
1582  }
1583  rcond = dmin / dmax;
1584  }
1585  else
1586  rcond = 1.0;
1587  }
1588  else
1589  (*current_liboctave_error_handler) ("incorrect matrix type");
1590  }
1591 
1592  return retval;
1593 }
1594 
1597  octave_idx_type& err, double& rcond,
1598  solve_singularity_handler sing_handler,
1599  bool calc_cond) const
1600 {
1601  ComplexMatrix retval;
1602 
1603  octave_idx_type nr = rows ();
1604  octave_idx_type nc = cols ();
1605  octave_idx_type nm = (nc > nr ? nc : nr);
1606  err = 0;
1607 
1608  if (nr != b.rows ())
1610  ("matrix dimension mismatch solution of linear equations");
1611  else if (nr == 0 || nc == 0 || b.cols () == 0)
1612  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1613  else
1614  {
1615  // Print spparms("spumoni") info if requested
1616  int typ = mattype.type ();
1617  mattype.info ();
1618 
1619  if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
1620  {
1621  double anorm = 0.;
1622  double ainvnorm = 0.;
1623  octave_idx_type b_nc = b.cols ();
1624  rcond = 1.;
1625 
1626  if (calc_cond)
1627  {
1628  // Calculate the 1-norm of matrix for rcond calculation
1629  for (octave_idx_type j = 0; j < nc; j++)
1630  {
1631  double atmp = 0.;
1632  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1633  atmp += std::abs (data (i));
1634  if (atmp > anorm)
1635  anorm = atmp;
1636  }
1637  }
1638 
1639  if (typ == MatrixType::Permuted_Upper)
1640  {
1641  retval.resize (nc, b_nc);
1642  octave_idx_type *perm = mattype.triangular_perm ();
1643  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1644 
1645  for (octave_idx_type j = 0; j < b_nc; j++)
1646  {
1647  for (octave_idx_type i = 0; i < nr; i++)
1648  work[i] = b(i,j);
1649  for (octave_idx_type i = nr; i < nc; i++)
1650  work[i] = 0.;
1651 
1652  for (octave_idx_type k = nc-1; k >= 0; k--)
1653  {
1654  octave_idx_type kidx = perm[k];
1655 
1656  if (work[k] != 0.)
1657  {
1658  if (ridx (cidx (kidx+1)-1) != k
1659  || data (cidx (kidx+1)-1) == 0.)
1660  {
1661  err = -2;
1662  goto triangular_error;
1663  }
1664 
1665  Complex tmp = work[k] / data (cidx (kidx+1)-1);
1666  work[k] = tmp;
1667  for (octave_idx_type i = cidx (kidx);
1668  i < cidx (kidx+1)-1; i++)
1669  {
1670  octave_idx_type iidx = ridx (i);
1671  work[iidx] = work[iidx] - tmp * data (i);
1672  }
1673  }
1674  }
1675 
1676  for (octave_idx_type i = 0; i < nc; i++)
1677  retval(perm[i], j) = work[i];
1678  }
1679 
1680  if (calc_cond)
1681  {
1682  // Calculation of 1-norm of inv(*this)
1683  for (octave_idx_type i = 0; i < nm; i++)
1684  work[i] = 0.;
1685 
1686  for (octave_idx_type j = 0; j < nr; j++)
1687  {
1688  work[j] = 1.;
1689 
1690  for (octave_idx_type k = j; k >= 0; k--)
1691  {
1692  octave_idx_type iidx = perm[k];
1693 
1694  if (work[k] != 0.)
1695  {
1696  Complex tmp = work[k] / data (cidx (iidx+1)-1);
1697  work[k] = tmp;
1698  for (octave_idx_type i = cidx (iidx);
1699  i < cidx (iidx+1)-1; i++)
1700  {
1701  octave_idx_type idx2 = ridx (i);
1702  work[idx2] = work[idx2] - tmp * data (i);
1703  }
1704  }
1705  }
1706  double atmp = 0;
1707  for (octave_idx_type i = 0; i < j+1; i++)
1708  {
1709  atmp += std::abs (work[i]);
1710  work[i] = 0.;
1711  }
1712  if (atmp > ainvnorm)
1713  ainvnorm = atmp;
1714  }
1715  rcond = 1. / ainvnorm / anorm;
1716  }
1717  }
1718  else
1719  {
1720  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1721  retval.resize (nc, b_nc);
1722 
1723  for (octave_idx_type j = 0; j < b_nc; j++)
1724  {
1725  for (octave_idx_type i = 0; i < nr; i++)
1726  work[i] = b(i,j);
1727  for (octave_idx_type i = nr; i < nc; i++)
1728  work[i] = 0.;
1729 
1730  for (octave_idx_type k = nc-1; k >= 0; k--)
1731  {
1732  if (work[k] != 0.)
1733  {
1734  if (ridx (cidx (k+1)-1) != k
1735  || data (cidx (k+1)-1) == 0.)
1736  {
1737  err = -2;
1738  goto triangular_error;
1739  }
1740 
1741  Complex tmp = work[k] / data (cidx (k+1)-1);
1742  work[k] = tmp;
1743  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1744  {
1745  octave_idx_type iidx = ridx (i);
1746  work[iidx] = work[iidx] - tmp * data (i);
1747  }
1748  }
1749  }
1750 
1751  for (octave_idx_type i = 0; i < nc; i++)
1752  retval.xelem (i, j) = work[i];
1753  }
1754 
1755  if (calc_cond)
1756  {
1757  // Calculation of 1-norm of inv(*this)
1758  for (octave_idx_type i = 0; i < nm; i++)
1759  work[i] = 0.;
1760 
1761  for (octave_idx_type j = 0; j < nr; j++)
1762  {
1763  work[j] = 1.;
1764 
1765  for (octave_idx_type k = j; k >= 0; k--)
1766  {
1767  if (work[k] != 0.)
1768  {
1769  Complex tmp = work[k] / data (cidx (k+1)-1);
1770  work[k] = tmp;
1771  for (octave_idx_type i = cidx (k);
1772  i < cidx (k+1)-1; i++)
1773  {
1774  octave_idx_type iidx = ridx (i);
1775  work[iidx] = work[iidx] - tmp * data (i);
1776  }
1777  }
1778  }
1779  double atmp = 0;
1780  for (octave_idx_type i = 0; i < j+1; i++)
1781  {
1782  atmp += std::abs (work[i]);
1783  work[i] = 0.;
1784  }
1785  if (atmp > ainvnorm)
1786  ainvnorm = atmp;
1787  }
1788  rcond = 1. / ainvnorm / anorm;
1789  }
1790  }
1791 
1792  triangular_error:
1793  if (err != 0)
1794  {
1795  if (sing_handler)
1796  {
1797  sing_handler (rcond);
1798  mattype.mark_as_rectangular ();
1799  }
1800  else
1801  gripe_singular_matrix (rcond);
1802  }
1803 
1804  volatile double rcond_plus_one = rcond + 1.0;
1805 
1806  if (rcond_plus_one == 1.0 || xisnan (rcond))
1807  {
1808  err = -2;
1809 
1810  if (sing_handler)
1811  {
1812  sing_handler (rcond);
1813  mattype.mark_as_rectangular ();
1814  }
1815  else
1816  gripe_singular_matrix (rcond);
1817  }
1818  }
1819  else
1820  (*current_liboctave_error_handler) ("incorrect matrix type");
1821  }
1822 
1823  return retval;
1824 }
1825 
1828  octave_idx_type& err, double& rcond,
1829  solve_singularity_handler sing_handler,
1830  bool calc_cond) const
1831 {
1832  SparseComplexMatrix retval;
1833 
1834  octave_idx_type nr = rows ();
1835  octave_idx_type nc = cols ();
1836  octave_idx_type nm = (nc > nr ? nc : nr);
1837  err = 0;
1838 
1839  if (nr != b.rows ())
1841  ("matrix dimension mismatch solution of linear equations");
1842  else if (nr == 0 || nc == 0 || b.cols () == 0)
1843  retval = SparseComplexMatrix (nc, b.cols ());
1844  else
1845  {
1846  // Print spparms("spumoni") info if requested
1847  int typ = mattype.type ();
1848  mattype.info ();
1849 
1850  if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
1851  {
1852  double anorm = 0.;
1853  double ainvnorm = 0.;
1854  rcond = 1.;
1855 
1856  if (calc_cond)
1857  {
1858  // Calculate the 1-norm of matrix for rcond calculation
1859  for (octave_idx_type j = 0; j < nc; j++)
1860  {
1861  double atmp = 0.;
1862  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1863  atmp += std::abs (data (i));
1864  if (atmp > anorm)
1865  anorm = atmp;
1866  }
1867  }
1868 
1869  octave_idx_type b_nc = b.cols ();
1870  octave_idx_type b_nz = b.nnz ();
1871  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1872  retval.xcidx (0) = 0;
1873  octave_idx_type ii = 0;
1874  octave_idx_type x_nz = b_nz;
1875 
1876  if (typ == MatrixType::Permuted_Upper)
1877  {
1878  octave_idx_type *perm = mattype.triangular_perm ();
1879  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1880 
1881  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1882  for (octave_idx_type i = 0; i < nc; i++)
1883  rperm[perm[i]] = i;
1884 
1885  for (octave_idx_type j = 0; j < b_nc; j++)
1886  {
1887  for (octave_idx_type i = 0; i < nm; i++)
1888  work[i] = 0.;
1889  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1890  work[b.ridx (i)] = b.data (i);
1891 
1892  for (octave_idx_type k = nc-1; k >= 0; k--)
1893  {
1894  octave_idx_type kidx = perm[k];
1895 
1896  if (work[k] != 0.)
1897  {
1898  if (ridx (cidx (kidx+1)-1) != k
1899  || data (cidx (kidx+1)-1) == 0.)
1900  {
1901  err = -2;
1902  goto triangular_error;
1903  }
1904 
1905  Complex tmp = work[k] / data (cidx (kidx+1)-1);
1906  work[k] = tmp;
1907  for (octave_idx_type i = cidx (kidx);
1908  i < cidx (kidx+1)-1; i++)
1909  {
1910  octave_idx_type iidx = ridx (i);
1911  work[iidx] = work[iidx] - tmp * data (i);
1912  }
1913  }
1914  }
1915 
1916  // Count nonzeros in work vector and adjust space in
1917  // retval if needed
1918  octave_idx_type new_nnz = 0;
1919  for (octave_idx_type i = 0; i < nc; i++)
1920  if (work[i] != 0.)
1921  new_nnz++;
1922 
1923  if (ii + new_nnz > x_nz)
1924  {
1925  // Resize the sparse matrix
1926  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1927  retval.change_capacity (sz);
1928  x_nz = sz;
1929  }
1930 
1931  for (octave_idx_type i = 0; i < nc; i++)
1932  if (work[rperm[i]] != 0.)
1933  {
1934  retval.xridx (ii) = i;
1935  retval.xdata (ii++) = work[rperm[i]];
1936  }
1937  retval.xcidx (j+1) = ii;
1938  }
1939 
1940  retval.maybe_compress ();
1941 
1942  if (calc_cond)
1943  {
1944  // Calculation of 1-norm of inv(*this)
1945  for (octave_idx_type i = 0; i < nm; i++)
1946  work[i] = 0.;
1947 
1948  for (octave_idx_type j = 0; j < nr; j++)
1949  {
1950  work[j] = 1.;
1951 
1952  for (octave_idx_type k = j; k >= 0; k--)
1953  {
1954  octave_idx_type iidx = perm[k];
1955 
1956  if (work[k] != 0.)
1957  {
1958  Complex tmp = work[k] / data (cidx (iidx+1)-1);
1959  work[k] = tmp;
1960  for (octave_idx_type i = cidx (iidx);
1961  i < cidx (iidx+1)-1; i++)
1962  {
1963  octave_idx_type idx2 = ridx (i);
1964  work[idx2] = work[idx2] - tmp * data (i);
1965  }
1966  }
1967  }
1968  double atmp = 0;
1969  for (octave_idx_type i = 0; i < j+1; i++)
1970  {
1971  atmp += std::abs (work[i]);
1972  work[i] = 0.;
1973  }
1974  if (atmp > ainvnorm)
1975  ainvnorm = atmp;
1976  }
1977  rcond = 1. / ainvnorm / anorm;
1978  }
1979  }
1980  else
1981  {
1982  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1983 
1984  for (octave_idx_type j = 0; j < b_nc; j++)
1985  {
1986  for (octave_idx_type i = 0; i < nm; i++)
1987  work[i] = 0.;
1988  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1989  work[b.ridx (i)] = b.data (i);
1990 
1991  for (octave_idx_type k = nc-1; k >= 0; k--)
1992  {
1993  if (work[k] != 0.)
1994  {
1995  if (ridx (cidx (k+1)-1) != k
1996  || data (cidx (k+1)-1) == 0.)
1997  {
1998  err = -2;
1999  goto triangular_error;
2000  }
2001 
2002  Complex tmp = work[k] / data (cidx (k+1)-1);
2003  work[k] = tmp;
2004  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2005  {
2006  octave_idx_type iidx = ridx (i);
2007  work[iidx] = work[iidx] - tmp * data (i);
2008  }
2009  }
2010  }
2011 
2012  // Count nonzeros in work vector and adjust space in
2013  // retval if needed
2014  octave_idx_type new_nnz = 0;
2015  for (octave_idx_type i = 0; i < nc; i++)
2016  if (work[i] != 0.)
2017  new_nnz++;
2018 
2019  if (ii + new_nnz > x_nz)
2020  {
2021  // Resize the sparse matrix
2022  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2023  retval.change_capacity (sz);
2024  x_nz = sz;
2025  }
2026 
2027  for (octave_idx_type i = 0; i < nc; i++)
2028  if (work[i] != 0.)
2029  {
2030  retval.xridx (ii) = i;
2031  retval.xdata (ii++) = work[i];
2032  }
2033  retval.xcidx (j+1) = ii;
2034  }
2035 
2036  retval.maybe_compress ();
2037 
2038  if (calc_cond)
2039  {
2040  // Calculation of 1-norm of inv(*this)
2041  for (octave_idx_type i = 0; i < nm; i++)
2042  work[i] = 0.;
2043 
2044  for (octave_idx_type j = 0; j < nr; j++)
2045  {
2046  work[j] = 1.;
2047 
2048  for (octave_idx_type k = j; k >= 0; k--)
2049  {
2050  if (work[k] != 0.)
2051  {
2052  Complex tmp = work[k] / data (cidx (k+1)-1);
2053  work[k] = tmp;
2054  for (octave_idx_type i = cidx (k);
2055  i < cidx (k+1)-1; i++)
2056  {
2057  octave_idx_type iidx = ridx (i);
2058  work[iidx] = work[iidx] - tmp * data (i);
2059  }
2060  }
2061  }
2062  double atmp = 0;
2063  for (octave_idx_type i = 0; i < j+1; i++)
2064  {
2065  atmp += std::abs (work[i]);
2066  work[i] = 0.;
2067  }
2068  if (atmp > ainvnorm)
2069  ainvnorm = atmp;
2070  }
2071  rcond = 1. / ainvnorm / anorm;
2072  }
2073  }
2074 
2075  triangular_error:
2076  if (err != 0)
2077  {
2078  if (sing_handler)
2079  {
2080  sing_handler (rcond);
2081  mattype.mark_as_rectangular ();
2082  }
2083  else
2084  gripe_singular_matrix (rcond);
2085  }
2086 
2087  volatile double rcond_plus_one = rcond + 1.0;
2088 
2089  if (rcond_plus_one == 1.0 || xisnan (rcond))
2090  {
2091  err = -2;
2092 
2093  if (sing_handler)
2094  {
2095  sing_handler (rcond);
2096  mattype.mark_as_rectangular ();
2097  }
2098  else
2099  gripe_singular_matrix (rcond);
2100  }
2101  }
2102  else
2103  (*current_liboctave_error_handler) ("incorrect matrix type");
2104  }
2105  return retval;
2106 }
2107 
2110  octave_idx_type& err, double& rcond,
2111  solve_singularity_handler sing_handler,
2112  bool calc_cond) const
2113 {
2114  ComplexMatrix retval;
2115 
2116  octave_idx_type nr = rows ();
2117  octave_idx_type nc = cols ();
2118  octave_idx_type nm = (nc > nr ? nc : nr);
2119  err = 0;
2120 
2121  if (nr != b.rows ())
2123  ("matrix dimension mismatch solution of linear equations");
2124  else if (nr == 0 || nc == 0 || b.cols () == 0)
2125  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2126  else
2127  {
2128  // Print spparms("spumoni") info if requested
2129  int typ = mattype.type ();
2130  mattype.info ();
2131 
2132  if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
2133  {
2134  double anorm = 0.;
2135  double ainvnorm = 0.;
2136  octave_idx_type b_nc = b.cols ();
2137  rcond = 1.;
2138 
2139  if (calc_cond)
2140  {
2141  // Calculate the 1-norm of matrix for rcond calculation
2142  for (octave_idx_type j = 0; j < nc; j++)
2143  {
2144  double atmp = 0.;
2145  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2146  atmp += std::abs (data (i));
2147  if (atmp > anorm)
2148  anorm = atmp;
2149  }
2150  }
2151 
2152  if (typ == MatrixType::Permuted_Upper)
2153  {
2154  retval.resize (nc, b_nc);
2155  octave_idx_type *perm = mattype.triangular_perm ();
2156  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2157 
2158  for (octave_idx_type j = 0; j < b_nc; j++)
2159  {
2160  for (octave_idx_type i = 0; i < nr; i++)
2161  work[i] = b(i,j);
2162  for (octave_idx_type i = nr; i < nc; i++)
2163  work[i] = 0.;
2164 
2165  for (octave_idx_type k = nc-1; k >= 0; k--)
2166  {
2167  octave_idx_type kidx = perm[k];
2168 
2169  if (work[k] != 0.)
2170  {
2171  if (ridx (cidx (kidx+1)-1) != k
2172  || data (cidx (kidx+1)-1) == 0.)
2173  {
2174  err = -2;
2175  goto triangular_error;
2176  }
2177 
2178  Complex tmp = work[k] / data (cidx (kidx+1)-1);
2179  work[k] = tmp;
2180  for (octave_idx_type i = cidx (kidx);
2181  i < cidx (kidx+1)-1; i++)
2182  {
2183  octave_idx_type iidx = ridx (i);
2184  work[iidx] = work[iidx] - tmp * data (i);
2185  }
2186  }
2187  }
2188 
2189  for (octave_idx_type i = 0; i < nc; i++)
2190  retval(perm[i], j) = work[i];
2191  }
2192 
2193  if (calc_cond)
2194  {
2195  // Calculation of 1-norm of inv(*this)
2196  for (octave_idx_type i = 0; i < nm; i++)
2197  work[i] = 0.;
2198 
2199  for (octave_idx_type j = 0; j < nr; j++)
2200  {
2201  work[j] = 1.;
2202 
2203  for (octave_idx_type k = j; k >= 0; k--)
2204  {
2205  octave_idx_type iidx = perm[k];
2206 
2207  if (work[k] != 0.)
2208  {
2209  Complex tmp = work[k] / data (cidx (iidx+1)-1);
2210  work[k] = tmp;
2211  for (octave_idx_type i = cidx (iidx);
2212  i < cidx (iidx+1)-1; i++)
2213  {
2214  octave_idx_type idx2 = ridx (i);
2215  work[idx2] = work[idx2] - tmp * data (i);
2216  }
2217  }
2218  }
2219  double atmp = 0;
2220  for (octave_idx_type i = 0; i < j+1; i++)
2221  {
2222  atmp += std::abs (work[i]);
2223  work[i] = 0.;
2224  }
2225  if (atmp > ainvnorm)
2226  ainvnorm = atmp;
2227  }
2228  rcond = 1. / ainvnorm / anorm;
2229  }
2230  }
2231  else
2232  {
2233  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2234  retval.resize (nc, b_nc);
2235 
2236  for (octave_idx_type j = 0; j < b_nc; j++)
2237  {
2238  for (octave_idx_type i = 0; i < nr; i++)
2239  work[i] = b(i,j);
2240  for (octave_idx_type i = nr; i < nc; i++)
2241  work[i] = 0.;
2242 
2243  for (octave_idx_type k = nc-1; k >= 0; k--)
2244  {
2245  if (work[k] != 0.)
2246  {
2247  if (ridx (cidx (k+1)-1) != k
2248  || data (cidx (k+1)-1) == 0.)
2249  {
2250  err = -2;
2251  goto triangular_error;
2252  }
2253 
2254  Complex tmp = work[k] / data (cidx (k+1)-1);
2255  work[k] = tmp;
2256  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2257  {
2258  octave_idx_type iidx = ridx (i);
2259  work[iidx] = work[iidx] - tmp * data (i);
2260  }
2261  }
2262  }
2263 
2264  for (octave_idx_type i = 0; i < nc; i++)
2265  retval.xelem (i, j) = work[i];
2266  }
2267 
2268  if (calc_cond)
2269  {
2270  // Calculation of 1-norm of inv(*this)
2271  for (octave_idx_type i = 0; i < nm; i++)
2272  work[i] = 0.;
2273 
2274  for (octave_idx_type j = 0; j < nr; j++)
2275  {
2276  work[j] = 1.;
2277 
2278  for (octave_idx_type k = j; k >= 0; k--)
2279  {
2280  if (work[k] != 0.)
2281  {
2282  Complex tmp = work[k] / data (cidx (k+1)-1);
2283  work[k] = tmp;
2284  for (octave_idx_type i = cidx (k);
2285  i < cidx (k+1)-1; i++)
2286  {
2287  octave_idx_type iidx = ridx (i);
2288  work[iidx] = work[iidx] - tmp * data (i);
2289  }
2290  }
2291  }
2292  double atmp = 0;
2293  for (octave_idx_type i = 0; i < j+1; i++)
2294  {
2295  atmp += std::abs (work[i]);
2296  work[i] = 0.;
2297  }
2298  if (atmp > ainvnorm)
2299  ainvnorm = atmp;
2300  }
2301  rcond = 1. / ainvnorm / anorm;
2302  }
2303  }
2304 
2305  triangular_error:
2306  if (err != 0)
2307  {
2308  if (sing_handler)
2309  {
2310  sing_handler (rcond);
2311  mattype.mark_as_rectangular ();
2312  }
2313  else
2314  gripe_singular_matrix (rcond);
2315  }
2316 
2317  volatile double rcond_plus_one = rcond + 1.0;
2318 
2319  if (rcond_plus_one == 1.0 || xisnan (rcond))
2320  {
2321  err = -2;
2322 
2323  if (sing_handler)
2324  {
2325  sing_handler (rcond);
2326  mattype.mark_as_rectangular ();
2327  }
2328  else
2329  gripe_singular_matrix (rcond);
2330  }
2331  }
2332  else
2333  (*current_liboctave_error_handler) ("incorrect matrix type");
2334  }
2335 
2336  return retval;
2337 }
2338 
2341  octave_idx_type& err, double& rcond,
2342  solve_singularity_handler sing_handler,
2343  bool calc_cond) const
2344 {
2345  SparseComplexMatrix retval;
2346 
2347  octave_idx_type nr = rows ();
2348  octave_idx_type nc = cols ();
2349  octave_idx_type nm = (nc > nr ? nc : nr);
2350  err = 0;
2351 
2352  if (nr != b.rows ())
2354  ("matrix dimension mismatch solution of linear equations");
2355  else if (nr == 0 || nc == 0 || b.cols () == 0)
2356  retval = SparseComplexMatrix (nc, b.cols ());
2357  else
2358  {
2359  // Print spparms("spumoni") info if requested
2360  int typ = mattype.type ();
2361  mattype.info ();
2362 
2363  if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper)
2364  {
2365  double anorm = 0.;
2366  double ainvnorm = 0.;
2367  rcond = 1.;
2368 
2369  if (calc_cond)
2370  {
2371  // Calculate the 1-norm of matrix for rcond calculation
2372  for (octave_idx_type j = 0; j < nc; j++)
2373  {
2374  double atmp = 0.;
2375  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2376  atmp += std::abs (data (i));
2377  if (atmp > anorm)
2378  anorm = atmp;
2379  }
2380  }
2381 
2382  octave_idx_type b_nc = b.cols ();
2383  octave_idx_type b_nz = b.nnz ();
2384  retval = SparseComplexMatrix (nc, b_nc, b_nz);
2385  retval.xcidx (0) = 0;
2386  octave_idx_type ii = 0;
2387  octave_idx_type x_nz = b_nz;
2388 
2389  if (typ == MatrixType::Permuted_Upper)
2390  {
2391  octave_idx_type *perm = mattype.triangular_perm ();
2392  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2393 
2394  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
2395  for (octave_idx_type i = 0; i < nc; i++)
2396  rperm[perm[i]] = i;
2397 
2398  for (octave_idx_type j = 0; j < b_nc; j++)
2399  {
2400  for (octave_idx_type i = 0; i < nm; i++)
2401  work[i] = 0.;
2402  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2403  work[b.ridx (i)] = b.data (i);
2404 
2405  for (octave_idx_type k = nc-1; k >= 0; k--)
2406  {
2407  octave_idx_type kidx = perm[k];
2408 
2409  if (work[k] != 0.)
2410  {
2411  if (ridx (cidx (kidx+1)-1) != k
2412  || data (cidx (kidx+1)-1) == 0.)
2413  {
2414  err = -2;
2415  goto triangular_error;
2416  }
2417 
2418  Complex tmp = work[k] / data (cidx (kidx+1)-1);
2419  work[k] = tmp;
2420  for (octave_idx_type i = cidx (kidx);
2421  i < cidx (kidx+1)-1; i++)
2422  {
2423  octave_idx_type iidx = ridx (i);
2424  work[iidx] = work[iidx] - tmp * data (i);
2425  }
2426  }
2427  }
2428 
2429  // Count nonzeros in work vector and adjust space in
2430  // retval if needed
2431  octave_idx_type new_nnz = 0;
2432  for (octave_idx_type i = 0; i < nc; i++)
2433  if (work[i] != 0.)
2434  new_nnz++;
2435 
2436  if (ii + new_nnz > x_nz)
2437  {
2438  // Resize the sparse matrix
2439  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2440  retval.change_capacity (sz);
2441  x_nz = sz;
2442  }
2443 
2444  for (octave_idx_type i = 0; i < nc; i++)
2445  if (work[rperm[i]] != 0.)
2446  {
2447  retval.xridx (ii) = i;
2448  retval.xdata (ii++) = work[rperm[i]];
2449  }
2450  retval.xcidx (j+1) = ii;
2451  }
2452 
2453  retval.maybe_compress ();
2454 
2455  if (calc_cond)
2456  {
2457  // Calculation of 1-norm of inv(*this)
2458  for (octave_idx_type i = 0; i < nm; i++)
2459  work[i] = 0.;
2460 
2461  for (octave_idx_type j = 0; j < nr; j++)
2462  {
2463  work[j] = 1.;
2464 
2465  for (octave_idx_type k = j; k >= 0; k--)
2466  {
2467  octave_idx_type iidx = perm[k];
2468 
2469  if (work[k] != 0.)
2470  {
2471  Complex tmp = work[k] / data (cidx (iidx+1)-1);
2472  work[k] = tmp;
2473  for (octave_idx_type i = cidx (iidx);
2474  i < cidx (iidx+1)-1; i++)
2475  {
2476  octave_idx_type idx2 = ridx (i);
2477  work[idx2] = work[idx2] - tmp * data (i);
2478  }
2479  }
2480  }
2481  double atmp = 0;
2482  for (octave_idx_type i = 0; i < j+1; i++)
2483  {
2484  atmp += std::abs (work[i]);
2485  work[i] = 0.;
2486  }
2487  if (atmp > ainvnorm)
2488  ainvnorm = atmp;
2489  }
2490  rcond = 1. / ainvnorm / anorm;
2491  }
2492  }
2493  else
2494  {
2495  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2496 
2497  for (octave_idx_type j = 0; j < b_nc; j++)
2498  {
2499  for (octave_idx_type i = 0; i < nm; i++)
2500  work[i] = 0.;
2501  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2502  work[b.ridx (i)] = b.data (i);
2503 
2504  for (octave_idx_type k = nr-1; k >= 0; k--)
2505  {
2506  if (work[k] != 0.)
2507  {
2508  if (ridx (cidx (k+1)-1) != k
2509  || data (cidx (k+1)-1) == 0.)
2510  {
2511  err = -2;
2512  goto triangular_error;
2513  }
2514 
2515  Complex tmp = work[k] / data (cidx (k+1)-1);
2516  work[k] = tmp;
2517  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2518  {
2519  octave_idx_type iidx = ridx (i);
2520  work[iidx] = work[iidx] - tmp * data (i);
2521  }
2522  }
2523  }
2524 
2525  // Count nonzeros in work vector and adjust space in
2526  // retval if needed
2527  octave_idx_type new_nnz = 0;
2528  for (octave_idx_type i = 0; i < nc; i++)
2529  if (work[i] != 0.)
2530  new_nnz++;
2531 
2532  if (ii + new_nnz > x_nz)
2533  {
2534  // Resize the sparse matrix
2535  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2536  retval.change_capacity (sz);
2537  x_nz = sz;
2538  }
2539 
2540  for (octave_idx_type i = 0; i < nc; i++)
2541  if (work[i] != 0.)
2542  {
2543  retval.xridx (ii) = i;
2544  retval.xdata (ii++) = work[i];
2545  }
2546  retval.xcidx (j+1) = ii;
2547  }
2548 
2549  retval.maybe_compress ();
2550 
2551  if (calc_cond)
2552  {
2553  // Calculation of 1-norm of inv(*this)
2554  for (octave_idx_type i = 0; i < nm; i++)
2555  work[i] = 0.;
2556 
2557  for (octave_idx_type j = 0; j < nr; j++)
2558  {
2559  work[j] = 1.;
2560 
2561  for (octave_idx_type k = j; k >= 0; k--)
2562  {
2563  if (work[k] != 0.)
2564  {
2565  Complex tmp = work[k] / data (cidx (k+1)-1);
2566  work[k] = tmp;
2567  for (octave_idx_type i = cidx (k);
2568  i < cidx (k+1)-1; i++)
2569  {
2570  octave_idx_type iidx = ridx (i);
2571  work[iidx] = work[iidx] - tmp * data (i);
2572  }
2573  }
2574  }
2575  double atmp = 0;
2576  for (octave_idx_type i = 0; i < j+1; i++)
2577  {
2578  atmp += std::abs (work[i]);
2579  work[i] = 0.;
2580  }
2581  if (atmp > ainvnorm)
2582  ainvnorm = atmp;
2583  }
2584  rcond = 1. / ainvnorm / anorm;
2585  }
2586  }
2587 
2588  triangular_error:
2589  if (err != 0)
2590  {
2591  if (sing_handler)
2592  {
2593  sing_handler (rcond);
2594  mattype.mark_as_rectangular ();
2595  }
2596  else
2597  gripe_singular_matrix (rcond);
2598  }
2599 
2600  volatile double rcond_plus_one = rcond + 1.0;
2601 
2602  if (rcond_plus_one == 1.0 || xisnan (rcond))
2603  {
2604  err = -2;
2605 
2606  if (sing_handler)
2607  {
2608  sing_handler (rcond);
2609  mattype.mark_as_rectangular ();
2610  }
2611  else
2612  gripe_singular_matrix (rcond);
2613  }
2614  }
2615  else
2616  (*current_liboctave_error_handler) ("incorrect matrix type");
2617  }
2618 
2619  return retval;
2620 }
2621 
2624  octave_idx_type& err, double& rcond,
2625  solve_singularity_handler sing_handler,
2626  bool calc_cond) const
2627 {
2628  ComplexMatrix retval;
2629 
2630  octave_idx_type nr = rows ();
2631  octave_idx_type nc = cols ();
2632  octave_idx_type nm = (nc > nr ? nc : nr);
2633  err = 0;
2634 
2635  if (nr != b.rows ())
2637  ("matrix dimension mismatch solution of linear equations");
2638  else if (nr == 0 || nc == 0 || b.cols () == 0)
2639  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2640  else
2641  {
2642  // Print spparms("spumoni") info if requested
2643  int typ = mattype.type ();
2644  mattype.info ();
2645 
2646  if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
2647  {
2648  double anorm = 0.;
2649  double ainvnorm = 0.;
2650  octave_idx_type b_nc = b.cols ();
2651  rcond = 1.;
2652 
2653  if (calc_cond)
2654  {
2655  // Calculate the 1-norm of matrix for rcond calculation
2656  for (octave_idx_type j = 0; j < nc; j++)
2657  {
2658  double atmp = 0.;
2659  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2660  atmp += std::abs (data (i));
2661  if (atmp > anorm)
2662  anorm = atmp;
2663  }
2664  }
2665 
2666  if (typ == MatrixType::Permuted_Lower)
2667  {
2668  retval.resize (nc, b_nc);
2669  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2670  octave_idx_type *perm = mattype.triangular_perm ();
2671 
2672  for (octave_idx_type j = 0; j < b_nc; j++)
2673  {
2674  for (octave_idx_type i = 0; i < nm; i++)
2675  work[i] = 0.;
2676  for (octave_idx_type i = 0; i < nr; i++)
2677  work[perm[i]] = b(i,j);
2678 
2679  for (octave_idx_type k = 0; k < nc; k++)
2680  {
2681  if (work[k] != 0.)
2682  {
2683  octave_idx_type minr = nr;
2684  octave_idx_type mini = 0;
2685 
2686  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2687  if (perm[ridx (i)] < minr)
2688  {
2689  minr = perm[ridx (i)];
2690  mini = i;
2691  }
2692 
2693  if (minr != k || data (mini) == 0.)
2694  {
2695  err = -2;
2696  goto triangular_error;
2697  }
2698 
2699  Complex tmp = work[k] / data (mini);
2700  work[k] = tmp;
2701  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2702  {
2703  if (i == mini)
2704  continue;
2705 
2706  octave_idx_type iidx = perm[ridx (i)];
2707  work[iidx] = work[iidx] - tmp * data (i);
2708  }
2709  }
2710  }
2711 
2712  for (octave_idx_type i = 0; i < nc; i++)
2713  retval(i, j) = work[i];
2714  }
2715 
2716  if (calc_cond)
2717  {
2718  // Calculation of 1-norm of inv(*this)
2719  for (octave_idx_type i = 0; i < nm; i++)
2720  work[i] = 0.;
2721 
2722  for (octave_idx_type j = 0; j < nr; j++)
2723  {
2724  work[j] = 1.;
2725 
2726  for (octave_idx_type k = 0; k < nc; k++)
2727  {
2728  if (work[k] != 0.)
2729  {
2730  octave_idx_type minr = nr;
2731  octave_idx_type mini = 0;
2732 
2733  for (octave_idx_type i = cidx (k);
2734  i < cidx (k+1); i++)
2735  if (perm[ridx (i)] < minr)
2736  {
2737  minr = perm[ridx (i)];
2738  mini = i;
2739  }
2740 
2741  Complex tmp = work[k] / data (mini);
2742  work[k] = tmp;
2743  for (octave_idx_type i = cidx (k);
2744  i < cidx (k+1); i++)
2745  {
2746  if (i == mini)
2747  continue;
2748 
2749  octave_idx_type iidx = perm[ridx (i)];
2750  work[iidx] = work[iidx] - tmp * data (i);
2751  }
2752  }
2753  }
2754 
2755  double atmp = 0;
2756  for (octave_idx_type i = j; i < nc; i++)
2757  {
2758  atmp += std::abs (work[i]);
2759  work[i] = 0.;
2760  }
2761  if (atmp > ainvnorm)
2762  ainvnorm = atmp;
2763  }
2764  rcond = 1. / ainvnorm / anorm;
2765  }
2766  }
2767  else
2768  {
2769  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2770  retval.resize (nc, b_nc, 0.);
2771 
2772  for (octave_idx_type j = 0; j < b_nc; j++)
2773  {
2774  for (octave_idx_type i = 0; i < nr; i++)
2775  work[i] = b(i,j);
2776  for (octave_idx_type i = nr; i < nc; i++)
2777  work[i] = 0.;
2778  for (octave_idx_type k = 0; k < nc; k++)
2779  {
2780  if (work[k] != 0.)
2781  {
2782  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2783  {
2784  err = -2;
2785  goto triangular_error;
2786  }
2787 
2788  Complex tmp = work[k] / data (cidx (k));
2789  work[k] = tmp;
2790  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2791  {
2792  octave_idx_type iidx = ridx (i);
2793  work[iidx] = work[iidx] - tmp * data (i);
2794  }
2795  }
2796  }
2797  for (octave_idx_type i = 0; i < nc; i++)
2798  retval.xelem (i, j) = work[i];
2799  }
2800 
2801  if (calc_cond)
2802  {
2803  // Calculation of 1-norm of inv(*this)
2804  for (octave_idx_type i = 0; i < nm; i++)
2805  work[i] = 0.;
2806 
2807  for (octave_idx_type j = 0; j < nr; j++)
2808  {
2809  work[j] = 1.;
2810 
2811  for (octave_idx_type k = j; k < nc; k++)
2812  {
2813 
2814  if (work[k] != 0.)
2815  {
2816  Complex tmp = work[k] / data (cidx (k));
2817  work[k] = tmp;
2818  for (octave_idx_type i = cidx (k)+1;
2819  i < cidx (k+1); i++)
2820  {
2821  octave_idx_type iidx = ridx (i);
2822  work[iidx] = work[iidx] - tmp * data (i);
2823  }
2824  }
2825  }
2826  double atmp = 0;
2827  for (octave_idx_type i = j; i < nc; i++)
2828  {
2829  atmp += std::abs (work[i]);
2830  work[i] = 0.;
2831  }
2832  if (atmp > ainvnorm)
2833  ainvnorm = atmp;
2834  }
2835  rcond = 1. / ainvnorm / anorm;
2836  }
2837  }
2838  triangular_error:
2839  if (err != 0)
2840  {
2841  if (sing_handler)
2842  {
2843  sing_handler (rcond);
2844  mattype.mark_as_rectangular ();
2845  }
2846  else
2847  gripe_singular_matrix (rcond);
2848  }
2849 
2850  volatile double rcond_plus_one = rcond + 1.0;
2851 
2852  if (rcond_plus_one == 1.0 || xisnan (rcond))
2853  {
2854  err = -2;
2855 
2856  if (sing_handler)
2857  {
2858  sing_handler (rcond);
2859  mattype.mark_as_rectangular ();
2860  }
2861  else
2862  gripe_singular_matrix (rcond);
2863  }
2864  }
2865  else
2866  (*current_liboctave_error_handler) ("incorrect matrix type");
2867  }
2868 
2869  return retval;
2870 }
2871 
2874  octave_idx_type& err, double& rcond,
2875  solve_singularity_handler sing_handler,
2876  bool calc_cond) const
2877 {
2878  SparseComplexMatrix retval;
2879 
2880  octave_idx_type nr = rows ();
2881  octave_idx_type nc = cols ();
2882  octave_idx_type nm = (nc > nr ? nc : nr);
2883 
2884  err = 0;
2885 
2886  if (nr != b.rows ())
2888  ("matrix dimension mismatch solution of linear equations");
2889  else if (nr == 0 || nc == 0 || b.cols () == 0)
2890  retval = SparseComplexMatrix (nc, b.cols ());
2891  else
2892  {
2893  // Print spparms("spumoni") info if requested
2894  int typ = mattype.type ();
2895  mattype.info ();
2896 
2897  if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
2898  {
2899  double anorm = 0.;
2900  double ainvnorm = 0.;
2901  rcond = 1.;
2902 
2903  if (calc_cond)
2904  {
2905  // Calculate the 1-norm of matrix for rcond calculation
2906  for (octave_idx_type j = 0; j < nc; j++)
2907  {
2908  double atmp = 0.;
2909  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2910  atmp += std::abs (data (i));
2911  if (atmp > anorm)
2912  anorm = atmp;
2913  }
2914  }
2915 
2916  octave_idx_type b_nc = b.cols ();
2917  octave_idx_type b_nz = b.nnz ();
2918  retval = SparseComplexMatrix (nc, b_nc, b_nz);
2919  retval.xcidx (0) = 0;
2920  octave_idx_type ii = 0;
2921  octave_idx_type x_nz = b_nz;
2922 
2923  if (typ == MatrixType::Permuted_Lower)
2924  {
2925  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2926  octave_idx_type *perm = mattype.triangular_perm ();
2927 
2928  for (octave_idx_type j = 0; j < b_nc; j++)
2929  {
2930  for (octave_idx_type i = 0; i < nm; i++)
2931  work[i] = 0.;
2932  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2933  work[perm[b.ridx (i)]] = b.data (i);
2934 
2935  for (octave_idx_type k = 0; k < nc; k++)
2936  {
2937  if (work[k] != 0.)
2938  {
2939  octave_idx_type minr = nr;
2940  octave_idx_type mini = 0;
2941 
2942  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2943  if (perm[ridx (i)] < minr)
2944  {
2945  minr = perm[ridx (i)];
2946  mini = i;
2947  }
2948 
2949  if (minr != k || data (mini) == 0.)
2950  {
2951  err = -2;
2952  goto triangular_error;
2953  }
2954 
2955  Complex tmp = work[k] / data (mini);
2956  work[k] = tmp;
2957  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2958  {
2959  if (i == mini)
2960  continue;
2961 
2962  octave_idx_type iidx = perm[ridx (i)];
2963  work[iidx] = work[iidx] - tmp * data (i);
2964  }
2965  }
2966  }
2967 
2968  // Count nonzeros in work vector and adjust space in
2969  // retval if needed
2970  octave_idx_type new_nnz = 0;
2971  for (octave_idx_type i = 0; i < nc; i++)
2972  if (work[i] != 0.)
2973  new_nnz++;
2974 
2975  if (ii + new_nnz > x_nz)
2976  {
2977  // Resize the sparse matrix
2978  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2979  retval.change_capacity (sz);
2980  x_nz = sz;
2981  }
2982 
2983  for (octave_idx_type i = 0; i < nc; i++)
2984  if (work[i] != 0.)
2985  {
2986  retval.xridx (ii) = i;
2987  retval.xdata (ii++) = work[i];
2988  }
2989  retval.xcidx (j+1) = ii;
2990  }
2991 
2992  retval.maybe_compress ();
2993 
2994  if (calc_cond)
2995  {
2996  // Calculation of 1-norm of inv(*this)
2997  for (octave_idx_type i = 0; i < nm; i++)
2998  work[i] = 0.;
2999 
3000  for (octave_idx_type j = 0; j < nr; j++)
3001  {
3002  work[j] = 1.;
3003 
3004  for (octave_idx_type k = 0; k < nc; k++)
3005  {
3006  if (work[k] != 0.)
3007  {
3008  octave_idx_type minr = nr;
3009  octave_idx_type mini = 0;
3010 
3011  for (octave_idx_type i = cidx (k);
3012  i < cidx (k+1); i++)
3013  if (perm[ridx (i)] < minr)
3014  {
3015  minr = perm[ridx (i)];
3016  mini = i;
3017  }
3018 
3019  Complex tmp = work[k] / data (mini);
3020  work[k] = tmp;
3021  for (octave_idx_type i = cidx (k);
3022  i < cidx (k+1); i++)
3023  {
3024  if (i == mini)
3025  continue;
3026 
3027  octave_idx_type iidx = perm[ridx (i)];
3028  work[iidx] = work[iidx] - tmp * data (i);
3029  }
3030  }
3031  }
3032 
3033  double atmp = 0;
3034  for (octave_idx_type i = j; i < nc; i++)
3035  {
3036  atmp += std::abs (work[i]);
3037  work[i] = 0.;
3038  }
3039  if (atmp > ainvnorm)
3040  ainvnorm = atmp;
3041  }
3042  rcond = 1. / ainvnorm / anorm;
3043  }
3044  }
3045  else
3046  {
3047  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3048 
3049  for (octave_idx_type j = 0; j < b_nc; j++)
3050  {
3051  for (octave_idx_type i = 0; i < nm; i++)
3052  work[i] = 0.;
3053  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3054  work[b.ridx (i)] = b.data (i);
3055 
3056  for (octave_idx_type k = 0; k < nc; k++)
3057  {
3058  if (work[k] != 0.)
3059  {
3060  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3061  {
3062  err = -2;
3063  goto triangular_error;
3064  }
3065 
3066  Complex tmp = work[k] / data (cidx (k));
3067  work[k] = tmp;
3068  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3069  {
3070  octave_idx_type iidx = ridx (i);
3071  work[iidx] = work[iidx] - tmp * data (i);
3072  }
3073  }
3074  }
3075 
3076  // Count nonzeros in work vector and adjust space in
3077  // retval if needed
3078  octave_idx_type new_nnz = 0;
3079  for (octave_idx_type i = 0; i < nc; i++)
3080  if (work[i] != 0.)
3081  new_nnz++;
3082 
3083  if (ii + new_nnz > x_nz)
3084  {
3085  // Resize the sparse matrix
3086  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3087  retval.change_capacity (sz);
3088  x_nz = sz;
3089  }
3090 
3091  for (octave_idx_type i = 0; i < nc; i++)
3092  if (work[i] != 0.)
3093  {
3094  retval.xridx (ii) = i;
3095  retval.xdata (ii++) = work[i];
3096  }
3097  retval.xcidx (j+1) = ii;
3098  }
3099 
3100  retval.maybe_compress ();
3101 
3102  if (calc_cond)
3103  {
3104  // Calculation of 1-norm of inv(*this)
3105  for (octave_idx_type i = 0; i < nm; i++)
3106  work[i] = 0.;
3107 
3108  for (octave_idx_type j = 0; j < nr; j++)
3109  {
3110  work[j] = 1.;
3111 
3112  for (octave_idx_type k = j; k < nc; k++)
3113  {
3114 
3115  if (work[k] != 0.)
3116  {
3117  Complex tmp = work[k] / data (cidx (k));
3118  work[k] = tmp;
3119  for (octave_idx_type i = cidx (k)+1;
3120  i < cidx (k+1); i++)
3121  {
3122  octave_idx_type iidx = ridx (i);
3123  work[iidx] = work[iidx] - tmp * data (i);
3124  }
3125  }
3126  }
3127  double atmp = 0;
3128  for (octave_idx_type i = j; i < nc; i++)
3129  {
3130  atmp += std::abs (work[i]);
3131  work[i] = 0.;
3132  }
3133  if (atmp > ainvnorm)
3134  ainvnorm = atmp;
3135  }
3136  rcond = 1. / ainvnorm / anorm;
3137  }
3138  }
3139 
3140  triangular_error:
3141  if (err != 0)
3142  {
3143  if (sing_handler)
3144  {
3145  sing_handler (rcond);
3146  mattype.mark_as_rectangular ();
3147  }
3148  else
3149  gripe_singular_matrix (rcond);
3150  }
3151 
3152  volatile double rcond_plus_one = rcond + 1.0;
3153 
3154  if (rcond_plus_one == 1.0 || xisnan (rcond))
3155  {
3156  err = -2;
3157 
3158  if (sing_handler)
3159  {
3160  sing_handler (rcond);
3161  mattype.mark_as_rectangular ();
3162  }
3163  else
3164  gripe_singular_matrix (rcond);
3165  }
3166  }
3167  else
3168  (*current_liboctave_error_handler) ("incorrect matrix type");
3169  }
3170 
3171  return retval;
3172 }
3173 
3176  octave_idx_type& err, double& rcond,
3177  solve_singularity_handler sing_handler,
3178  bool calc_cond) const
3179 {
3180  ComplexMatrix retval;
3181 
3182  octave_idx_type nr = rows ();
3183  octave_idx_type nc = cols ();
3184  octave_idx_type nm = (nc > nr ? nc : nr);
3185  err = 0;
3186 
3187  if (nr != b.rows ())
3189  ("matrix dimension mismatch solution of linear equations");
3190  else if (nr == 0 || nc == 0 || b.cols () == 0)
3191  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3192  else
3193  {
3194  // Print spparms("spumoni") info if requested
3195  int typ = mattype.type ();
3196  mattype.info ();
3197 
3198  if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
3199  {
3200  double anorm = 0.;
3201  double ainvnorm = 0.;
3202  octave_idx_type b_nc = b.cols ();
3203  rcond = 1.;
3204 
3205  if (calc_cond)
3206  {
3207  // Calculate the 1-norm of matrix for rcond calculation
3208  for (octave_idx_type j = 0; j < nc; j++)
3209  {
3210  double atmp = 0.;
3211  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3212  atmp += std::abs (data (i));
3213  if (atmp > anorm)
3214  anorm = atmp;
3215  }
3216  }
3217 
3218  if (typ == MatrixType::Permuted_Lower)
3219  {
3220  retval.resize (nc, b_nc);
3221  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3222  octave_idx_type *perm = mattype.triangular_perm ();
3223 
3224  for (octave_idx_type j = 0; j < b_nc; j++)
3225  {
3226  for (octave_idx_type i = 0; i < nm; i++)
3227  work[i] = 0.;
3228  for (octave_idx_type i = 0; i < nr; i++)
3229  work[perm[i]] = b(i,j);
3230 
3231  for (octave_idx_type k = 0; k < nc; k++)
3232  {
3233  if (work[k] != 0.)
3234  {
3235  octave_idx_type minr = nr;
3236  octave_idx_type mini = 0;
3237 
3238  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3239  if (perm[ridx (i)] < minr)
3240  {
3241  minr = perm[ridx (i)];
3242  mini = i;
3243  }
3244 
3245  if (minr != k || data (mini) == 0.)
3246  {
3247  err = -2;
3248  goto triangular_error;
3249  }
3250 
3251  Complex tmp = work[k] / data (mini);
3252  work[k] = tmp;
3253  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3254  {
3255  if (i == mini)
3256  continue;
3257 
3258  octave_idx_type iidx = perm[ridx (i)];
3259  work[iidx] = work[iidx] - tmp * data (i);
3260  }
3261  }
3262  }
3263 
3264  for (octave_idx_type i = 0; i < nc; i++)
3265  retval(i, j) = work[i];
3266  }
3267 
3268  if (calc_cond)
3269  {
3270  // Calculation of 1-norm of inv(*this)
3271  for (octave_idx_type i = 0; i < nm; i++)
3272  work[i] = 0.;
3273 
3274  for (octave_idx_type j = 0; j < nr; j++)
3275  {
3276  work[j] = 1.;
3277 
3278  for (octave_idx_type k = 0; k < nc; k++)
3279  {
3280  if (work[k] != 0.)
3281  {
3282  octave_idx_type minr = nr;
3283  octave_idx_type mini = 0;
3284 
3285  for (octave_idx_type i = cidx (k);
3286  i < cidx (k+1); i++)
3287  if (perm[ridx (i)] < minr)
3288  {
3289  minr = perm[ridx (i)];
3290  mini = i;
3291  }
3292 
3293  Complex tmp = work[k] / data (mini);
3294  work[k] = tmp;
3295  for (octave_idx_type i = cidx (k);
3296  i < cidx (k+1); i++)
3297  {
3298  if (i == mini)
3299  continue;
3300 
3301  octave_idx_type iidx = perm[ridx (i)];
3302  work[iidx] = work[iidx] - tmp * data (i);
3303  }
3304  }
3305  }
3306 
3307  double atmp = 0;
3308  for (octave_idx_type i = j; i < nc; i++)
3309  {
3310  atmp += std::abs (work[i]);
3311  work[i] = 0.;
3312  }
3313  if (atmp > ainvnorm)
3314  ainvnorm = atmp;
3315  }
3316  rcond = 1. / ainvnorm / anorm;
3317  }
3318  }
3319  else
3320  {
3321  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3322  retval.resize (nc, b_nc, 0.);
3323 
3324 
3325  for (octave_idx_type j = 0; j < b_nc; j++)
3326  {
3327  for (octave_idx_type i = 0; i < nr; i++)
3328  work[i] = b(i,j);
3329  for (octave_idx_type i = nr; i < nc; i++)
3330  work[i] = 0.;
3331 
3332  for (octave_idx_type k = 0; k < nc; k++)
3333  {
3334  if (work[k] != 0.)
3335  {
3336  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3337  {
3338  err = -2;
3339  goto triangular_error;
3340  }
3341 
3342  Complex tmp = work[k] / data (cidx (k));
3343  work[k] = tmp;
3344  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3345  {
3346  octave_idx_type iidx = ridx (i);
3347  work[iidx] = work[iidx] - tmp * data (i);
3348  }
3349  }
3350  }
3351 
3352  for (octave_idx_type i = 0; i < nc; i++)
3353  retval.xelem (i, j) = work[i];
3354  }
3355 
3356  if (calc_cond)
3357  {
3358  // Calculation of 1-norm of inv(*this)
3359  for (octave_idx_type i = 0; i < nm; i++)
3360  work[i] = 0.;
3361 
3362  for (octave_idx_type j = 0; j < nr; j++)
3363  {
3364  work[j] = 1.;
3365 
3366  for (octave_idx_type k = j; k < nc; k++)
3367  {
3368 
3369  if (work[k] != 0.)
3370  {
3371  Complex tmp = work[k] / data (cidx (k));
3372  work[k] = tmp;
3373  for (octave_idx_type i = cidx (k)+1;
3374  i < cidx (k+1); i++)
3375  {
3376  octave_idx_type iidx = ridx (i);
3377  work[iidx] = work[iidx] - tmp * data (i);
3378  }
3379  }
3380  }
3381  double atmp = 0;
3382  for (octave_idx_type i = j; i < nc; i++)
3383  {
3384  atmp += std::abs (work[i]);
3385  work[i] = 0.;
3386  }
3387  if (atmp > ainvnorm)
3388  ainvnorm = atmp;
3389  }
3390  rcond = 1. / ainvnorm / anorm;
3391  }
3392  }
3393 
3394  triangular_error:
3395  if (err != 0)
3396  {
3397  if (sing_handler)
3398  {
3399  sing_handler (rcond);
3400  mattype.mark_as_rectangular ();
3401  }
3402  else
3403  gripe_singular_matrix (rcond);
3404  }
3405 
3406  volatile double rcond_plus_one = rcond + 1.0;
3407 
3408  if (rcond_plus_one == 1.0 || xisnan (rcond))
3409  {
3410  err = -2;
3411 
3412  if (sing_handler)
3413  {
3414  sing_handler (rcond);
3415  mattype.mark_as_rectangular ();
3416  }
3417  else
3418  gripe_singular_matrix (rcond);
3419  }
3420  }
3421  else
3422  (*current_liboctave_error_handler) ("incorrect matrix type");
3423  }
3424 
3425  return retval;
3426 }
3427 
3430  octave_idx_type& err, double& rcond,
3431  solve_singularity_handler sing_handler,
3432  bool calc_cond) const
3433 {
3434  SparseComplexMatrix retval;
3435 
3436  octave_idx_type nr = rows ();
3437  octave_idx_type nc = cols ();
3438  octave_idx_type nm = (nc > nr ? nc : nr);
3439  err = 0;
3440 
3441  if (nr != b.rows ())
3443  ("matrix dimension mismatch solution of linear equations");
3444  else if (nr == 0 || nc == 0 || b.cols () == 0)
3445  retval = SparseComplexMatrix (nc, b.cols ());
3446  else
3447  {
3448  // Print spparms("spumoni") info if requested
3449  int typ = mattype.type ();
3450  mattype.info ();
3451 
3452  if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower)
3453  {
3454  double anorm = 0.;
3455  double ainvnorm = 0.;
3456  rcond = 1.;
3457 
3458  if (calc_cond)
3459  {
3460  // Calculate the 1-norm of matrix for rcond calculation
3461  for (octave_idx_type j = 0; j < nc; j++)
3462  {
3463  double atmp = 0.;
3464  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3465  atmp += std::abs (data (i));
3466  if (atmp > anorm)
3467  anorm = atmp;
3468  }
3469  }
3470 
3471  octave_idx_type b_nc = b.cols ();
3472  octave_idx_type b_nz = b.nnz ();
3473  retval = SparseComplexMatrix (nc, b_nc, b_nz);
3474  retval.xcidx (0) = 0;
3475  octave_idx_type ii = 0;
3476  octave_idx_type x_nz = b_nz;
3477 
3478  if (typ == MatrixType::Permuted_Lower)
3479  {
3480  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3481  octave_idx_type *perm = mattype.triangular_perm ();
3482 
3483  for (octave_idx_type j = 0; j < b_nc; j++)
3484  {
3485  for (octave_idx_type i = 0; i < nm; i++)
3486  work[i] = 0.;
3487  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3488  work[perm[b.ridx (i)]] = b.data (i);
3489 
3490  for (octave_idx_type k = 0; k < nc; k++)
3491  {
3492  if (work[k] != 0.)
3493  {
3494  octave_idx_type minr = nr;
3495  octave_idx_type mini = 0;
3496 
3497  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3498  if (perm[ridx (i)] < minr)
3499  {
3500  minr = perm[ridx (i)];
3501  mini = i;
3502  }
3503 
3504  if (minr != k || data (mini) == 0.)
3505  {
3506  err = -2;
3507  goto triangular_error;
3508  }
3509 
3510  Complex tmp = work[k] / data (mini);
3511  work[k] = tmp;
3512  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3513  {
3514  if (i == mini)
3515  continue;
3516 
3517  octave_idx_type iidx = perm[ridx (i)];
3518  work[iidx] = work[iidx] - tmp * data (i);
3519  }
3520  }
3521  }
3522 
3523  // Count nonzeros in work vector and adjust space in
3524  // retval if needed
3525  octave_idx_type new_nnz = 0;
3526  for (octave_idx_type i = 0; i < nc; i++)
3527  if (work[i] != 0.)
3528  new_nnz++;
3529 
3530  if (ii + new_nnz > x_nz)
3531  {
3532  // Resize the sparse matrix
3533  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3534  retval.change_capacity (sz);
3535  x_nz = sz;
3536  }
3537 
3538  for (octave_idx_type i = 0; i < nc; i++)
3539  if (work[i] != 0.)
3540  {
3541  retval.xridx (ii) = i;
3542  retval.xdata (ii++) = work[i];
3543  }
3544  retval.xcidx (j+1) = ii;
3545  }
3546 
3547  retval.maybe_compress ();
3548 
3549  if (calc_cond)
3550  {
3551  // Calculation of 1-norm of inv(*this)
3552  for (octave_idx_type i = 0; i < nm; i++)
3553  work[i] = 0.;
3554 
3555  for (octave_idx_type j = 0; j < nr; j++)
3556  {
3557  work[j] = 1.;
3558 
3559  for (octave_idx_type k = 0; k < nc; k++)
3560  {
3561  if (work[k] != 0.)
3562  {
3563  octave_idx_type minr = nr;
3564  octave_idx_type mini = 0;
3565 
3566  for (octave_idx_type i = cidx (k);
3567  i < cidx (k+1); i++)
3568  if (perm[ridx (i)] < minr)
3569  {
3570  minr = perm[ridx (i)];
3571  mini = i;
3572  }
3573 
3574  Complex tmp = work[k] / data (mini);
3575  work[k] = tmp;
3576  for (octave_idx_type i = cidx (k);
3577  i < cidx (k+1); i++)
3578  {
3579  if (i == mini)
3580  continue;
3581 
3582  octave_idx_type iidx = perm[ridx (i)];
3583  work[iidx] = work[iidx] - tmp * data (i);
3584  }
3585  }
3586  }
3587 
3588  double atmp = 0;
3589  for (octave_idx_type i = j; i < nc; i++)
3590  {
3591  atmp += std::abs (work[i]);
3592  work[i] = 0.;
3593  }
3594  if (atmp > ainvnorm)
3595  ainvnorm = atmp;
3596  }
3597  rcond = 1. / ainvnorm / anorm;
3598  }
3599  }
3600  else
3601  {
3602  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3603 
3604  for (octave_idx_type j = 0; j < b_nc; j++)
3605  {
3606  for (octave_idx_type i = 0; i < nm; i++)
3607  work[i] = 0.;
3608  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3609  work[b.ridx (i)] = b.data (i);
3610 
3611  for (octave_idx_type k = 0; k < nc; k++)
3612  {
3613  if (work[k] != 0.)
3614  {
3615  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3616  {
3617  err = -2;
3618  goto triangular_error;
3619  }
3620 
3621  Complex tmp = work[k] / data (cidx (k));
3622  work[k] = tmp;
3623  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3624  {
3625  octave_idx_type iidx = ridx (i);
3626  work[iidx] = work[iidx] - tmp * data (i);
3627  }
3628  }
3629  }
3630 
3631  // Count nonzeros in work vector and adjust space in
3632  // retval if needed
3633  octave_idx_type new_nnz = 0;
3634  for (octave_idx_type i = 0; i < nc; i++)
3635  if (work[i] != 0.)
3636  new_nnz++;
3637 
3638  if (ii + new_nnz > x_nz)
3639  {
3640  // Resize the sparse matrix
3641  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3642  retval.change_capacity (sz);
3643  x_nz = sz;
3644  }
3645 
3646  for (octave_idx_type i = 0; i < nc; i++)
3647  if (work[i] != 0.)
3648  {
3649  retval.xridx (ii) = i;
3650  retval.xdata (ii++) = work[i];
3651  }
3652  retval.xcidx (j+1) = ii;
3653  }
3654 
3655  retval.maybe_compress ();
3656 
3657  if (calc_cond)
3658  {
3659  // Calculation of 1-norm of inv(*this)
3660  for (octave_idx_type i = 0; i < nm; i++)
3661  work[i] = 0.;
3662 
3663  for (octave_idx_type j = 0; j < nr; j++)
3664  {
3665  work[j] = 1.;
3666 
3667  for (octave_idx_type k = j; k < nc; k++)
3668  {
3669 
3670  if (work[k] != 0.)
3671  {
3672  Complex tmp = work[k] / data (cidx (k));
3673  work[k] = tmp;
3674  for (octave_idx_type i = cidx (k)+1;
3675  i < cidx (k+1); i++)
3676  {
3677  octave_idx_type iidx = ridx (i);
3678  work[iidx] = work[iidx] - tmp * data (i);
3679  }
3680  }
3681  }
3682  double atmp = 0;
3683  for (octave_idx_type i = j; i < nc; i++)
3684  {
3685  atmp += std::abs (work[i]);
3686  work[i] = 0.;
3687  }
3688  if (atmp > ainvnorm)
3689  ainvnorm = atmp;
3690  }
3691  rcond = 1. / ainvnorm / anorm;
3692  }
3693  }
3694 
3695  triangular_error:
3696  if (err != 0)
3697  {
3698  if (sing_handler)
3699  {
3700  sing_handler (rcond);
3701  mattype.mark_as_rectangular ();
3702  }
3703  else
3704  gripe_singular_matrix (rcond);
3705  }
3706 
3707  volatile double rcond_plus_one = rcond + 1.0;
3708 
3709  if (rcond_plus_one == 1.0 || xisnan (rcond))
3710  {
3711  err = -2;
3712 
3713  if (sing_handler)
3714  {
3715  sing_handler (rcond);
3716  mattype.mark_as_rectangular ();
3717  }
3718  else
3719  gripe_singular_matrix (rcond);
3720  }
3721  }
3722  else
3723  (*current_liboctave_error_handler) ("incorrect matrix type");
3724  }
3725 
3726  return retval;
3727 }
3728 
3731  octave_idx_type& err, double& rcond,
3732  solve_singularity_handler sing_handler,
3733  bool calc_cond) const
3734 {
3735  ComplexMatrix retval;
3736 
3737  octave_idx_type nr = rows ();
3738  octave_idx_type nc = cols ();
3739  err = 0;
3740 
3741  if (nr != nc || nr != b.rows ())
3743  ("matrix dimension mismatch solution of linear equations");
3744  else if (nr == 0 || b.cols () == 0)
3745  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3746  else if (calc_cond)
3747  (*current_liboctave_error_handler)
3748  ("calculation of condition number not implemented");
3749  else
3750  {
3751  // Print spparms("spumoni") info if requested
3752  volatile int typ = mattype.type ();
3753  mattype.info ();
3754 
3756  {
3757  OCTAVE_LOCAL_BUFFER (double, D, nr);
3758  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3759 
3760  if (mattype.is_dense ())
3761  {
3762  octave_idx_type ii = 0;
3763 
3764  for (octave_idx_type j = 0; j < nc-1; j++)
3765  {
3766  D[j] = std::real (data (ii++));
3767  DL[j] = data (ii);
3768  ii += 2;
3769  }
3770  D[nc-1] = std::real (data (ii));
3771  }
3772  else
3773  {
3774  D[0] = 0.;
3775  for (octave_idx_type i = 0; i < nr - 1; i++)
3776  {
3777  D[i+1] = 0.;
3778  DL[i] = 0.;
3779  }
3780 
3781  for (octave_idx_type j = 0; j < nc; j++)
3782  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3783  {
3784  if (ridx (i) == j)
3785  D[j] = std::real (data (i));
3786  else if (ridx (i) == j + 1)
3787  DL[j] = data (i);
3788  }
3789  }
3790 
3791  octave_idx_type b_nc = b.cols ();
3792  retval = ComplexMatrix (b);
3793  Complex *result = retval.fortran_vec ();
3794 
3795  F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
3796  b.rows (), err));
3797 
3798  if (err != 0)
3799  {
3800  err = 0;
3801  mattype.mark_as_unsymmetric ();
3803  }
3804  else
3805  rcond = 1.;
3806  }
3807 
3808  if (typ == MatrixType::Tridiagonal)
3809  {
3810  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3811  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3812  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3813 
3814  if (mattype.is_dense ())
3815  {
3816  octave_idx_type ii = 0;
3817 
3818  for (octave_idx_type j = 0; j < nc-1; j++)
3819  {
3820  D[j] = data (ii++);
3821  DL[j] = data (ii++);
3822  DU[j] = data (ii++);
3823  }
3824  D[nc-1] = data (ii);
3825  }
3826  else
3827  {
3828  D[0] = 0.;
3829  for (octave_idx_type i = 0; i < nr - 1; i++)
3830  {
3831  D[i+1] = 0.;
3832  DL[i] = 0.;
3833  DU[i] = 0.;
3834  }
3835 
3836  for (octave_idx_type j = 0; j < nc; j++)
3837  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3838  {
3839  if (ridx (i) == j)
3840  D[j] = data (i);
3841  else if (ridx (i) == j + 1)
3842  DL[j] = data (i);
3843  else if (ridx (i) == j - 1)
3844  DU[j-1] = data (i);
3845  }
3846  }
3847 
3848  octave_idx_type b_nc = b.cols ();
3849  retval = ComplexMatrix (b);
3850  Complex *result = retval.fortran_vec ();
3851 
3852  F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
3853  b.rows (), err));
3854 
3855  if (err != 0)
3856  {
3857  rcond = 0.;
3858  err = -2;
3859 
3860  if (sing_handler)
3861  {
3862  sing_handler (rcond);
3863  mattype.mark_as_rectangular ();
3864  }
3865  else
3867 
3868  }
3869  else
3870  rcond = 1.;
3871  }
3872  else if (typ != MatrixType::Tridiagonal_Hermitian)
3873  (*current_liboctave_error_handler) ("incorrect matrix type");
3874  }
3875 
3876  return retval;
3877 }
3878 
3881  octave_idx_type& err, double& rcond,
3882  solve_singularity_handler sing_handler,
3883  bool calc_cond) const
3884 {
3885  SparseComplexMatrix retval;
3886 
3887  octave_idx_type nr = rows ();
3888  octave_idx_type nc = cols ();
3889  err = 0;
3890 
3891  if (nr != nc || nr != b.rows ())
3893  ("matrix dimension mismatch solution of linear equations");
3894  else if (nr == 0 || b.cols () == 0)
3895  retval = SparseComplexMatrix (nc, b.cols ());
3896  else if (calc_cond)
3897  (*current_liboctave_error_handler)
3898  ("calculation of condition number not implemented");
3899  else
3900  {
3901  // Print spparms("spumoni") info if requested
3902  int typ = mattype.type ();
3903  mattype.info ();
3904 
3905  // Note can't treat symmetric case as there is no dpttrf function
3906  if (typ == MatrixType::Tridiagonal
3908  {
3909  OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
3910  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3911  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3912  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3913  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
3914  octave_idx_type *pipvt = ipvt.fortran_vec ();
3915 
3916  if (mattype.is_dense ())
3917  {
3918  octave_idx_type ii = 0;
3919 
3920  for (octave_idx_type j = 0; j < nc-1; j++)
3921  {
3922  D[j] = data (ii++);
3923  DL[j] = data (ii++);
3924  DU[j] = data (ii++);
3925  }
3926  D[nc-1] = data (ii);
3927  }
3928  else
3929  {
3930  D[0] = 0.;
3931  for (octave_idx_type i = 0; i < nr - 1; i++)
3932  {
3933  D[i+1] = 0.;
3934  DL[i] = 0.;
3935  DU[i] = 0.;
3936  }
3937 
3938  for (octave_idx_type j = 0; j < nc; j++)
3939  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3940  {
3941  if (ridx (i) == j)
3942  D[j] = data (i);
3943  else if (ridx (i) == j + 1)
3944  DL[j] = data (i);
3945  else if (ridx (i) == j - 1)
3946  DU[j-1] = data (i);
3947  }
3948  }
3949 
3950  F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
3951 
3952  if (err != 0)
3953  {
3954  err = -2;
3955  rcond = 0.0;
3956 
3957  if (sing_handler)
3958  {
3959  sing_handler (rcond);
3960  mattype.mark_as_rectangular ();
3961  }
3962  else
3964  }
3965  else
3966  {
3967  char job = 'N';
3968  volatile octave_idx_type x_nz = b.nnz ();
3969  octave_idx_type b_nc = b.cols ();
3970  retval = SparseComplexMatrix (nr, b_nc, x_nz);
3971  retval.xcidx (0) = 0;
3972  volatile octave_idx_type ii = 0;
3973  rcond = 1.0;
3974 
3975  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
3976 
3977  for (volatile octave_idx_type j = 0; j < b_nc; j++)
3978  {
3979  for (octave_idx_type i = 0; i < nr; i++)
3980  work[i] = 0.;
3981  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3982  work[b.ridx (i)] = b.data (i);
3983 
3984  F77_XFCN (zgttrs, ZGTTRS,
3985  (F77_CONST_CHAR_ARG2 (&job, 1),
3986  nr, 1, DL, D, DU, DU2, pipvt,
3987  work, b.rows (), err
3988  F77_CHAR_ARG_LEN (1)));
3989 
3990  // Count nonzeros in work vector and adjust
3991  // space in retval if needed
3992  octave_idx_type new_nnz = 0;
3993  for (octave_idx_type i = 0; i < nr; i++)
3994  if (work[i] != 0.)
3995  new_nnz++;
3996 
3997  if (ii + new_nnz > x_nz)
3998  {
3999  // Resize the sparse matrix
4000  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4001  retval.change_capacity (sz);
4002  x_nz = sz;
4003  }
4004 
4005  for (octave_idx_type i = 0; i < nr; i++)
4006  if (work[i] != 0.)
4007  {
4008  retval.xridx (ii) = i;
4009  retval.xdata (ii++) = work[i];
4010  }
4011  retval.xcidx (j+1) = ii;
4012  }
4013 
4014  retval.maybe_compress ();
4015  }
4016  }
4017  else if (typ != MatrixType::Tridiagonal_Hermitian)
4018  (*current_liboctave_error_handler) ("incorrect matrix type");
4019  }
4020 
4021  return retval;
4022 }
4023 
4026  octave_idx_type& err, double& rcond,
4027  solve_singularity_handler sing_handler,
4028  bool calc_cond) const
4029 {
4030  ComplexMatrix retval;
4031 
4032  octave_idx_type nr = rows ();
4033  octave_idx_type nc = cols ();
4034  err = 0;
4035 
4036  if (nr != nc || nr != b.rows ())
4038  ("matrix dimension mismatch solution of linear equations");
4039  else if (nr == 0 || b.cols () == 0)
4040  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4041  else if (calc_cond)
4042  (*current_liboctave_error_handler)
4043  ("calculation of condition number not implemented");
4044  else
4045  {
4046  // Print spparms("spumoni") info if requested
4047  volatile int typ = mattype.type ();
4048  mattype.info ();
4049 
4051  {
4052  OCTAVE_LOCAL_BUFFER (double, D, nr);
4053  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4054 
4055  if (mattype.is_dense ())
4056  {
4057  octave_idx_type ii = 0;
4058 
4059  for (octave_idx_type j = 0; j < nc-1; j++)
4060  {
4061  D[j] = std::real (data (ii++));
4062  DL[j] = data (ii);
4063  ii += 2;
4064  }
4065  D[nc-1] = std::real (data (ii));
4066  }
4067  else
4068  {
4069  D[0] = 0.;
4070  for (octave_idx_type i = 0; i < nr - 1; i++)
4071  {
4072  D[i+1] = 0.;
4073  DL[i] = 0.;
4074  }
4075 
4076  for (octave_idx_type j = 0; j < nc; j++)
4077  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4078  {
4079  if (ridx (i) == j)
4080  D[j] = std::real (data (i));
4081  else if (ridx (i) == j + 1)
4082  DL[j] = data (i);
4083  }
4084  }
4085 
4086  octave_idx_type b_nr = b.rows ();
4087  octave_idx_type b_nc = b.cols ();
4088  rcond = 1.;
4089 
4090  retval = ComplexMatrix (b);
4091  Complex *result = retval.fortran_vec ();
4092 
4093  F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
4094  b_nr, err));
4095 
4096  if (err != 0)
4097  {
4098  err = 0;
4099  mattype.mark_as_unsymmetric ();
4101  }
4102  }
4103 
4104  if (typ == MatrixType::Tridiagonal)
4105  {
4106  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4107  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4108  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4109 
4110  if (mattype.is_dense ())
4111  {
4112  octave_idx_type ii = 0;
4113 
4114  for (octave_idx_type j = 0; j < nc-1; j++)
4115  {
4116  D[j] = data (ii++);
4117  DL[j] = data (ii++);
4118  DU[j] = data (ii++);
4119  }
4120  D[nc-1] = data (ii);
4121  }
4122  else
4123  {
4124  D[0] = 0.;
4125  for (octave_idx_type i = 0; i < nr - 1; i++)
4126  {
4127  D[i+1] = 0.;
4128  DL[i] = 0.;
4129  DU[i] = 0.;
4130  }
4131 
4132  for (octave_idx_type j = 0; j < nc; j++)
4133  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4134  {
4135  if (ridx (i) == j)
4136  D[j] = data (i);
4137  else if (ridx (i) == j + 1)
4138  DL[j] = data (i);
4139  else if (ridx (i) == j - 1)
4140  DU[j-1] = data (i);
4141  }
4142  }
4143 
4144  octave_idx_type b_nr = b.rows ();
4145  octave_idx_type b_nc = b.cols ();
4146  rcond = 1.;
4147 
4148  retval = ComplexMatrix (b);
4149  Complex *result = retval.fortran_vec ();
4150 
4151  F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
4152  b_nr, err));
4153 
4154  if (err != 0)
4155  {
4156  rcond = 0.;
4157  err = -2;
4158 
4159  if (sing_handler)
4160  {
4161  sing_handler (rcond);
4162  mattype.mark_as_rectangular ();
4163  }
4164  else
4166  }
4167  }
4168  else if (typ != MatrixType::Tridiagonal_Hermitian)
4169  (*current_liboctave_error_handler) ("incorrect matrix type");
4170  }
4171 
4172  return retval;
4173 }
4174 
4177  const SparseComplexMatrix& b,
4178  octave_idx_type& err, double& rcond,
4179  solve_singularity_handler sing_handler,
4180  bool calc_cond) const
4181 {
4182  SparseComplexMatrix retval;
4183 
4184  octave_idx_type nr = rows ();
4185  octave_idx_type nc = cols ();
4186  err = 0;
4187 
4188  if (nr != nc || nr != b.rows ())
4190  ("matrix dimension mismatch solution of linear equations");
4191  else if (nr == 0 || b.cols () == 0)
4192  retval = SparseComplexMatrix (nc, b.cols ());
4193  else if (calc_cond)
4194  (*current_liboctave_error_handler)
4195  ("calculation of condition number not implemented");
4196  else
4197  {
4198  // Print spparms("spumoni") info if requested
4199  int typ = mattype.type ();
4200  mattype.info ();
4201 
4202  // Note can't treat symmetric case as there is no dpttrf function
4203  if (typ == MatrixType::Tridiagonal
4205  {
4206  OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
4207  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4208  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4209  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4210  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4211  octave_idx_type *pipvt = ipvt.fortran_vec ();
4212 
4213  if (mattype.is_dense ())
4214  {
4215  octave_idx_type ii = 0;
4216 
4217  for (octave_idx_type j = 0; j < nc-1; j++)
4218  {
4219  D[j] = data (ii++);
4220  DL[j] = data (ii++);
4221  DU[j] = data (ii++);
4222  }
4223  D[nc-1] = data (ii);
4224  }
4225  else
4226  {
4227  D[0] = 0.;
4228  for (octave_idx_type i = 0; i < nr - 1; i++)
4229  {
4230  D[i+1] = 0.;
4231  DL[i] = 0.;
4232  DU[i] = 0.;
4233  }
4234 
4235  for (octave_idx_type j = 0; j < nc; j++)
4236  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4237  {
4238  if (ridx (i) == j)
4239  D[j] = data (i);
4240  else if (ridx (i) == j + 1)
4241  DL[j] = data (i);
4242  else if (ridx (i) == j - 1)
4243  DU[j-1] = data (i);
4244  }
4245  }
4246 
4247  F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4248 
4249  if (err != 0)
4250  {
4251  rcond = 0.0;
4252  err = -2;
4253 
4254  if (sing_handler)
4255  {
4256  sing_handler (rcond);
4257  mattype.mark_as_rectangular ();
4258  }
4259  else
4261  }
4262  else
4263  {
4264  rcond = 1.;
4265  char job = 'N';
4266  octave_idx_type b_nr = b.rows ();
4267  octave_idx_type b_nc = b.cols ();
4268  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4269 
4270  // Take a first guess that the number of nonzero terms
4271  // will be as many as in b
4272  volatile octave_idx_type x_nz = b.nnz ();
4273  volatile octave_idx_type ii = 0;
4274  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4275 
4276  retval.xcidx (0) = 0;
4277  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4278  {
4279 
4280  for (octave_idx_type i = 0; i < b_nr; i++)
4281  Bx[i] = b(i,j);
4282 
4283  F77_XFCN (zgttrs, ZGTTRS,
4284  (F77_CONST_CHAR_ARG2 (&job, 1),
4285  nr, 1, DL, D, DU, DU2, pipvt,
4286  Bx, b_nr, err
4287  F77_CHAR_ARG_LEN (1)));
4288 
4289  if (err != 0)
4290  {
4291  (*current_liboctave_error_handler)
4292  ("SparseComplexMatrix::solve solve failed");
4293 
4294  err = -1;
4295  break;
4296  }
4297 
4298  // Count nonzeros in work vector and adjust
4299  // space in retval if needed
4300  octave_idx_type new_nnz = 0;
4301  for (octave_idx_type i = 0; i < nr; i++)
4302  if (Bx[i] != 0.)
4303  new_nnz++;
4304 
4305  if (ii + new_nnz > x_nz)
4306  {
4307  // Resize the sparse matrix
4308  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4309  retval.change_capacity (sz);
4310  x_nz = sz;
4311  }
4312 
4313  for (octave_idx_type i = 0; i < nr; i++)
4314  if (Bx[i] != 0.)
4315  {
4316  retval.xridx (ii) = i;
4317  retval.xdata (ii++) = Bx[i];
4318  }
4319 
4320  retval.xcidx (j+1) = ii;
4321  }
4322 
4323  retval.maybe_compress ();
4324  }
4325  }
4326  else if (typ != MatrixType::Tridiagonal_Hermitian)
4327  (*current_liboctave_error_handler) ("incorrect matrix type");
4328  }
4329 
4330  return retval;
4331 }
4332 
4335  octave_idx_type& err, double& rcond,
4336  solve_singularity_handler sing_handler,
4337  bool calc_cond) const
4338 {
4339  ComplexMatrix retval;
4340 
4341  octave_idx_type nr = rows ();
4342  octave_idx_type nc = cols ();
4343  err = 0;
4344 
4345  if (nr != nc || nr != b.rows ())
4347  ("matrix dimension mismatch solution of linear equations");
4348  else if (nr == 0 || b.cols () == 0)
4349  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4350  else
4351  {
4352  // Print spparms("spumoni") info if requested
4353  volatile int typ = mattype.type ();
4354  mattype.info ();
4355 
4356  if (typ == MatrixType::Banded_Hermitian)
4357  {
4358  octave_idx_type n_lower = mattype.nlower ();
4359  octave_idx_type ldm = n_lower + 1;
4360  ComplexMatrix m_band (ldm, nc);
4361  Complex *tmp_data = m_band.fortran_vec ();
4362 
4363  if (! mattype.is_dense ())
4364  {
4365  octave_idx_type ii = 0;
4366 
4367  for (octave_idx_type j = 0; j < ldm; j++)
4368  for (octave_idx_type i = 0; i < nc; i++)
4369  tmp_data[ii++] = 0.;
4370  }
4371 
4372  for (octave_idx_type j = 0; j < nc; j++)
4373  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4374  {
4375  octave_idx_type ri = ridx (i);
4376  if (ri >= j)
4377  m_band(ri - j, j) = data (i);
4378  }
4379 
4380  // Calculate the norm of the matrix, for later use.
4381  double anorm;
4382  if (calc_cond)
4383  anorm = m_band.abs ().sum ().row (0).max ();
4384 
4385  char job = 'L';
4386  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4387  nr, n_lower, tmp_data, ldm, err
4388  F77_CHAR_ARG_LEN (1)));
4389 
4390  if (err != 0)
4391  {
4392  rcond = 0.0;
4393  // Matrix is not positive definite!! Fall through to
4394  // unsymmetric banded solver.
4395  mattype.mark_as_unsymmetric ();
4396  typ = MatrixType::Banded;
4397  err = 0;
4398  }
4399  else
4400  {
4401  if (calc_cond)
4402  {
4403  Array<Complex> z (dim_vector (2 * nr, 1));
4404  Complex *pz = z.fortran_vec ();
4405  Array<double> iz (dim_vector (nr, 1));
4406  double *piz = iz.fortran_vec ();
4407 
4408  F77_XFCN (zpbcon, ZPBCON,
4409  (F77_CONST_CHAR_ARG2 (&job, 1),
4410  nr, n_lower, tmp_data, ldm,
4411  anorm, rcond, pz, piz, err
4412  F77_CHAR_ARG_LEN (1)));
4413 
4414  if (err != 0)
4415  err = -2;
4416 
4417  volatile double rcond_plus_one = rcond + 1.0;
4418 
4419  if (rcond_plus_one == 1.0 || xisnan (rcond))
4420  {
4421  err = -2;
4422 
4423  if (sing_handler)
4424  {
4425  sing_handler (rcond);
4426  mattype.mark_as_rectangular ();
4427  }
4428  else
4429  gripe_singular_matrix (rcond);
4430  }
4431  }
4432  else
4433  rcond = 1.0;
4434 
4435  if (err == 0)
4436  {
4437  retval = ComplexMatrix (b);
4438  Complex *result = retval.fortran_vec ();
4439 
4440  octave_idx_type b_nc = b.cols ();
4441 
4442  F77_XFCN (zpbtrs, ZPBTRS,
4443  (F77_CONST_CHAR_ARG2 (&job, 1),
4444  nr, n_lower, b_nc, tmp_data,
4445  ldm, result, b.rows (), err
4446  F77_CHAR_ARG_LEN (1)));
4447 
4448  if (err != 0)
4449  {
4450  (*current_liboctave_error_handler)
4451  ("SparseMatrix::solve solve failed");
4452  err = -1;
4453  }
4454  }
4455  }
4456  }
4457 
4458  if (typ == MatrixType::Banded)
4459  {
4460  // Create the storage for the banded form of the sparse matrix
4461  octave_idx_type n_upper = mattype.nupper ();
4462  octave_idx_type n_lower = mattype.nlower ();
4463  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4464 
4465  ComplexMatrix m_band (ldm, nc);
4466  Complex *tmp_data = m_band.fortran_vec ();
4467 
4468  if (! mattype.is_dense ())
4469  {
4470  octave_idx_type ii = 0;
4471 
4472  for (octave_idx_type j = 0; j < ldm; j++)
4473  for (octave_idx_type i = 0; i < nc; i++)
4474  tmp_data[ii++] = 0.;
4475  }
4476 
4477  for (octave_idx_type j = 0; j < nc; j++)
4478  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4479  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4480 
4481  // Calculate the norm of the matrix, for later use.
4482  double anorm = 0.0;
4483  if (calc_cond)
4484  {
4485  for (octave_idx_type j = 0; j < nr; j++)
4486  {
4487  double atmp = 0.;
4488  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4489  atmp += std::abs (data (i));
4490  if (atmp > anorm)
4491  anorm = atmp;
4492  }
4493  }
4494 
4495  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4496  octave_idx_type *pipvt = ipvt.fortran_vec ();
4497 
4498  F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data,
4499  ldm, pipvt, err));
4500 
4501  // Throw-away extra info LAPACK gives so as to not
4502  // change output.
4503  if (err != 0)
4504  {
4505  rcond = 0.0;
4506  err = -2;
4507 
4508  if (sing_handler)
4509  {
4510  sing_handler (rcond);
4511  mattype.mark_as_rectangular ();
4512  }
4513  else
4515  }
4516  else
4517  {
4518  if (calc_cond)
4519  {
4520  char job = '1';
4521  Array<Complex> z (dim_vector (2 * nr, 1));
4522  Complex *pz = z.fortran_vec ();
4523  Array<double> iz (dim_vector (nr, 1));
4524  double *piz = iz.fortran_vec ();
4525 
4526  F77_XFCN (zgbcon, ZGBCON,
4527  (F77_CONST_CHAR_ARG2 (&job, 1),
4528  nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4529  anorm, rcond, pz, piz, err
4530  F77_CHAR_ARG_LEN (1)));
4531 
4532  if (err != 0)
4533  err = -2;
4534 
4535  volatile double rcond_plus_one = rcond + 1.0;
4536 
4537  if (rcond_plus_one == 1.0 || xisnan (rcond))
4538  {
4539  err = -2;
4540 
4541  if (sing_handler)
4542  {
4543  sing_handler (rcond);
4544  mattype.mark_as_rectangular ();
4545  }
4546  else
4547  gripe_singular_matrix (rcond);
4548  }
4549  }
4550  else
4551  rcond = 1.;
4552 
4553  if (err == 0)
4554  {
4555  retval = ComplexMatrix (b);
4556  Complex *result = retval.fortran_vec ();
4557 
4558  octave_idx_type b_nc = b.cols ();
4559 
4560  char job = 'N';
4561  F77_XFCN (zgbtrs, ZGBTRS,
4562  (F77_CONST_CHAR_ARG2 (&job, 1),
4563  nr, n_lower, n_upper, b_nc, tmp_data,
4564  ldm, pipvt, result, b.rows (), err
4565  F77_CHAR_ARG_LEN (1)));
4566  }
4567  }
4568  }
4569  else if (typ != MatrixType::Banded_Hermitian)
4570  (*current_liboctave_error_handler) ("incorrect matrix type");
4571  }
4572 
4573  return retval;
4574 }
4575 
4578  octave_idx_type& err, double& rcond,
4579  solve_singularity_handler sing_handler,
4580  bool calc_cond) const
4581 {
4582  SparseComplexMatrix retval;
4583 
4584  octave_idx_type nr = rows ();
4585  octave_idx_type nc = cols ();
4586  err = 0;
4587 
4588  if (nr != nc || nr != b.rows ())
4590  ("matrix dimension mismatch solution of linear equations");
4591  else if (nr == 0 || b.cols () == 0)
4592  retval = SparseComplexMatrix (nc, b.cols ());
4593  else
4594  {
4595  // Print spparms("spumoni") info if requested
4596  volatile int typ = mattype.type ();
4597  mattype.info ();
4598 
4599  if (typ == MatrixType::Banded_Hermitian)
4600  {
4601  octave_idx_type n_lower = mattype.nlower ();
4602  octave_idx_type ldm = n_lower + 1;
4603 
4604  ComplexMatrix m_band (ldm, nc);
4605  Complex *tmp_data = m_band.fortran_vec ();
4606 
4607  if (! mattype.is_dense ())
4608  {
4609  octave_idx_type ii = 0;
4610 
4611  for (octave_idx_type j = 0; j < ldm; j++)
4612  for (octave_idx_type i = 0; i < nc; i++)
4613  tmp_data[ii++] = 0.;
4614  }
4615 
4616  for (octave_idx_type j = 0; j < nc; j++)
4617  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4618  {
4619  octave_idx_type ri = ridx (i);
4620  if (ri >= j)
4621  m_band(ri - j, j) = data (i);
4622  }
4623 
4624  // Calculate the norm of the matrix, for later use.
4625  double anorm;
4626  if (calc_cond)
4627  anorm = m_band.abs ().sum ().row (0).max ();
4628 
4629  char job = 'L';
4630  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4631  nr, n_lower, tmp_data, ldm, err
4632  F77_CHAR_ARG_LEN (1)));
4633 
4634  if (err != 0)
4635  {
4636  rcond = 0.0;
4637  mattype.mark_as_unsymmetric ();
4638  typ = MatrixType::Banded;
4639  err = 0;
4640  }
4641  else
4642  {
4643  if (calc_cond)
4644  {
4645  Array<Complex> z (dim_vector (2 * nr, 1));
4646  Complex *pz = z.fortran_vec ();
4647  Array<double> iz (dim_vector (nr, 1));
4648  double *piz = iz.fortran_vec ();
4649 
4650  F77_XFCN (zpbcon, ZPBCON,
4651  (F77_CONST_CHAR_ARG2 (&job, 1),
4652  nr, n_lower, tmp_data, ldm,
4653  anorm, rcond, pz, piz, err
4654  F77_CHAR_ARG_LEN (1)));
4655 
4656  if (err != 0)
4657  err = -2;
4658 
4659  volatile double rcond_plus_one = rcond + 1.0;
4660 
4661  if (rcond_plus_one == 1.0 || xisnan (rcond))
4662  {
4663  err = -2;
4664 
4665  if (sing_handler)
4666  {
4667  sing_handler (rcond);
4668  mattype.mark_as_rectangular ();
4669  }
4670  else
4671  gripe_singular_matrix (rcond);
4672  }
4673  }
4674  else
4675  rcond = 1.0;
4676 
4677  if (err == 0)
4678  {
4679  octave_idx_type b_nr = b.rows ();
4680  octave_idx_type b_nc = b.cols ();
4681  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4682 
4683  // Take a first guess that the number of nonzero terms
4684  // will be as many as in b
4685  volatile octave_idx_type x_nz = b.nnz ();
4686  volatile octave_idx_type ii = 0;
4687  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4688 
4689  retval.xcidx (0) = 0;
4690  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4691  {
4692  for (octave_idx_type i = 0; i < b_nr; i++)
4693  Bx[i] = b.elem (i, j);
4694 
4695  F77_XFCN (zpbtrs, ZPBTRS,
4696  (F77_CONST_CHAR_ARG2 (&job, 1),
4697  nr, n_lower, 1, tmp_data,
4698  ldm, Bx, b_nr, err
4699  F77_CHAR_ARG_LEN (1)));
4700 
4701  if (err != 0)
4702  {
4703  (*current_liboctave_error_handler)
4704  ("SparseComplexMatrix::solve solve failed");
4705  err = -1;
4706  break;
4707  }
4708 
4709  for (octave_idx_type i = 0; i < b_nr; i++)
4710  {
4711  Complex tmp = Bx[i];
4712  if (tmp != 0.0)
4713  {
4714  if (ii == x_nz)
4715  {
4716  // Resize the sparse matrix
4717  octave_idx_type sz = x_nz *
4718  (b_nc - j) / b_nc;
4719  sz = (sz > 10 ? sz : 10) + x_nz;
4720  retval.change_capacity (sz);
4721  x_nz = sz;
4722  }
4723  retval.xdata (ii) = tmp;
4724  retval.xridx (ii++) = i;
4725  }
4726  }
4727  retval.xcidx (j+1) = ii;
4728  }
4729 
4730  retval.maybe_compress ();
4731  }
4732  }
4733  }
4734 
4735  if (typ == MatrixType::Banded)
4736  {
4737  // Create the storage for the banded form of the sparse matrix
4738  octave_idx_type n_upper = mattype.nupper ();
4739  octave_idx_type n_lower = mattype.nlower ();
4740  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4741 
4742  ComplexMatrix m_band (ldm, nc);
4743  Complex *tmp_data = m_band.fortran_vec ();
4744 
4745  if (! mattype.is_dense ())
4746  {
4747  octave_idx_type ii = 0;
4748 
4749  for (octave_idx_type j = 0; j < ldm; j++)
4750  for (octave_idx_type i = 0; i < nc; i++)
4751  tmp_data[ii++] = 0.;
4752  }
4753 
4754  for (octave_idx_type j = 0; j < nc; j++)
4755  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4756  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4757 
4758  // Calculate the norm of the matrix, for later use.
4759  double anorm = 0.0;
4760  if (calc_cond)
4761  {
4762  for (octave_idx_type j = 0; j < nr; j++)
4763  {
4764  double atmp = 0.;
4765  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4766  atmp += std::abs (data (i));
4767  if (atmp > anorm)
4768  anorm = atmp;
4769  }
4770  }
4771 
4772  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4773  octave_idx_type *pipvt = ipvt.fortran_vec ();
4774 
4775  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
4776  ldm, pipvt, err));
4777 
4778  if (err != 0)
4779  {
4780  rcond = 0.0;
4781  err = -2;
4782 
4783  if (sing_handler)
4784  {
4785  sing_handler (rcond);
4786  mattype.mark_as_rectangular ();
4787  }
4788  else
4790  }
4791  else
4792  {
4793  if (calc_cond)
4794  {
4795  char job = '1';
4796  Array<Complex> z (dim_vector (2 * nr, 1));
4797  Complex *pz = z.fortran_vec ();
4798  Array<double> iz (dim_vector (nr, 1));
4799  double *piz = iz.fortran_vec ();
4800 
4801  F77_XFCN (zgbcon, ZGBCON,
4802  (F77_CONST_CHAR_ARG2 (&job, 1),
4803  nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4804  anorm, rcond, pz, piz, err
4805  F77_CHAR_ARG_LEN (1)));
4806 
4807  if (err != 0)
4808  err = -2;
4809 
4810  volatile double rcond_plus_one = rcond + 1.0;
4811 
4812  if (rcond_plus_one == 1.0 || xisnan (rcond))
4813  {
4814  err = -2;
4815 
4816  if (sing_handler)
4817  {
4818  sing_handler (rcond);
4819  mattype.mark_as_rectangular ();
4820  }
4821  else
4822  gripe_singular_matrix (rcond);
4823  }
4824  }
4825  else
4826  rcond = 1.;
4827 
4828  if (err == 0)
4829  {
4830  char job = 'N';
4831  volatile octave_idx_type x_nz = b.nnz ();
4832  octave_idx_type b_nc = b.cols ();
4833  retval = SparseComplexMatrix (nr, b_nc, x_nz);
4834  retval.xcidx (0) = 0;
4835  volatile octave_idx_type ii = 0;
4836 
4837  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4838 
4839  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4840  {
4841  for (octave_idx_type i = 0; i < nr; i++)
4842  work[i] = 0.;
4843  for (octave_idx_type i = b.cidx (j);
4844  i < b.cidx (j+1); i++)
4845  work[b.ridx (i)] = b.data (i);
4846 
4847  F77_XFCN (zgbtrs, ZGBTRS,
4848  (F77_CONST_CHAR_ARG2 (&job, 1),
4849  nr, n_lower, n_upper, 1, tmp_data,
4850  ldm, pipvt, work, b.rows (), err
4851  F77_CHAR_ARG_LEN (1)));
4852 
4853  // Count nonzeros in work vector and adjust
4854  // space in retval if needed
4855  octave_idx_type new_nnz = 0;
4856  for (octave_idx_type i = 0; i < nr; i++)
4857  if (work[i] != 0.)
4858  new_nnz++;
4859 
4860  if (ii + new_nnz > x_nz)
4861  {
4862  // Resize the sparse matrix
4863  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4864  retval.change_capacity (sz);
4865  x_nz = sz;
4866  }
4867 
4868  for (octave_idx_type i = 0; i < nr; i++)
4869  if (work[i] != 0.)
4870  {
4871  retval.xridx (ii) = i;
4872  retval.xdata (ii++) = work[i];
4873  }
4874  retval.xcidx (j+1) = ii;
4875  }
4876 
4877  retval.maybe_compress ();
4878  }
4879  }
4880  }
4881  else if (typ != MatrixType::Banded_Hermitian)
4882  (*current_liboctave_error_handler) ("incorrect matrix type");
4883  }
4884 
4885  return retval;
4886 }
4887 
4890  octave_idx_type& err, double& rcond,
4891  solve_singularity_handler sing_handler,
4892  bool calc_cond) const
4893 {
4894  ComplexMatrix retval;
4895 
4896  octave_idx_type nr = rows ();
4897  octave_idx_type nc = cols ();
4898  err = 0;
4899 
4900  if (nr != nc || nr != b.rows ())
4902  ("matrix dimension mismatch solution of linear equations");
4903  else if (nr == 0 || b.cols () == 0)
4904  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4905  else
4906  {
4907  // Print spparms("spumoni") info if requested
4908  volatile int typ = mattype.type ();
4909  mattype.info ();
4910 
4911  if (typ == MatrixType::Banded_Hermitian)
4912  {
4913  octave_idx_type n_lower = mattype.nlower ();
4914  octave_idx_type ldm = n_lower + 1;
4915 
4916  ComplexMatrix m_band (ldm, nc);
4917  Complex *tmp_data = m_band.fortran_vec ();
4918 
4919  if (! mattype.is_dense ())
4920  {
4921  octave_idx_type ii = 0;
4922 
4923  for (octave_idx_type j = 0; j < ldm; j++)
4924  for (octave_idx_type i = 0; i < nc; i++)
4925  tmp_data[ii++] = 0.;
4926  }
4927 
4928  for (octave_idx_type j = 0; j < nc; j++)
4929  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4930  {
4931  octave_idx_type ri = ridx (i);
4932  if (ri >= j)
4933  m_band(ri - j, j) = data (i);
4934  }
4935 
4936  // Calculate the norm of the matrix, for later use.
4937  double anorm;
4938  if (calc_cond)
4939  anorm = m_band.abs ().sum ().row (0).max ();
4940 
4941  char job = 'L';
4942  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4943  nr, n_lower, tmp_data, ldm, err
4944  F77_CHAR_ARG_LEN (1)));
4945 
4946  if (err != 0)
4947  {
4948  // Matrix is not positive definite!! Fall through to
4949  // unsymmetric banded solver.
4950  rcond = 0.0;
4951  mattype.mark_as_unsymmetric ();
4952  typ = MatrixType::Banded;
4953  err = 0;
4954  }
4955  else
4956  {
4957  if (calc_cond)
4958  {
4959  Array<Complex> z (dim_vector (2 * nr, 1));
4960  Complex *pz = z.fortran_vec ();
4961  Array<double> iz (dim_vector (nr, 1));
4962  double *piz = iz.fortran_vec ();
4963 
4964  F77_XFCN (zpbcon, ZPBCON,
4965  (F77_CONST_CHAR_ARG2 (&job, 1),
4966  nr, n_lower, tmp_data, ldm,
4967  anorm, rcond, pz, piz, err
4968  F77_CHAR_ARG_LEN (1)));
4969 
4970  if (err != 0)
4971  err = -2;
4972 
4973  volatile double rcond_plus_one = rcond + 1.0;
4974 
4975  if (rcond_plus_one == 1.0 || xisnan (rcond))
4976  {
4977  err = -2;
4978 
4979  if (sing_handler)
4980  {
4981  sing_handler (rcond);
4982  mattype.mark_as_rectangular ();
4983  }
4984  else
4985  gripe_singular_matrix (rcond);
4986  }
4987  }
4988  else
4989  rcond = 1.0;
4990 
4991  if (err == 0)
4992  {
4993  octave_idx_type b_nr = b.rows ();
4994  octave_idx_type b_nc = b.cols ();
4995  retval = ComplexMatrix (b);
4996  Complex *result = retval.fortran_vec ();
4997 
4998  F77_XFCN (zpbtrs, ZPBTRS,
4999  (F77_CONST_CHAR_ARG2 (&job, 1),
5000  nr, n_lower, b_nc, tmp_data,
5001  ldm, result, b_nr, err
5002  F77_CHAR_ARG_LEN (1)));
5003 
5004  if (err != 0)
5005  {
5006  (*current_liboctave_error_handler)
5007  ("SparseComplexMatrix::solve solve failed");
5008  err = -1;
5009  }
5010  }
5011  }
5012  }
5013 
5014  if (typ == MatrixType::Banded)
5015  {
5016  // Create the storage for the banded form of the sparse matrix
5017  octave_idx_type n_upper = mattype.nupper ();
5018  octave_idx_type n_lower = mattype.nlower ();
5019  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
5020 
5021  ComplexMatrix m_band (ldm, nc);
5022  Complex *tmp_data = m_band.fortran_vec ();
5023 
5024  if (! mattype.is_dense ())
5025  {
5026  octave_idx_type ii = 0;
5027 
5028  for (octave_idx_type j = 0; j < ldm; j++)
5029  for (octave_idx_type i = 0; i < nc; i++)
5030  tmp_data[ii++] = 0.;
5031  }
5032 
5033  for (octave_idx_type j = 0; j < nc; j++)
5034  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5035  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5036 
5037  // Calculate the norm of the matrix, for later use.
5038  double anorm = 0.0;
5039  if (calc_cond)
5040  {
5041  for (octave_idx_type j = 0; j < nr; j++)
5042  {
5043  double atmp = 0.;
5044  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5045  atmp += std::abs (data (i));
5046  if (atmp > anorm)
5047  anorm = atmp;
5048  }
5049  }
5050 
5051  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
5052  octave_idx_type *pipvt = ipvt.fortran_vec ();
5053 
5054  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5055  ldm, pipvt, err));
5056 
5057  if (err != 0)
5058  {
5059  err = -2;
5060  rcond = 0.0;
5061 
5062  if (sing_handler)
5063  {
5064  sing_handler (rcond);
5065  mattype.mark_as_rectangular ();
5066  }
5067  else
5069  }
5070  else
5071  {
5072  if (calc_cond)
5073  {
5074  char job = '1';
5075  Array<Complex> z (dim_vector (2 * nr, 1));
5076  Complex *pz = z.fortran_vec ();
5077  Array<double> iz (dim_vector (nr, 1));
5078  double *piz = iz.fortran_vec ();
5079 
5080  F77_XFCN (zgbcon, ZGBCON,
5081  (F77_CONST_CHAR_ARG2 (&job, 1),
5082  nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5083  anorm, rcond, pz, piz, err
5084  F77_CHAR_ARG_LEN (1)));
5085 
5086  if (err != 0)
5087  err = -2;
5088 
5089  volatile double rcond_plus_one = rcond + 1.0;
5090 
5091  if (rcond_plus_one == 1.0 || xisnan (rcond))
5092  {
5093  err = -2;
5094 
5095  if (sing_handler)
5096  {
5097  sing_handler (rcond);
5098  mattype.mark_as_rectangular ();
5099  }
5100  else
5101  gripe_singular_matrix (rcond);
5102  }
5103  }
5104  else
5105  rcond = 1.;
5106 
5107  if (err == 0)
5108  {
5109  char job = 'N';
5110  octave_idx_type b_nc = b.cols ();
5111  retval = ComplexMatrix (b);
5112  Complex *result = retval.fortran_vec ();
5113 
5114  F77_XFCN (zgbtrs, ZGBTRS,
5115  (F77_CONST_CHAR_ARG2 (&job, 1),
5116  nr, n_lower, n_upper, b_nc, tmp_data,
5117  ldm, pipvt, result, b.rows (), err
5118  F77_CHAR_ARG_LEN (1)));
5119  }
5120  }
5121  }
5122  else if (typ != MatrixType::Banded_Hermitian)
5123  (*current_liboctave_error_handler) ("incorrect matrix type");
5124  }
5125 
5126  return retval;
5127 }
5128 
5131  octave_idx_type& err, double& rcond,
5132  solve_singularity_handler sing_handler,
5133  bool calc_cond) const
5134 {
5135  SparseComplexMatrix retval;
5136 
5137  octave_idx_type nr = rows ();
5138  octave_idx_type nc = cols ();
5139  err = 0;
5140 
5141  if (nr != nc || nr != b.rows ())
5143  ("matrix dimension mismatch solution of linear equations");
5144  else if (nr == 0 || b.cols () == 0)
5145  retval = SparseComplexMatrix (nc, b.cols ());
5146  else
5147  {
5148  // Print spparms("spumoni") info if requested
5149  volatile int typ = mattype.type ();
5150  mattype.info ();
5151 
5152  if (typ == MatrixType::Banded_Hermitian)
5153  {
5154  octave_idx_type n_lower = mattype.nlower ();
5155  octave_idx_type ldm = n_lower + 1;
5156 
5157  ComplexMatrix m_band (ldm, nc);
5158  Complex *tmp_data = m_band.fortran_vec ();
5159 
5160  if (! mattype.is_dense ())
5161  {
5162  octave_idx_type ii = 0;
5163 
5164  for (octave_idx_type j = 0; j < ldm; j++)
5165  for (octave_idx_type i = 0; i < nc; i++)
5166  tmp_data[ii++] = 0.;
5167  }
5168 
5169  for (octave_idx_type j = 0; j < nc; j++)
5170  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5171  {
5172  octave_idx_type ri = ridx (i);
5173  if (ri >= j)
5174  m_band(ri - j, j) = data (i);
5175  }
5176 
5177  // Calculate the norm of the matrix, for later use.
5178  double anorm;
5179  if (calc_cond)
5180  anorm = m_band.abs ().sum ().row (0).max ();
5181 
5182  char job = 'L';
5183  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5184  nr, n_lower, tmp_data, ldm, err
5185  F77_CHAR_ARG_LEN (1)));
5186 
5187  if (err != 0)
5188  {
5189  // Matrix is not positive definite!! Fall through to
5190  // unsymmetric banded solver.
5191  mattype.mark_as_unsymmetric ();
5192  typ = MatrixType::Banded;
5193 
5194  rcond = 0.0;
5195  err = 0;
5196  }
5197  else
5198  {
5199  if (calc_cond)
5200  {
5201  Array<Complex> z (dim_vector (2 * nr, 1));
5202  Complex *pz = z.fortran_vec ();
5203  Array<double> iz (dim_vector (nr, 1));
5204  double *piz = iz.fortran_vec ();
5205 
5206  F77_XFCN (zpbcon, ZPBCON,
5207  (F77_CONST_CHAR_ARG2 (&job, 1),
5208  nr, n_lower, tmp_data, ldm,
5209  anorm, rcond, pz, piz, err
5210  F77_CHAR_ARG_LEN (1)));
5211 
5212  if (err != 0)
5213  err = -2;
5214 
5215  volatile double rcond_plus_one = rcond + 1.0;
5216 
5217  if (rcond_plus_one == 1.0 || xisnan (rcond))
5218  {
5219  err = -2;
5220 
5221  if (sing_handler)
5222  {
5223  sing_handler (rcond);
5224  mattype.mark_as_rectangular ();
5225  }
5226  else
5227  gripe_singular_matrix (rcond);
5228  }
5229  }
5230  else
5231  rcond = 1.0;
5232 
5233  if (err == 0)
5234  {
5235  octave_idx_type b_nr = b.rows ();
5236  octave_idx_type b_nc = b.cols ();
5237  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
5238 
5239  // Take a first guess that the number of nonzero terms
5240  // will be as many as in b
5241  volatile octave_idx_type x_nz = b.nnz ();
5242  volatile octave_idx_type ii = 0;
5243  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5244 
5245  retval.xcidx (0) = 0;
5246  for (volatile octave_idx_type j = 0; j < b_nc; j++)
5247  {
5248 
5249  for (octave_idx_type i = 0; i < b_nr; i++)
5250  Bx[i] = b(i,j);
5251 
5252  F77_XFCN (zpbtrs, ZPBTRS,
5253  (F77_CONST_CHAR_ARG2 (&job, 1),
5254  nr, n_lower, 1, tmp_data,
5255  ldm, Bx, b_nr, err
5256  F77_CHAR_ARG_LEN (1)));
5257 
5258  if (err != 0)
5259  {
5260  (*current_liboctave_error_handler)
5261  ("SparseMatrix::solve solve failed");
5262  err = -1;
5263  break;
5264  }
5265 
5266  // Count nonzeros in work vector and adjust
5267  // space in retval if needed
5268  octave_idx_type new_nnz = 0;
5269  for (octave_idx_type i = 0; i < nr; i++)
5270  if (Bx[i] != 0.)
5271  new_nnz++;
5272 
5273  if (ii + new_nnz > x_nz)
5274  {
5275  // Resize the sparse matrix
5276  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5277  retval.change_capacity (sz);
5278  x_nz = sz;
5279  }
5280 
5281  for (octave_idx_type i = 0; i < nr; i++)
5282  if (Bx[i] != 0.)
5283  {
5284  retval.xridx (ii) = i;
5285  retval.xdata (ii++) = Bx[i];
5286  }
5287 
5288  retval.xcidx (j+1) = ii;
5289  }
5290 
5291  retval.maybe_compress ();
5292  }
5293  }
5294  }
5295 
5296  if (typ == MatrixType::Banded)
5297  {
5298  // Create the storage for the banded form of the sparse matrix
5299  octave_idx_type n_upper = mattype.nupper ();
5300  octave_idx_type n_lower = mattype.nlower ();
5301  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
5302 
5303  ComplexMatrix m_band (ldm, nc);
5304  Complex *tmp_data = m_band.fortran_vec ();
5305 
5306  if (! mattype.is_dense ())
5307  {
5308  octave_idx_type ii = 0;
5309 
5310  for (octave_idx_type j = 0; j < ldm; j++)
5311  for (octave_idx_type i = 0; i < nc; i++)
5312  tmp_data[ii++] = 0.;
5313  }
5314 
5315  for (octave_idx_type j = 0; j < nc; j++)
5316  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5317  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5318 
5319  // Calculate the norm of the matrix, for later use.
5320  double anorm = 0.0;
5321  if (calc_cond)
5322  {
5323  for (octave_idx_type j = 0; j < nr; j++)
5324  {
5325  double atmp = 0.;
5326  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5327  atmp += std::abs (data (i));
5328  if (atmp > anorm)
5329  anorm = atmp;
5330  }
5331  }
5332 
5333  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
5334  octave_idx_type *pipvt = ipvt.fortran_vec ();
5335 
5336  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5337  ldm, pipvt, err));
5338 
5339  if (err != 0)
5340  {
5341  err = -2;
5342  rcond = 0.0;
5343 
5344  if (sing_handler)
5345  {
5346  sing_handler (rcond);
5347  mattype.mark_as_rectangular ();
5348  }
5349  else
5351  }
5352  else
5353  {
5354  if (calc_cond)
5355  {
5356  char job = '1';
5357  Array<Complex> z (dim_vector (2 * nr, 1));
5358  Complex *pz = z.fortran_vec ();
5359  Array<double> iz (dim_vector (nr, 1));
5360  double *piz = iz.fortran_vec ();
5361 
5362  F77_XFCN (zgbcon, ZGBCON,
5363  (F77_CONST_CHAR_ARG2 (&job, 1),
5364  nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5365  anorm, rcond, pz, piz, err
5366  F77_CHAR_ARG_LEN (1)));
5367 
5368  if (err != 0)
5369  err = -2;
5370 
5371  volatile double rcond_plus_one = rcond + 1.0;
5372 
5373  if (rcond_plus_one == 1.0 || xisnan (rcond))
5374  {
5375  err = -2;
5376 
5377  if (sing_handler)
5378  {
5379  sing_handler (rcond);
5380  mattype.mark_as_rectangular ();
5381  }
5382  else
5383  gripe_singular_matrix (rcond);
5384  }
5385  }
5386  else
5387  rcond = 1.;
5388 
5389  if (err == 0)
5390  {
5391  char job = 'N';
5392  volatile octave_idx_type x_nz = b.nnz ();
5393  octave_idx_type b_nc = b.cols ();
5394  retval = SparseComplexMatrix (nr, b_nc, x_nz);
5395  retval.xcidx (0) = 0;
5396  volatile octave_idx_type ii = 0;
5397 
5398  OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
5399 
5400  for (volatile octave_idx_type j = 0; j < b_nc; j++)
5401  {
5402  for (octave_idx_type i = 0; i < nr; i++)
5403  Bx[i] = 0.;
5404 
5405  for (octave_idx_type i = b.cidx (j);
5406  i < b.cidx (j+1); i++)
5407  Bx[b.ridx (i)] = b.data (i);
5408 
5409  F77_XFCN (zgbtrs, ZGBTRS,
5410  (F77_CONST_CHAR_ARG2 (&job, 1),
5411  nr, n_lower, n_upper, 1, tmp_data,
5412  ldm, pipvt, Bx, b.rows (), err
5413  F77_CHAR_ARG_LEN (1)));
5414 
5415  // Count nonzeros in work vector and adjust
5416  // space in retval if needed
5417  octave_idx_type new_nnz = 0;
5418  for (octave_idx_type i = 0; i < nr; i++)
5419  if (Bx[i] != 0.)
5420  new_nnz++;
5421 
5422  if (ii + new_nnz > x_nz)
5423  {
5424  // Resize the sparse matrix
5425  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5426  retval.change_capacity (sz);
5427  x_nz = sz;
5428  }
5429 
5430  for (octave_idx_type i = 0; i < nr; i++)
5431  if (Bx[i] != 0.)
5432  {
5433  retval.xridx (ii) = i;
5434  retval.xdata (ii++) = Bx[i];
5435  }
5436  retval.xcidx (j+1) = ii;
5437  }
5438 
5439  retval.maybe_compress ();
5440  }
5441  }
5442  }
5443  else if (typ != MatrixType::Banded_Hermitian)
5444  (*current_liboctave_error_handler) ("incorrect matrix type");
5445  }
5446 
5447  return retval;
5448 }
5449 
5450 void *
5452  Matrix &Control, Matrix &Info,
5453  solve_singularity_handler sing_handler,
5454  bool calc_cond) const
5455 {
5456  // The return values
5457  void *Numeric = 0;
5458  err = 0;
5459 
5460 #ifdef HAVE_UMFPACK
5461  // Setup the control parameters
5462  Control = Matrix (UMFPACK_CONTROL, 1);
5463  double *control = Control.fortran_vec ();
5464  UMFPACK_ZNAME (defaults) (control);
5465 
5466  double tmp = octave_sparse_params::get_key ("spumoni");
5467  if (!xisnan (tmp))
5468  Control (UMFPACK_PRL) = tmp;
5469  tmp = octave_sparse_params::get_key ("piv_tol");
5470  if (!xisnan (tmp))
5471  {
5472  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5473  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5474  }
5475 
5476  // Set whether we are allowed to modify Q or not
5477  tmp = octave_sparse_params::get_key ("autoamd");
5478  if (!xisnan (tmp))
5479  Control (UMFPACK_FIXQ) = tmp;
5480 
5481  UMFPACK_ZNAME (report_control) (control);
5482 
5483  const octave_idx_type *Ap = cidx ();
5484  const octave_idx_type *Ai = ridx ();
5485  const Complex *Ax = data ();
5486  octave_idx_type nr = rows ();
5487  octave_idx_type nc = cols ();
5488 
5489  UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
5490  reinterpret_cast<const double *> (Ax),
5491  0, 1, control);
5492 
5493  void *Symbolic;
5494  Info = Matrix (1, UMFPACK_INFO);
5495  double *info = Info.fortran_vec ();
5496  int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
5497  reinterpret_cast<const double *> (Ax),
5498  0, 0, &Symbolic, control, info);
5499 
5500  if (status < 0)
5501  {
5502  (*current_liboctave_error_handler)
5503  ("SparseComplexMatrix::solve symbolic factorization failed");
5504  err = -1;
5505 
5506  UMFPACK_ZNAME (report_status) (control, status);
5507  UMFPACK_ZNAME (report_info) (control, info);
5508 
5509  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5510  }
5511  else
5512  {
5513  UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
5514 
5515  status = UMFPACK_ZNAME (numeric) (Ap, Ai,
5516  reinterpret_cast<const double *> (Ax),
5517  0, Symbolic, &Numeric, control, info);
5518  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5519 
5520  if (calc_cond)
5521  rcond = Info (UMFPACK_RCOND);
5522  else
5523  rcond = 1.;
5524  volatile double rcond_plus_one = rcond + 1.0;
5525 
5526  if (status == UMFPACK_WARNING_singular_matrix
5527  || rcond_plus_one == 1.0 || xisnan (rcond))
5528  {
5529  UMFPACK_ZNAME (report_numeric) (Numeric, control);
5530 
5531  err = -2;
5532 
5533  if (sing_handler)
5534  sing_handler (rcond);
5535  else
5536  gripe_singular_matrix (rcond);
5537  }
5538  else if (status < 0)
5539  {
5540  (*current_liboctave_error_handler)
5541  ("SparseComplexMatrix::solve numeric factorization failed");
5542 
5543  UMFPACK_ZNAME (report_status) (control, status);
5544  UMFPACK_ZNAME (report_info) (control, info);
5545 
5546  err = -1;
5547  }
5548  else
5549  {
5550  UMFPACK_ZNAME (report_numeric) (Numeric, control);
5551  }
5552  }
5553 
5554  if (err != 0)
5555  UMFPACK_ZNAME (free_numeric) (&Numeric);
5556 #else
5557  (*current_liboctave_error_handler) ("UMFPACK not installed");
5558 #endif
5559 
5560  return Numeric;
5561 }
5562 
5565  octave_idx_type& err, double& rcond,
5566  solve_singularity_handler sing_handler,
5567  bool calc_cond) const
5568 {
5569  ComplexMatrix retval;
5570 
5571  octave_idx_type nr = rows ();
5572  octave_idx_type nc = cols ();
5573  err = 0;
5574 
5575  if (nr != nc || nr != b.rows ())
5577  ("matrix dimension mismatch solution of linear equations");
5578  else if (nr == 0 || b.cols () == 0)
5579  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
5580  else
5581  {
5582  // Print spparms("spumoni") info if requested
5583  volatile int typ = mattype.type ();
5584  mattype.info ();
5585 
5586  if (typ == MatrixType::Hermitian)
5587  {
5588 #ifdef HAVE_CHOLMOD
5589  cholmod_common Common;
5590  cholmod_common *cm = &Common;
5591 
5592  // Setup initial parameters
5593  CHOLMOD_NAME(start) (cm);
5594  cm->prefer_zomplex = false;
5595 
5596  double spu = octave_sparse_params::get_key ("spumoni");
5597  if (spu == 0.)
5598  {
5599  cm->print = -1;
5600  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
5601  }
5602  else
5603  {
5604  cm->print = static_cast<int> (spu) + 2;
5605  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
5606  }
5607 
5608  cm->error_handler = &SparseCholError;
5609  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5610  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5611 
5612  cm->final_ll = true;
5613 
5614  cholmod_sparse Astore;
5615  cholmod_sparse *A = &Astore;
5616  double dummy;
5617  A->nrow = nr;
5618  A->ncol = nc;
5619 
5620  A->p = cidx ();
5621  A->i = ridx ();
5622  A->nzmax = nnz ();
5623  A->packed = true;
5624  A->sorted = true;
5625  A->nz = 0;
5626 #ifdef USE_64_BIT_IDX_T
5627  A->itype = CHOLMOD_LONG;
5628 #else
5629  A->itype = CHOLMOD_INT;
5630 #endif
5631  A->dtype = CHOLMOD_DOUBLE;
5632  A->stype = 1;
5633  A->xtype = CHOLMOD_COMPLEX;
5634 
5635  if (nr < 1)
5636  A->x = &dummy;
5637  else
5638  A->x = data ();
5639 
5640  cholmod_dense Bstore;
5641  cholmod_dense *B = &Bstore;
5642  B->nrow = b.rows ();
5643  B->ncol = b.cols ();
5644  B->d = B->nrow;
5645  B->nzmax = B->nrow * B->ncol;
5646  B->dtype = CHOLMOD_DOUBLE;
5647  B->xtype = CHOLMOD_REAL;
5648  if (nc < 1 || b.cols () < 1)
5649  B->x = &dummy;
5650  else
5651  // We won't alter it, honest :-)
5652  B->x = const_cast<double *>(b.fortran_vec ());
5653 
5654  cholmod_factor *L;
5656  L = CHOLMOD_NAME(analyze) (A, cm);
5657  CHOLMOD_NAME(factorize) (A, L, cm);
5658  if (calc_cond)
5659  rcond = CHOLMOD_NAME(rcond)(L, cm);
5660  else
5661  rcond = 1.;
5663 
5664  if (rcond == 0.0)
5665  {
5666  // Either its indefinite or singular. Try UMFPACK
5667  mattype.mark_as_unsymmetric ();
5668  typ = MatrixType::Full;
5669  }
5670  else
5671  {
5672  volatile double rcond_plus_one = rcond + 1.0;
5673 
5674  if (rcond_plus_one == 1.0 || xisnan (rcond))
5675  {
5676  err = -2;
5677 
5678  if (sing_handler)
5679  {
5680  sing_handler (rcond);
5681  mattype.mark_as_rectangular ();
5682  }
5683  else
5684  gripe_singular_matrix (rcond);
5685 
5686  return retval;
5687  }
5688 
5689  cholmod_dense *X;
5691  X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5693 
5694  retval.resize (b.rows (), b.cols ());
5695  for (octave_idx_type j = 0; j < b.cols (); j++)
5696  {
5697  octave_idx_type jr = j * b.rows ();
5698  for (octave_idx_type i = 0; i < b.rows (); i++)
5699  retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i];
5700  }
5701 
5703  CHOLMOD_NAME(free_dense) (&X, cm);
5704  CHOLMOD_NAME(free_factor) (&L, cm);
5705  CHOLMOD_NAME(finish) (cm);
5706  static char tmp[] = " ";
5707  CHOLMOD_NAME(print_common) (tmp, cm);
5709  }
5710 #else
5711  (*current_liboctave_warning_with_id_handler)
5712  ("Octave:missing-dependency", "CHOLMOD not installed");
5713 
5714  mattype.mark_as_unsymmetric ();
5715  typ = MatrixType::Full;
5716 #endif
5717  }
5718 
5719  if (typ == MatrixType::Full)
5720  {
5721 #ifdef HAVE_UMFPACK
5722  Matrix Control, Info;
5723  void *Numeric = factorize (err, rcond, Control, Info,
5724  sing_handler, calc_cond);
5725 
5726  if (err == 0)
5727  {
5728  octave_idx_type b_nr = b.rows ();
5729  octave_idx_type b_nc = b.cols ();
5730  int status = 0;
5731  double *control = Control.fortran_vec ();
5732  double *info = Info.fortran_vec ();
5733  const octave_idx_type *Ap = cidx ();
5734  const octave_idx_type *Ai = ridx ();
5735  const Complex *Ax = data ();
5736 #ifdef UMFPACK_SEPARATE_SPLIT
5737  const double *Bx = b.fortran_vec ();
5738  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
5739  for (octave_idx_type i = 0; i < b_nr; i++)
5740  Bz[i] = 0.;
5741 #else
5742  OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
5743 #endif
5744  retval.resize (b_nr, b_nc);
5745  Complex *Xx = retval.fortran_vec ();
5746 
5747  for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
5748  {
5749 #ifdef UMFPACK_SEPARATE_SPLIT
5750  status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5751  Ai,
5752  reinterpret_cast<const double *> (Ax),
5753  0,
5754  reinterpret_cast<double *> (&Xx[iidx]),
5755  0,
5756  &Bx[iidx], Bz, Numeric,
5757  control, info);
5758 #else
5759  for (octave_idx_type i = 0; i < b_nr; i++)
5760  Bz[i] = b.elem (i, j);
5761 
5762  status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5763  Ai,
5764  reinterpret_cast<const double *> (Ax),
5765  0,
5766  reinterpret_cast<double *> (&Xx[iidx]),
5767  0,
5768  reinterpret_cast<const double *> (Bz),
5769  0, Numeric,
5770  control, info);
5771 #endif
5772 
5773  if (status < 0)
5774  {
5775  (*current_liboctave_error_handler)
5776  ("SparseComplexMatrix::solve solve failed");
5777 
5778  UMFPACK_ZNAME (report_status) (control, status);
5779 
5780  err = -1;
5781 
5782  break;
5783  }
5784  }
5785 
5786  UMFPACK_ZNAME (report_info) (control, info);
5787 
5788  UMFPACK_ZNAME (free_numeric) (&Numeric);
5789  }
5790  else
5791  mattype.mark_as_rectangular ();
5792 
5793 #else
5794  (*current_liboctave_error_handler) ("UMFPACK not installed");
5795 #endif
5796  }
5797  else if (typ != MatrixType::Hermitian)
5798  (*current_liboctave_error_handler) ("incorrect matrix type");
5799  }
5800 
5801  return retval;
5802 }
5803 
5806  octave_idx_type& err, double& rcond,
5807  solve_singularity_handler sing_handler,
5808  bool calc_cond) const
5809 {
5810  SparseComplexMatrix retval;
5811 
5812  octave_idx_type nr = rows ();
5813  octave_idx_type nc = cols ();
5814  err = 0;
5815 
5816  if (nr != nc || nr != b.rows ())
5818  ("matrix dimension mismatch solution of linear equations");
5819  else if (nr == 0 || b.cols () == 0)
5820  retval = SparseComplexMatrix (nc, b.cols ());
5821  else
5822  {
5823  // Print spparms("spumoni") info if requested
5824  volatile int typ = mattype.type ();
5825  mattype.info ();
5826 
5827  if (typ == MatrixType::Hermitian)
5828  {
5829 #ifdef HAVE_CHOLMOD
5830  cholmod_common Common;
5831  cholmod_common *cm = &Common;
5832 
5833  // Setup initial parameters
5834  CHOLMOD_NAME(start) (cm);
5835  cm->prefer_zomplex = false;
5836 
5837  double spu = octave_sparse_params::get_key ("spumoni");
5838  if (spu == 0.)
5839  {
5840  cm->print = -1;
5841  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
5842  }
5843  else
5844  {
5845  cm->print = static_cast<int> (spu) + 2;
5846  SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
5847  }
5848 
5849  cm->error_handler = &SparseCholError;
5850  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
5851  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
5852 
5853  cm->final_ll = true;
5854 
5855  cholmod_sparse Astore;
5856  cholmod_sparse *A = &Astore;
5857  double dummy;
5858  A->nrow = nr;
5859  A->ncol = nc;
5860 
5861  A->p = cidx ();
5862  A->i = ridx ();
5863  A->nzmax = nnz ();
5864  A->packed = true;
5865  A->sorted = true;
5866  A->nz = 0;
5867 #ifdef USE_64_BIT_IDX_T
5868  A->itype = CHOLMOD_LONG;
5869 #else
5870  A->itype = CHOLMOD_INT;
5871 #endif
5872  A->dtype = CHOLMOD_DOUBLE;
5873  A->stype = 1;
5874  A->xtype = CHOLMOD_COMPLEX;
5875 
5876  if (nr < 1)
5877  A->x = &dummy;
5878  else
5879  A->x = data ();
5880 
5881  cholmod_sparse Bstore;
5882  cholmod_sparse *B = &Bstore;
5883  B->nrow = b.rows ();
5884  B->ncol = b.cols ();
5885  B->p = b.cidx ();
5886  B->i = b.ridx ();
5887  B->nzmax = b.nnz ();
5888  B->packed = true;
5889  B->sorted = true;
5890  B->nz = 0;
5891 #ifdef USE_64_BIT_IDX_T
5892  B->itype = CHOLMOD_LONG;
5893 #else
5894  B->itype = CHOLMOD_INT;
5895 #endif
5896  B->dtype = CHOLMOD_DOUBLE;
5897  B->stype = 0;
5898  B->xtype = CHOLMOD_REAL;
5899 
5900  if (b.rows () < 1 || b.cols () < 1)
5901  B->x = &dummy;
5902  else
5903  B->x = b.data ();
5904 
5905  cholmod_factor *L;
5907  L = CHOLMOD_NAME(analyze) (A, cm);
5908  CHOLMOD_NAME(factorize) (A, L, cm);
5909  if (calc_cond)
5910  rcond = CHOLMOD_NAME(rcond)(L, cm);
5911  else
5912  rcond = 1.;
5914 
5915  if (rcond == 0.0)
5916  {
5917  // Either its indefinite or singular. Try UMFPACK
5918  mattype.mark_as_unsymmetric ();
5919  typ = MatrixType::Full;
5920  }
5921  else
5922  {
5923  volatile double rcond_plus_one = rcond + 1.0;
5924 
5925  if (rcond_plus_one == 1.0 || xisnan (rcond))
5926  {
5927  err = -2;
5928 
5929  if (sing_handler)
5930  {
5931  sing_handler (rcond);
5932  mattype.mark_as_rectangular ();
5933  }
5934  else
5935  gripe_singular_matrix (rcond);
5936 
5937  return retval;
5938  }
5939 
5940  cholmod_sparse *X;
5942  X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
5944 
5945  retval = SparseComplexMatrix
5946  (static_cast<octave_idx_type>(X->nrow),
5947  static_cast<octave_idx_type>(X->ncol),
5948  static_cast<octave_idx_type>(X->nzmax));
5949  for (octave_idx_type j = 0;
5950  j <= static_cast<octave_idx_type>(X->ncol); j++)
5951  retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
5952  for (octave_idx_type j = 0;
5953  j < static_cast<octave_idx_type>(X->nzmax); j++)
5954  {
5955  retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
5956  retval.xdata (j) = static_cast<Complex *>(X->x)[j];
5957  }
5958 
5960  CHOLMOD_NAME(free_sparse) (&X, cm);
5961  CHOLMOD_NAME(free_factor) (&L, cm);
5962  CHOLMOD_NAME(finish) (cm);
5963  static char tmp[] = " ";
5964  CHOLMOD_NAME(print_common) (tmp, cm);
5966  }
5967 #else
5968  (*current_liboctave_warning_with_id_handler)
5969  ("Octave:missing-dependency", "CHOLMOD not installed");
5970 
5971  mattype.mark_as_unsymmetric ();
5972  typ = MatrixType::Full;
5973 #endif
5974  }
5975 
5976  if (typ == MatrixType::Full)
5977  {
5978 #ifdef HAVE_UMFPACK
5979  Matrix Control, Info;
5980  void *Numeric = factorize (err, rcond, Control, Info,
5981  sing_handler, calc_cond);
5982 
5983  if (err == 0)
5984  {
5985  octave_idx_type b_nr = b.rows ();
5986