GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CSparse.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2017 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5 Copyright (C) 2010 VZLU Prague
6 
7 This file is part of Octave.
8 
9 Octave is free software; you can redistribute it and/or modify it
10 under the terms of the GNU General Public License as published by the
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13 
14 Octave is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along with Octave; see the file COPYING. If not, see
21 <http://www.gnu.org/licenses/>.
22 
23 */
24 
25 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28 
29 #include <cfloat>
30 
31 #include <iostream>
32 #include <vector>
33 
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 "lo-lapack-proto.h"
40 #include "mx-m-cs.h"
41 #include "mx-cs-m.h"
42 #include "mx-cm-s.h"
43 #include "mx-fcm-fs.h"
44 #include "mx-s-cm.h"
45 #include "mx-fs-fcm.h"
46 #include "oct-locbuf.h"
47 
48 #include "dDiagMatrix.h"
49 #include "CDiagMatrix.h"
50 #include "CSparse.h"
51 #include "boolSparse.h"
52 #include "dSparse.h"
53 #include "functor.h"
54 #include "oct-spparms.h"
55 #include "sparse-lu.h"
56 #include "oct-sparse.h"
57 #include "sparse-util.h"
58 #include "sparse-chol.h"
59 #include "sparse-qr.h"
60 
61 #include "Sparse-op-defs.h"
62 
63 #include "Sparse-diag-op-defs.h"
64 
65 #include "Sparse-perm-op-defs.h"
66 
67 // Define whether to use a basic QR solver or one that uses a Dulmange
68 // Mendelsohn factorization to seperate the problem into under-determined,
69 // well-determined and over-determined parts and solves them seperately
70 #if ! defined (USE_QRSOLVE)
71 # include "sparse-dmsolve.h"
72 #endif
73 
75  : MSparse<Complex> (a)
76 { }
77 
79  : MSparse<Complex> (a.rows (), a.cols (), a.nnz ())
80 {
81  octave_idx_type nc = cols ();
82  octave_idx_type nz = a.nnz ();
83 
84  for (octave_idx_type i = 0; i < nc + 1; i++)
85  cidx (i) = a.cidx (i);
86 
87  for (octave_idx_type i = 0; i < nz; i++)
88  {
89  data (i) = Complex (a.data (i));
90  ridx (i) = a.ridx (i);
91  }
92 }
93 
95  : MSparse<Complex> (a.rows (), a.cols (), a.length ())
96 {
97  octave_idx_type j = 0;
98  octave_idx_type l = a.length ();
99  for (octave_idx_type i = 0; i < l; i++)
100  {
101  cidx (i) = j;
102  if (a(i, i) != 0.0)
103  {
104  data (j) = a(i, i);
105  ridx (j) = i;
106  j++;
107  }
108  }
109  for (octave_idx_type i = l; i <= a.cols (); i++)
110  cidx (i) = j;
111 }
112 bool
114 {
115  octave_idx_type nr = rows ();
116  octave_idx_type nc = cols ();
117  octave_idx_type nz = nnz ();
118  octave_idx_type nr_a = a.rows ();
119  octave_idx_type nc_a = a.cols ();
120  octave_idx_type nz_a = a.nnz ();
121 
122  if (nr != nr_a || nc != nc_a || nz != nz_a)
123  return false;
124 
125  for (octave_idx_type i = 0; i < nc + 1; i++)
126  if (cidx (i) != a.cidx (i))
127  return false;
128 
129  for (octave_idx_type i = 0; i < nz; i++)
130  if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
131  return false;
132 
133  return true;
134 }
135 
136 bool
138 {
139  return !(*this == a);
140 }
141 
142 bool
144 {
145  octave_idx_type nr = rows ();
146  octave_idx_type nc = cols ();
147 
148  if (nr == nc && nr > 0)
149  {
150  for (octave_idx_type j = 0; j < nc; j++)
151  {
152  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
153  {
154  octave_idx_type ri = ridx (i);
155 
156  if (ri != j)
157  {
158  bool found = false;
159 
160  for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
161  {
162  if (ridx (k) == j)
163  {
164  if (data (i) == conj (data (k)))
165  found = true;
166  break;
167  }
168  }
169 
170  if (! found)
171  return false;
172  }
173  }
174  }
175 
176  return true;
177  }
178 
179  return false;
180 }
181 
182 static const Complex
185 
188 {
189  Array<octave_idx_type> dummy_idx;
190  return max (dummy_idx, dim);
191 }
192 
195 {
197  dim_vector dv = dims ();
198  octave_idx_type nr = dv(0);
199  octave_idx_type nc = dv(1);
200 
201  if (dim >= dv.ndims ())
202  {
203  idx_arg.resize (dim_vector (nr, nc), 0);
204  return *this;
205  }
206 
207  if (dim < 0)
208  dim = dv.first_non_singleton ();
209 
210  if (dim == 0)
211  {
212  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
213 
214  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
215  return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
216 
217  octave_idx_type nel = 0;
218  for (octave_idx_type j = 0; j < nc; j++)
219  {
220  Complex tmp_max;
221  double abs_max = octave::numeric_limits<double>::NaN ();
222  octave_idx_type idx_j = 0;
223  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
224  {
225  if (ridx (i) != idx_j)
226  break;
227  else
228  idx_j++;
229  }
230 
231  if (idx_j != nr)
232  {
233  tmp_max = 0.;
234  abs_max = 0.;
235  }
236 
237  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
238  {
239  Complex tmp = data (i);
240 
241  if (octave::math::isnan (tmp))
242  continue;
243 
244  double abs_tmp = std::abs (tmp);
245 
246  if (octave::math::isnan (abs_max) || abs_tmp > abs_max)
247  {
248  idx_j = ridx (i);
249  tmp_max = tmp;
250  abs_max = abs_tmp;
251  }
252  }
253 
254  idx_arg.elem (j) = octave::math::isnan (tmp_max) ? 0 : idx_j;
255  if (abs_max != 0.)
256  nel++;
257  }
258 
259  result = SparseComplexMatrix (1, nc, nel);
260 
261  octave_idx_type ii = 0;
262  result.xcidx (0) = 0;
263  for (octave_idx_type j = 0; j < nc; j++)
264  {
265  Complex tmp = elem (idx_arg(j), j);
266  if (tmp != 0.)
267  {
268  result.xdata (ii) = tmp;
269  result.xridx (ii++) = 0;
270  }
271  result.xcidx (j+1) = ii;
272  }
273  }
274  else
275  {
276  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
277 
278  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
279  return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
280 
282 
283  for (octave_idx_type i = 0; i < nr; i++)
284  found[i] = 0;
285 
286  for (octave_idx_type j = 0; j < nc; j++)
287  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
288  if (found[ridx (i)] == -j)
289  found[ridx (i)] = -j - 1;
290 
291  for (octave_idx_type i = 0; i < nr; i++)
292  if (found[i] > -nc && found[i] < 0)
293  idx_arg.elem (i) = -found[i];
294 
295  for (octave_idx_type j = 0; j < nc; j++)
296  {
297  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
298  {
299  octave_idx_type ir = ridx (i);
300  octave_idx_type ix = idx_arg.elem (ir);
301  Complex tmp = data (i);
302 
303  if (octave::math::isnan (tmp))
304  continue;
305  else if (ix == -1 || std::abs (tmp) > std::abs (elem (ir, ix)))
306  idx_arg.elem (ir) = j;
307  }
308  }
309 
310  octave_idx_type nel = 0;
311  for (octave_idx_type j = 0; j < nr; j++)
312  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
313  nel++;
314 
315  result = SparseComplexMatrix (nr, 1, nel);
316 
317  octave_idx_type ii = 0;
318  result.xcidx (0) = 0;
319  result.xcidx (1) = nel;
320  for (octave_idx_type j = 0; j < nr; j++)
321  {
322  if (idx_arg(j) == -1)
323  {
324  idx_arg(j) = 0;
325  result.xdata (ii) = Complex_NaN_result;
326  result.xridx (ii++) = j;
327  }
328  else
329  {
330  Complex tmp = elem (j, idx_arg(j));
331  if (tmp != 0.)
332  {
333  result.xdata (ii) = tmp;
334  result.xridx (ii++) = j;
335  }
336  }
337  }
338  }
339 
340  return result;
341 }
342 
345 {
346  Array<octave_idx_type> dummy_idx;
347  return min (dummy_idx, dim);
348 }
349 
352 {
354  dim_vector dv = dims ();
355  octave_idx_type nr = dv(0);
356  octave_idx_type nc = dv(1);
357 
358  if (dim >= dv.ndims ())
359  {
360  idx_arg.resize (dim_vector (nr, nc), 0);
361  return *this;
362  }
363 
364  if (dim < 0)
365  dim = dv.first_non_singleton ();
366 
367  if (dim == 0)
368  {
369  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
370 
371  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
372  return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
373 
374  octave_idx_type nel = 0;
375  for (octave_idx_type j = 0; j < nc; j++)
376  {
377  Complex tmp_min;
378  double abs_min = octave::numeric_limits<double>::NaN ();
379  octave_idx_type idx_j = 0;
380  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
381  {
382  if (ridx (i) != idx_j)
383  break;
384  else
385  idx_j++;
386  }
387 
388  if (idx_j != nr)
389  {
390  tmp_min = 0.;
391  abs_min = 0.;
392  }
393 
394  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
395  {
396  Complex tmp = data (i);
397 
398  if (octave::math::isnan (tmp))
399  continue;
400 
401  double abs_tmp = std::abs (tmp);
402 
403  if (octave::math::isnan (abs_min) || abs_tmp < abs_min)
404  {
405  idx_j = ridx (i);
406  tmp_min = tmp;
407  abs_min = abs_tmp;
408  }
409  }
410 
411  idx_arg.elem (j) = octave::math::isnan (tmp_min) ? 0 : idx_j;
412  if (abs_min != 0.)
413  nel++;
414  }
415 
416  result = SparseComplexMatrix (1, nc, nel);
417 
418  octave_idx_type ii = 0;
419  result.xcidx (0) = 0;
420  for (octave_idx_type j = 0; j < nc; j++)
421  {
422  Complex tmp = elem (idx_arg(j), j);
423  if (tmp != 0.)
424  {
425  result.xdata (ii) = tmp;
426  result.xridx (ii++) = 0;
427  }
428  result.xcidx (j+1) = ii;
429  }
430  }
431  else
432  {
433  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
434 
435  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
436  return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
437 
439 
440  for (octave_idx_type i = 0; i < nr; i++)
441  found[i] = 0;
442 
443  for (octave_idx_type j = 0; j < nc; j++)
444  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
445  if (found[ridx (i)] == -j)
446  found[ridx (i)] = -j - 1;
447 
448  for (octave_idx_type i = 0; i < nr; i++)
449  if (found[i] > -nc && found[i] < 0)
450  idx_arg.elem (i) = -found[i];
451 
452  for (octave_idx_type j = 0; j < nc; j++)
453  {
454  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
455  {
456  octave_idx_type ir = ridx (i);
457  octave_idx_type ix = idx_arg.elem (ir);
458  Complex tmp = data (i);
459 
460  if (octave::math::isnan (tmp))
461  continue;
462  else if (ix == -1 || std::abs (tmp) < std::abs (elem (ir, ix)))
463  idx_arg.elem (ir) = j;
464  }
465  }
466 
467  octave_idx_type nel = 0;
468  for (octave_idx_type j = 0; j < nr; j++)
469  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
470  nel++;
471 
472  result = SparseComplexMatrix (nr, 1, nel);
473 
474  octave_idx_type ii = 0;
475  result.xcidx (0) = 0;
476  result.xcidx (1) = nel;
477  for (octave_idx_type j = 0; j < nr; j++)
478  {
479  if (idx_arg(j) == -1)
480  {
481  idx_arg(j) = 0;
482  result.xdata (ii) = Complex_NaN_result;
483  result.xridx (ii++) = j;
484  }
485  else
486  {
487  Complex tmp = elem (j, idx_arg(j));
488  if (tmp != 0.)
489  {
490  result.xdata (ii) = tmp;
491  result.xridx (ii++) = j;
492  }
493  }
494  }
495  }
496 
497  return result;
498 }
499 
500 /*
501 
502 %!assert (max (max (speye (65536) * 1i)), sparse (1i))
503 %!assert (min (min (speye (65536) * 1i)), sparse (0))
504 %!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
505 %!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
506 %!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
507 %!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
508 %!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
509 %!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
510 %!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
511 %!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
512 
513 */
514 
517 {
518  octave_idx_type nc = columns ();
519  ComplexRowVector retval (nc, 0);
520 
521  for (octave_idx_type j = 0; j < nc; j++)
522  for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
523  {
524  if (ridx (k) == i)
525  {
526  retval(j) = data (k);
527  break;
528  }
529  }
530 
531  return retval;
532 }
533 
536 {
537  octave_idx_type nr = rows ();
538  ComplexColumnVector retval (nr, 0);
539 
540  for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
541  retval(ridx (k)) = data (k);
542 
543  return retval;
544 }
545 
546 // destructive insert/delete/reorder operations
547 
551 {
553  return insert (tmp /*a*/, r, c);
554 }
555 
559 {
560  MSparse<Complex>::insert (a, r, c);
561  return *this;
562 }
563 
566  const Array<octave_idx_type>& indx)
567 {
569  return insert (tmp /*a*/, indx);
570 }
571 
574  const Array<octave_idx_type>& indx)
575 {
576  MSparse<Complex>::insert (a, indx);
577  return *this;
578 }
579 
583 {
584  // Don't use numel to avoid all possiblity of an overflow
585  if (rb.rows () > 0 && rb.cols () > 0)
586  insert (rb, ra_idx(0), ra_idx(1));
587  return *this;
588 }
589 
593 {
595  if (rb.rows () > 0 && rb.cols () > 0)
596  insert (tmp, ra_idx(0), ra_idx(1));
597  return *this;
598 }
599 
602 {
604 }
605 
608 {
609  octave_idx_type nr = rows ();
610  octave_idx_type nc = cols ();
611  octave_idx_type nz = nnz ();
612  SparseComplexMatrix retval (nc, nr, nz);
613 
614  for (octave_idx_type i = 0; i < nz; i++)
615  retval.xcidx (ridx (i) + 1)++;
616  // retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
617  nz = 0;
618  for (octave_idx_type i = 1; i <= nr; i++)
619  {
620  const octave_idx_type tmp = retval.xcidx (i);
621  retval.xcidx (i) = nz;
622  nz += tmp;
623  }
624  // retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
625 
626  for (octave_idx_type j = 0; j < nc; j++)
627  for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
628  {
629  octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
630  retval.xridx (q) = j;
631  retval.xdata (q) = conj (data (k));
632  }
633  assert (nnz () == retval.xcidx (nr));
634  // retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
635  // and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
636 
637  return retval;
638 }
639 
642 {
643  octave_idx_type nr = a.rows ();
644  octave_idx_type nc = a.cols ();
645  octave_idx_type nz = a.nnz ();
646  SparseComplexMatrix retval (nc, nr, nz);
647 
648  for (octave_idx_type i = 0; i < nc + 1; i++)
649  retval.cidx (i) = a.cidx (i);
650 
651  for (octave_idx_type i = 0; i < nz; i++)
652  {
653  retval.data (i) = conj (a.data (i));
654  retval.ridx (i) = a.ridx (i);
655  }
656 
657  return retval;
658 }
659 
662 {
663  octave_idx_type info;
664  double rcond;
665  MatrixType mattype (*this);
666  return inverse (mattype, info, rcond, 0, 0);
667 }
668 
671 {
672  octave_idx_type info;
673  double rcond;
674  return inverse (mattype, info, rcond, 0, 0);
675 }
676 
679 {
680  double rcond;
681  return inverse (mattype, info, rcond, 0, 0);
682 }
683 
686  double& rcond, const bool,
687  const bool calccond) const
688 {
690 
691  octave_idx_type nr = rows ();
692  octave_idx_type nc = cols ();
693  info = 0;
694 
695  if (nr == 0 || nc == 0 || nr != nc)
696  (*current_liboctave_error_handler) ("inverse requires square matrix");
697 
698  // Print spparms("spumoni") info if requested
699  int typ = mattyp.type ();
700  mattyp.info ();
701 
703  (*current_liboctave_error_handler) ("incorrect matrix type");
704 
706  retval = transpose ();
707  else
708  retval = *this;
709 
710  // Force make_unique to be called
711  Complex *v = retval.data ();
712 
713  if (calccond)
714  {
715  double dmax = 0.;
716  double dmin = octave::numeric_limits<double>::Inf ();
717  for (octave_idx_type i = 0; i < nr; i++)
718  {
719  double tmp = std::abs (v[i]);
720  if (tmp > dmax)
721  dmax = tmp;
722  if (tmp < dmin)
723  dmin = tmp;
724  }
725  rcond = dmin / dmax;
726  }
727 
728  for (octave_idx_type i = 0; i < nr; i++)
729  v[i] = 1.0 / v[i];
730 
731  return retval;
732 }
733 
736  double& rcond, const bool,
737  const bool calccond) const
738 {
740 
741  octave_idx_type nr = rows ();
742  octave_idx_type nc = cols ();
743  info = 0;
744 
745  if (nr == 0 || nc == 0 || nr != nc)
746  (*current_liboctave_error_handler) ("inverse requires square matrix");
747 
748  // Print spparms("spumoni") info if requested
749  int typ = mattyp.type ();
750  mattyp.info ();
751 
752  if (typ != MatrixType::Upper && typ != MatrixType::Permuted_Upper
753  && typ != MatrixType::Lower && typ != MatrixType::Permuted_Lower)
754  (*current_liboctave_error_handler) ("incorrect matrix type");
755 
756  double anorm = 0.;
757  double ainvnorm = 0.;
758 
759  if (calccond)
760  {
761  // Calculate the 1-norm of matrix for rcond calculation
762  for (octave_idx_type j = 0; j < nr; j++)
763  {
764  double atmp = 0.;
765  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
766  atmp += std::abs (data (i));
767  if (atmp > anorm)
768  anorm = atmp;
769  }
770  }
771 
772  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
773  {
774  octave_idx_type nz = nnz ();
775  octave_idx_type cx = 0;
776  octave_idx_type nz2 = nz;
777  retval = SparseComplexMatrix (nr, nc, nz2);
778 
779  for (octave_idx_type i = 0; i < nr; i++)
780  {
781  octave_quit ();
782  // place the 1 in the identity position
783  octave_idx_type cx_colstart = cx;
784 
785  if (cx == nz2)
786  {
787  nz2 *= 2;
788  retval.change_capacity (nz2);
789  }
790 
791  retval.xcidx (i) = cx;
792  retval.xridx (cx) = i;
793  retval.xdata (cx) = 1.0;
794  cx++;
795 
796  // iterate accross columns of input matrix
797  for (octave_idx_type j = i+1; j < nr; j++)
798  {
799  Complex v = 0.;
800  // iterate to calculate sum
801  octave_idx_type colXp = retval.xcidx (i);
802  octave_idx_type colUp = cidx (j);
803  octave_idx_type rpX, rpU;
804 
805  if (cidx (j) == cidx (j+1))
806  (*current_liboctave_error_handler) ("division by zero");
807 
808  do
809  {
810  octave_quit ();
811  rpX = retval.xridx (colXp);
812  rpU = ridx (colUp);
813 
814  if (rpX < rpU)
815  colXp++;
816  else if (rpX > rpU)
817  colUp++;
818  else
819  {
820  v -= retval.xdata (colXp) * data (colUp);
821  colXp++;
822  colUp++;
823  }
824  }
825  while (rpX < j && rpU < j && colXp < cx && colUp < nz);
826 
827  // get A(m,m)
828  if (typ == MatrixType::Upper)
829  colUp = cidx (j+1) - 1;
830  else
831  colUp = cidx (j);
832  Complex pivot = data (colUp);
833  if (pivot == 0. || ridx (colUp) != j)
834  (*current_liboctave_error_handler) ("division by zero");
835 
836  if (v != 0.)
837  {
838  if (cx == nz2)
839  {
840  nz2 *= 2;
841  retval.change_capacity (nz2);
842  }
843 
844  retval.xridx (cx) = j;
845  retval.xdata (cx) = v / pivot;
846  cx++;
847  }
848  }
849 
850  // get A(m,m)
851  octave_idx_type colUp;
852  if (typ == MatrixType::Upper)
853  colUp = cidx (i+1) - 1;
854  else
855  colUp = cidx (i);
856  Complex pivot = data (colUp);
857  if (pivot == 0. || ridx (colUp) != i)
858  (*current_liboctave_error_handler) ("division by zero");
859 
860  if (pivot != 1.0)
861  for (octave_idx_type j = cx_colstart; j < cx; j++)
862  retval.xdata (j) /= pivot;
863  }
864  retval.xcidx (nr) = cx;
865  retval.maybe_compress ();
866  }
867  else
868  {
869  octave_idx_type nz = nnz ();
870  octave_idx_type cx = 0;
871  octave_idx_type nz2 = nz;
872  retval = SparseComplexMatrix (nr, nc, nz2);
873 
874  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
876 
877  octave_idx_type *perm = mattyp.triangular_perm ();
878  if (typ == MatrixType::Permuted_Upper)
879  {
880  for (octave_idx_type i = 0; i < nr; i++)
881  rperm[perm[i]] = i;
882  }
883  else
884  {
885  for (octave_idx_type i = 0; i < nr; i++)
886  rperm[i] = perm[i];
887  for (octave_idx_type i = 0; i < nr; i++)
888  perm[rperm[i]] = i;
889  }
890 
891  for (octave_idx_type i = 0; i < nr; i++)
892  {
893  octave_quit ();
894  octave_idx_type iidx = rperm[i];
895 
896  for (octave_idx_type j = 0; j < nr; j++)
897  work[j] = 0.;
898 
899  // place the 1 in the identity position
900  work[iidx] = 1.0;
901 
902  // iterate accross columns of input matrix
903  for (octave_idx_type j = iidx+1; j < nr; j++)
904  {
905  Complex v = 0.;
906  octave_idx_type jidx = perm[j];
907  // iterate to calculate sum
908  for (octave_idx_type k = cidx (jidx);
909  k < cidx (jidx+1); k++)
910  {
911  octave_quit ();
912  v -= work[ridx (k)] * data (k);
913  }
914 
915  // get A(m,m)
916  Complex pivot;
917  if (typ == MatrixType::Permuted_Upper)
918  pivot = data (cidx (jidx+1) - 1);
919  else
920  pivot = data (cidx (jidx));
921  if (pivot == 0.)
922  (*current_liboctave_error_handler) ("division by zero");
923 
924  work[j] = v / pivot;
925  }
926 
927  // get A(m,m)
928  octave_idx_type colUp;
929  if (typ == MatrixType::Permuted_Upper)
930  colUp = cidx (perm[iidx]+1) - 1;
931  else
932  colUp = cidx (perm[iidx]);
933 
934  Complex pivot = data (colUp);
935  if (pivot == 0.)
936  (*current_liboctave_error_handler) ("division by zero");
937 
938  octave_idx_type new_cx = cx;
939  for (octave_idx_type j = iidx; j < nr; j++)
940  if (work[j] != 0.0)
941  {
942  new_cx++;
943  if (pivot != 1.0)
944  work[j] /= pivot;
945  }
946 
947  if (cx < new_cx)
948  {
949  nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
950  retval.change_capacity (nz2);
951  }
952 
953  retval.xcidx (i) = cx;
954  for (octave_idx_type j = iidx; j < nr; j++)
955  if (work[j] != 0.)
956  {
957  retval.xridx (cx) = j;
958  retval.xdata (cx++) = work[j];
959  }
960  }
961 
962  retval.xcidx (nr) = cx;
963  retval.maybe_compress ();
964  }
965 
966  if (calccond)
967  {
968  // Calculate the 1-norm of inverse matrix for rcond calculation
969  for (octave_idx_type j = 0; j < nr; j++)
970  {
971  double atmp = 0.;
972  for (octave_idx_type i = retval.cidx (j);
973  i < retval.cidx (j+1); i++)
974  atmp += std::abs (retval.data (i));
975  if (atmp > ainvnorm)
976  ainvnorm = atmp;
977  }
978 
979  rcond = 1. / ainvnorm / anorm;
980  }
981 
982  return retval;
983 }
984 
987  double& rcond, bool, bool calc_cond) const
988 {
989  int typ = mattype.type (false);
991 
992  if (typ == MatrixType::Unknown)
993  typ = mattype.type (*this);
994 
996  ret = dinverse (mattype, info, rcond, true, calc_cond);
997  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
998  ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
999  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1000  {
1001  MatrixType newtype = mattype.transpose ();
1002  ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
1003  }
1004  else
1005  {
1006  if (mattype.is_hermitian ())
1007  {
1008  MatrixType tmp_typ (MatrixType::Upper);
1009  octave::math::sparse_chol<SparseComplexMatrix> fact (*this, info, false);
1010  rcond = fact.rcond ();
1011  if (info == 0)
1012  {
1013  double rcond2;
1014  SparseMatrix Q = fact.Q ();
1015  SparseComplexMatrix InvL = fact.L ().transpose ().
1016  tinverse (tmp_typ, info, rcond2,
1017  true, false);
1018  ret = Q * InvL.hermitian () * InvL * Q.transpose ();
1019  }
1020  else
1021  {
1022  // Matrix is either singular or not positive definite
1023  mattype.mark_as_unsymmetric ();
1024  }
1025  }
1026 
1027  if (! mattype.is_hermitian ())
1028  {
1029  octave_idx_type n = rows ();
1030  ColumnVector Qinit(n);
1031  for (octave_idx_type i = 0; i < n; i++)
1032  Qinit(i) = i;
1033 
1034  MatrixType tmp_typ (MatrixType::Upper);
1036  Qinit, Matrix (),
1037  false, false);
1038  rcond = fact.rcond ();
1039  double rcond2;
1040  SparseComplexMatrix InvL = fact.L ().transpose ().
1041  tinverse (tmp_typ, info, rcond2,
1042  true, false);
1043  SparseComplexMatrix InvU = fact.U ().
1044  tinverse (tmp_typ, info, rcond2,
1045  true, false).transpose ();
1046  ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1047  }
1048  }
1049 
1050  return ret;
1051 }
1052 
1053 ComplexDET
1055 {
1056  octave_idx_type info;
1057  double rcond;
1058  return determinant (info, rcond, 0);
1059 }
1060 
1061 ComplexDET
1063 {
1064  double rcond;
1065  return determinant (info, rcond, 0);
1066 }
1067 
1068 ComplexDET
1070  bool) const
1071 {
1073 
1074 #if defined (HAVE_UMFPACK)
1075 
1076  octave_idx_type nr = rows ();
1077  octave_idx_type nc = cols ();
1078 
1079  if (nr == 0 || nc == 0 || nr != nc)
1080  {
1081  retval = ComplexDET (1.0);
1082  }
1083  else
1084  {
1085  err = 0;
1086 
1087  // Setup the control parameters
1088  Matrix Control (UMFPACK_CONTROL, 1);
1089  double *control = Control.fortran_vec ();
1090  UMFPACK_ZNAME (defaults) (control);
1091 
1092  double tmp = octave_sparse_params::get_key ("spumoni");
1093  if (! octave::math::isnan (tmp))
1094  Control (UMFPACK_PRL) = tmp;
1095 
1096  tmp = octave_sparse_params::get_key ("piv_tol");
1097  if (! octave::math::isnan (tmp))
1098  {
1099  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1100  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1101  }
1102 
1103  // Set whether we are allowed to modify Q or not
1104  tmp = octave_sparse_params::get_key ("autoamd");
1105  if (! octave::math::isnan (tmp))
1106  Control (UMFPACK_FIXQ) = tmp;
1107 
1108  // Turn-off UMFPACK scaling for LU
1109  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1110 
1111  UMFPACK_ZNAME (report_control) (control);
1112 
1113  const octave_idx_type *Ap = cidx ();
1114  const octave_idx_type *Ai = ridx ();
1115  const Complex *Ax = data ();
1116 
1117  UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
1118  reinterpret_cast<const double *> (Ax),
1119  0, 1, control);
1120 
1121  void *Symbolic;
1122  Matrix Info (1, UMFPACK_INFO);
1123  double *info = Info.fortran_vec ();
1124  int status = UMFPACK_ZNAME (qsymbolic)
1125  (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0,
1126  0, &Symbolic, control, info);
1127 
1128  if (status < 0)
1129  {
1130  UMFPACK_ZNAME (report_status) (control, status);
1131  UMFPACK_ZNAME (report_info) (control, info);
1132 
1133  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1134 
1135  (*current_liboctave_error_handler)
1136  ("SparseComplexMatrix::determinant symbolic factorization failed");
1137  }
1138  else
1139  {
1140  UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
1141 
1142  void *Numeric;
1143  status
1144  = UMFPACK_ZNAME (numeric) (Ap, Ai,
1145  reinterpret_cast<const double *> (Ax),
1146  0, Symbolic, &Numeric, control, info);
1147  UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1148 
1149  rcond = Info (UMFPACK_RCOND);
1150 
1151  if (status < 0)
1152  {
1153  UMFPACK_ZNAME (report_status) (control, status);
1154  UMFPACK_ZNAME (report_info) (control, info);
1155 
1156  UMFPACK_ZNAME (free_numeric) (&Numeric);
1157 
1158  (*current_liboctave_error_handler)
1159  ("SparseComplexMatrix::determinant numeric factorization failed");
1160  }
1161  else
1162  {
1163  UMFPACK_ZNAME (report_numeric) (Numeric, control);
1164 
1165  double c10[2], e10;
1166 
1167  status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10,
1168  Numeric, info);
1169 
1170  if (status < 0)
1171  {
1172  UMFPACK_ZNAME (report_status) (control, status);
1173  UMFPACK_ZNAME (report_info) (control, info);
1174 
1175  (*current_liboctave_error_handler)
1176  ("SparseComplexMatrix::determinant error calculating determinant");
1177  }
1178  else
1179  retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10);
1180 
1181  UMFPACK_ZNAME (free_numeric) (&Numeric);
1182  }
1183  }
1184  }
1185 
1186 #else
1187 
1188  octave_unused_parameter (err);
1189  octave_unused_parameter (rcond);
1190 
1191  (*current_liboctave_error_handler)
1192  ("support for UMFPACK was unavailable or disabled when liboctave was built");
1193 
1194 #endif
1195 
1196  return retval;
1197 }
1198 
1201  octave_idx_type& err, double& rcond,
1202  solve_singularity_handler, bool calc_cond) const
1203 {
1205 
1206  octave_idx_type nr = rows ();
1207  octave_idx_type nc = cols ();
1208  octave_idx_type nm = (nc < nr ? nc : nr);
1209  err = 0;
1210 
1211  if (nr != b.rows ())
1213  ("matrix dimension mismatch solution of linear equations");
1214 
1215  if (nr == 0 || nc == 0 || b.cols () == 0)
1216  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1217  else
1218  {
1219  // Print spparms("spumoni") info if requested
1220  int typ = mattype.type ();
1221  mattype.info ();
1222 
1224  (*current_liboctave_error_handler) ("incorrect matrix type");
1225 
1226  retval.resize (nc, b.cols (), Complex (0.,0.));
1227  if (typ == MatrixType::Diagonal)
1228  for (octave_idx_type j = 0; j < b.cols (); j++)
1229  for (octave_idx_type i = 0; i < nm; i++)
1230  retval(i,j) = b(i,j) / data (i);
1231  else
1232  for (octave_idx_type j = 0; j < b.cols (); j++)
1233  for (octave_idx_type k = 0; k < nc; k++)
1234  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1235  retval(k,j) = b(ridx (i),j) / data (i);
1236 
1237  if (calc_cond)
1238  {
1239  double dmax = 0.;
1240  double dmin = octave::numeric_limits<double>::Inf ();
1241  for (octave_idx_type i = 0; i < nm; i++)
1242  {
1243  double tmp = std::abs (data (i));
1244  if (tmp > dmax)
1245  dmax = tmp;
1246  if (tmp < dmin)
1247  dmin = tmp;
1248  }
1249  rcond = dmin / dmax;
1250  }
1251  else
1252  rcond = 1.0;
1253  }
1254 
1255  return retval;
1256 }
1257 
1260  octave_idx_type& err, double& rcond,
1261  solve_singularity_handler,
1262  bool calc_cond) const
1263 {
1265 
1266  octave_idx_type nr = rows ();
1267  octave_idx_type nc = cols ();
1268  octave_idx_type nm = (nc < nr ? nc : nr);
1269  err = 0;
1270 
1271  if (nr != b.rows ())
1273  ("matrix dimension mismatch solution of linear equations");
1274 
1275  if (nr == 0 || nc == 0 || b.cols () == 0)
1276  retval = SparseComplexMatrix (nc, b.cols ());
1277  else
1278  {
1279  // Print spparms("spumoni") info if requested
1280  int typ = mattype.type ();
1281  mattype.info ();
1282 
1284  (*current_liboctave_error_handler) ("incorrect matrix type");
1285 
1286  octave_idx_type b_nc = b.cols ();
1287  octave_idx_type b_nz = b.nnz ();
1288  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1289 
1290  retval.xcidx (0) = 0;
1291  octave_idx_type ii = 0;
1292  if (typ == MatrixType::Diagonal)
1293  for (octave_idx_type j = 0; j < b.cols (); j++)
1294  {
1295  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1296  {
1297  if (b.ridx (i) >= nm)
1298  break;
1299  retval.xridx (ii) = b.ridx (i);
1300  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1301  }
1302  retval.xcidx (j+1) = ii;
1303  }
1304  else
1305  for (octave_idx_type j = 0; j < b.cols (); j++)
1306  {
1307  for (octave_idx_type l = 0; l < nc; l++)
1308  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1309  {
1310  bool found = false;
1312  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1313  if (ridx (i) == b.ridx (k))
1314  {
1315  found = true;
1316  break;
1317  }
1318  if (found)
1319  {
1320  retval.xridx (ii) = l;
1321  retval.xdata (ii++) = b.data (k) / data (i);
1322  }
1323  }
1324  retval.xcidx (j+1) = ii;
1325  }
1326 
1327  if (calc_cond)
1328  {
1329  double dmax = 0.;
1330  double dmin = octave::numeric_limits<double>::Inf ();
1331  for (octave_idx_type i = 0; i < nm; i++)
1332  {
1333  double tmp = std::abs (data (i));
1334  if (tmp > dmax)
1335  dmax = tmp;
1336  if (tmp < dmin)
1337  dmin = tmp;
1338  }
1339  rcond = dmin / dmax;
1340  }
1341  else
1342  rcond = 1.0;
1343  }
1344 
1345  return retval;
1346 }
1347 
1350  octave_idx_type& err, double& rcond,
1351  solve_singularity_handler,
1352  bool calc_cond) const
1353 {
1355 
1356  octave_idx_type nr = rows ();
1357  octave_idx_type nc = cols ();
1358  octave_idx_type nm = (nc < nr ? nc : nr);
1359  err = 0;
1360 
1361  if (nr != b.rows ())
1363  ("matrix dimension mismatch solution of linear equations");
1364 
1365  if (nr == 0 || nc == 0 || b.cols () == 0)
1366  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1367  else
1368  {
1369  // Print spparms("spumoni") info if requested
1370  int typ = mattype.type ();
1371  mattype.info ();
1372 
1374  (*current_liboctave_error_handler) ("incorrect matrix type");
1375 
1376  retval.resize (nc, b.cols (), Complex (0.,0.));
1377  if (typ == MatrixType::Diagonal)
1378  for (octave_idx_type j = 0; j < b.cols (); j++)
1379  for (octave_idx_type i = 0; i < nm; i++)
1380  retval(i,j) = b(i,j) / data (i);
1381  else
1382  for (octave_idx_type j = 0; j < b.cols (); j++)
1383  for (octave_idx_type k = 0; k < nc; k++)
1384  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1385  retval(k,j) = b(ridx (i),j) / data (i);
1386 
1387  if (calc_cond)
1388  {
1389  double dmax = 0.;
1390  double dmin = octave::numeric_limits<double>::Inf ();
1391  for (octave_idx_type i = 0; i < nr; i++)
1392  {
1393  double tmp = std::abs (data (i));
1394  if (tmp > dmax)
1395  dmax = tmp;
1396  if (tmp < dmin)
1397  dmin = tmp;
1398  }
1399  rcond = dmin / dmax;
1400  }
1401  else
1402  rcond = 1.0;
1403  }
1404 
1405  return retval;
1406 }
1407 
1410  octave_idx_type& err, double& rcond,
1411  solve_singularity_handler,
1412  bool calc_cond) const
1413 {
1415 
1416  octave_idx_type nr = rows ();
1417  octave_idx_type nc = cols ();
1418  octave_idx_type nm = (nc < nr ? nc : nr);
1419  err = 0;
1420 
1421  if (nr != b.rows ())
1423  ("matrix dimension mismatch solution of linear equations");
1424 
1425  if (nr == 0 || nc == 0 || b.cols () == 0)
1426  retval = SparseComplexMatrix (nc, b.cols ());
1427  else
1428  {
1429  // Print spparms("spumoni") info if requested
1430  int typ = mattype.type ();
1431  mattype.info ();
1432 
1434  (*current_liboctave_error_handler) ("incorrect matrix type");
1435 
1436  octave_idx_type b_nc = b.cols ();
1437  octave_idx_type b_nz = b.nnz ();
1438  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1439 
1440  retval.xcidx (0) = 0;
1441  octave_idx_type ii = 0;
1442  if (typ == MatrixType::Diagonal)
1443  for (octave_idx_type j = 0; j < b.cols (); j++)
1444  {
1445  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1446  {
1447  if (b.ridx (i) >= nm)
1448  break;
1449  retval.xridx (ii) = b.ridx (i);
1450  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1451  }
1452  retval.xcidx (j+1) = ii;
1453  }
1454  else
1455  for (octave_idx_type j = 0; j < b.cols (); j++)
1456  {
1457  for (octave_idx_type l = 0; l < nc; l++)
1458  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1459  {
1460  bool found = false;
1462  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1463  if (ridx (i) == b.ridx (k))
1464  {
1465  found = true;
1466  break;
1467  }
1468  if (found)
1469  {
1470  retval.xridx (ii) = l;
1471  retval.xdata (ii++) = b.data (k) / data (i);
1472  }
1473  }
1474  retval.xcidx (j+1) = ii;
1475  }
1476 
1477  if (calc_cond)
1478  {
1479  double dmax = 0.;
1480  double dmin = octave::numeric_limits<double>::Inf ();
1481  for (octave_idx_type i = 0; i < nm; i++)
1482  {
1483  double tmp = std::abs (data (i));
1484  if (tmp > dmax)
1485  dmax = tmp;
1486  if (tmp < dmin)
1487  dmin = tmp;
1488  }
1489  rcond = dmin / dmax;
1490  }
1491  else
1492  rcond = 1.0;
1493  }
1494 
1495  return retval;
1496 }
1497 
1500  octave_idx_type& err, double& rcond,
1501  solve_singularity_handler sing_handler,
1502  bool calc_cond) const
1503 {
1505 
1506  octave_idx_type nr = rows ();
1507  octave_idx_type nc = cols ();
1508  octave_idx_type nm = (nc > nr ? nc : nr);
1509  err = 0;
1510 
1511  if (nr != b.rows ())
1513  ("matrix dimension mismatch solution of linear equations");
1514 
1515  if (nr == 0 || nc == 0 || b.cols () == 0)
1516  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1517  else
1518  {
1519  // Print spparms("spumoni") info if requested
1520  int typ = mattype.type ();
1521  mattype.info ();
1522 
1523  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1524  (*current_liboctave_error_handler) ("incorrect matrix type");
1525 
1526  double anorm = 0.;
1527  double ainvnorm = 0.;
1528  octave_idx_type b_nc = b.cols ();
1529  rcond = 1.;
1530 
1531  if (calc_cond)
1532  {
1533  // Calculate the 1-norm of matrix for rcond calculation
1534  for (octave_idx_type j = 0; j < nc; j++)
1535  {
1536  double atmp = 0.;
1537  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1538  atmp += std::abs (data (i));
1539  if (atmp > anorm)
1540  anorm = atmp;
1541  }
1542  }
1543 
1544  if (typ == MatrixType::Permuted_Upper)
1545  {
1546  retval.resize (nc, b_nc);
1547  octave_idx_type *perm = mattype.triangular_perm ();
1548  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1549 
1550  for (octave_idx_type j = 0; j < b_nc; j++)
1551  {
1552  for (octave_idx_type i = 0; i < nr; i++)
1553  work[i] = b(i,j);
1554  for (octave_idx_type i = nr; i < nc; i++)
1555  work[i] = 0.;
1556 
1557  for (octave_idx_type k = nc-1; k >= 0; k--)
1558  {
1559  octave_idx_type kidx = perm[k];
1560 
1561  if (work[k] != 0.)
1562  {
1563  if (ridx (cidx (kidx+1)-1) != k
1564  || data (cidx (kidx+1)-1) == 0.)
1565  {
1566  err = -2;
1567  goto triangular_error;
1568  }
1569 
1570  Complex tmp = work[k] / data (cidx (kidx+1)-1);
1571  work[k] = tmp;
1572  for (octave_idx_type i = cidx (kidx);
1573  i < cidx (kidx+1)-1; i++)
1574  {
1575  octave_idx_type iidx = ridx (i);
1576  work[iidx] = work[iidx] - tmp * data (i);
1577  }
1578  }
1579  }
1580 
1581  for (octave_idx_type i = 0; i < nc; i++)
1582  retval(perm[i], j) = work[i];
1583  }
1584 
1585  if (calc_cond)
1586  {
1587  // Calculation of 1-norm of inv(*this)
1588  for (octave_idx_type i = 0; i < nm; i++)
1589  work[i] = 0.;
1590 
1591  for (octave_idx_type j = 0; j < nr; j++)
1592  {
1593  work[j] = 1.;
1594 
1595  for (octave_idx_type k = j; k >= 0; k--)
1596  {
1597  octave_idx_type iidx = perm[k];
1598 
1599  if (work[k] != 0.)
1600  {
1601  Complex tmp = work[k] / data (cidx (iidx+1)-1);
1602  work[k] = tmp;
1603  for (octave_idx_type i = cidx (iidx);
1604  i < cidx (iidx+1)-1; i++)
1605  {
1606  octave_idx_type idx2 = ridx (i);
1607  work[idx2] = work[idx2] - tmp * data (i);
1608  }
1609  }
1610  }
1611  double atmp = 0;
1612  for (octave_idx_type i = 0; i < j+1; i++)
1613  {
1614  atmp += std::abs (work[i]);
1615  work[i] = 0.;
1616  }
1617  if (atmp > ainvnorm)
1618  ainvnorm = atmp;
1619  }
1620  rcond = 1. / ainvnorm / anorm;
1621  }
1622  }
1623  else
1624  {
1625  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1626  retval.resize (nc, b_nc);
1627 
1628  for (octave_idx_type j = 0; j < b_nc; j++)
1629  {
1630  for (octave_idx_type i = 0; i < nr; i++)
1631  work[i] = b(i,j);
1632  for (octave_idx_type i = nr; i < nc; i++)
1633  work[i] = 0.;
1634 
1635  for (octave_idx_type k = nc-1; k >= 0; k--)
1636  {
1637  if (work[k] != 0.)
1638  {
1639  if (ridx (cidx (k+1)-1) != k
1640  || data (cidx (k+1)-1) == 0.)
1641  {
1642  err = -2;
1643  goto triangular_error;
1644  }
1645 
1646  Complex tmp = work[k] / data (cidx (k+1)-1);
1647  work[k] = tmp;
1648  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1649  {
1650  octave_idx_type iidx = ridx (i);
1651  work[iidx] = work[iidx] - tmp * data (i);
1652  }
1653  }
1654  }
1655 
1656  for (octave_idx_type i = 0; i < nc; i++)
1657  retval.xelem (i, j) = work[i];
1658  }
1659 
1660  if (calc_cond)
1661  {
1662  // Calculation of 1-norm of inv(*this)
1663  for (octave_idx_type i = 0; i < nm; i++)
1664  work[i] = 0.;
1665 
1666  for (octave_idx_type j = 0; j < nr; j++)
1667  {
1668  work[j] = 1.;
1669 
1670  for (octave_idx_type k = j; k >= 0; k--)
1671  {
1672  if (work[k] != 0.)
1673  {
1674  Complex tmp = work[k] / data (cidx (k+1)-1);
1675  work[k] = tmp;
1676  for (octave_idx_type i = cidx (k);
1677  i < cidx (k+1)-1; i++)
1678  {
1679  octave_idx_type iidx = ridx (i);
1680  work[iidx] = work[iidx] - tmp * data (i);
1681  }
1682  }
1683  }
1684  double atmp = 0;
1685  for (octave_idx_type i = 0; i < j+1; i++)
1686  {
1687  atmp += std::abs (work[i]);
1688  work[i] = 0.;
1689  }
1690  if (atmp > ainvnorm)
1691  ainvnorm = atmp;
1692  }
1693  rcond = 1. / ainvnorm / anorm;
1694  }
1695  }
1696 
1697  triangular_error:
1698  if (err != 0)
1699  {
1700  if (sing_handler)
1701  {
1702  sing_handler (rcond);
1703  mattype.mark_as_rectangular ();
1704  }
1705  else
1707  }
1708 
1709  volatile double rcond_plus_one = rcond + 1.0;
1710 
1711  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
1712  {
1713  err = -2;
1714 
1715  if (sing_handler)
1716  {
1717  sing_handler (rcond);
1718  mattype.mark_as_rectangular ();
1719  }
1720  else
1722  }
1723  }
1724 
1725  return retval;
1726 }
1727 
1730  octave_idx_type& err, double& rcond,
1731  solve_singularity_handler sing_handler,
1732  bool calc_cond) const
1733 {
1735 
1736  octave_idx_type nr = rows ();
1737  octave_idx_type nc = cols ();
1738  octave_idx_type nm = (nc > nr ? nc : nr);
1739  err = 0;
1740 
1741  if (nr != b.rows ())
1743  ("matrix dimension mismatch solution of linear equations");
1744 
1745  if (nr == 0 || nc == 0 || b.cols () == 0)
1746  retval = SparseComplexMatrix (nc, b.cols ());
1747  else
1748  {
1749  // Print spparms("spumoni") info if requested
1750  int typ = mattype.type ();
1751  mattype.info ();
1752 
1753  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1754  (*current_liboctave_error_handler) ("incorrect matrix type");
1755 
1756  double anorm = 0.;
1757  double ainvnorm = 0.;
1758  rcond = 1.;
1759 
1760  if (calc_cond)
1761  {
1762  // Calculate the 1-norm of matrix for rcond calculation
1763  for (octave_idx_type j = 0; j < nc; j++)
1764  {
1765  double atmp = 0.;
1766  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1767  atmp += std::abs (data (i));
1768  if (atmp > anorm)
1769  anorm = atmp;
1770  }
1771  }
1772 
1773  octave_idx_type b_nc = b.cols ();
1774  octave_idx_type b_nz = b.nnz ();
1775  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1776  retval.xcidx (0) = 0;
1777  octave_idx_type ii = 0;
1778  octave_idx_type x_nz = b_nz;
1779 
1780  if (typ == MatrixType::Permuted_Upper)
1781  {
1782  octave_idx_type *perm = mattype.triangular_perm ();
1783  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1784 
1785  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1786  for (octave_idx_type i = 0; i < nc; i++)
1787  rperm[perm[i]] = i;
1788 
1789  for (octave_idx_type j = 0; j < b_nc; j++)
1790  {
1791  for (octave_idx_type i = 0; i < nm; i++)
1792  work[i] = 0.;
1793  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1794  work[b.ridx (i)] = b.data (i);
1795 
1796  for (octave_idx_type k = nc-1; k >= 0; k--)
1797  {
1798  octave_idx_type kidx = perm[k];
1799 
1800  if (work[k] != 0.)
1801  {
1802  if (ridx (cidx (kidx+1)-1) != k
1803  || data (cidx (kidx+1)-1) == 0.)
1804  {
1805  err = -2;
1806  goto triangular_error;
1807  }
1808 
1809  Complex tmp = work[k] / data (cidx (kidx+1)-1);
1810  work[k] = tmp;
1811  for (octave_idx_type i = cidx (kidx);
1812  i < cidx (kidx+1)-1; i++)
1813  {
1814  octave_idx_type iidx = ridx (i);
1815  work[iidx] = work[iidx] - tmp * data (i);
1816  }
1817  }
1818  }
1819 
1820  // Count nonzeros in work vector and adjust space in
1821  // retval if needed
1822  octave_idx_type new_nnz = 0;
1823  for (octave_idx_type i = 0; i < nc; i++)
1824  if (work[i] != 0.)
1825  new_nnz++;
1826 
1827  if (ii + new_nnz > x_nz)
1828  {
1829  // Resize the sparse matrix
1830  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1831  retval.change_capacity (sz);
1832  x_nz = sz;
1833  }
1834 
1835  for (octave_idx_type i = 0; i < nc; i++)
1836  if (work[rperm[i]] != 0.)
1837  {
1838  retval.xridx (ii) = i;
1839  retval.xdata (ii++) = work[rperm[i]];
1840  }
1841  retval.xcidx (j+1) = ii;
1842  }
1843 
1844  retval.maybe_compress ();
1845 
1846  if (calc_cond)
1847  {
1848  // Calculation of 1-norm of inv(*this)
1849  for (octave_idx_type i = 0; i < nm; i++)
1850  work[i] = 0.;
1851 
1852  for (octave_idx_type j = 0; j < nr; j++)
1853  {
1854  work[j] = 1.;
1855 
1856  for (octave_idx_type k = j; k >= 0; k--)
1857  {
1858  octave_idx_type iidx = perm[k];
1859 
1860  if (work[k] != 0.)
1861  {
1862  Complex tmp = work[k] / data (cidx (iidx+1)-1);
1863  work[k] = tmp;
1864  for (octave_idx_type i = cidx (iidx);
1865  i < cidx (iidx+1)-1; i++)
1866  {
1867  octave_idx_type idx2 = ridx (i);
1868  work[idx2] = work[idx2] - tmp * data (i);
1869  }
1870  }
1871  }
1872  double atmp = 0;
1873  for (octave_idx_type i = 0; i < j+1; i++)
1874  {
1875  atmp += std::abs (work[i]);
1876  work[i] = 0.;
1877  }
1878  if (atmp > ainvnorm)
1879  ainvnorm = atmp;
1880  }
1881  rcond = 1. / ainvnorm / anorm;
1882  }
1883  }
1884  else
1885  {
1886  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1887 
1888  for (octave_idx_type j = 0; j < b_nc; j++)
1889  {
1890  for (octave_idx_type i = 0; i < nm; i++)
1891  work[i] = 0.;
1892  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1893  work[b.ridx (i)] = b.data (i);
1894 
1895  for (octave_idx_type k = nc-1; k >= 0; k--)
1896  {
1897  if (work[k] != 0.)
1898  {
1899  if (ridx (cidx (k+1)-1) != k
1900  || data (cidx (k+1)-1) == 0.)
1901  {
1902  err = -2;
1903  goto triangular_error;
1904  }
1905 
1906  Complex tmp = work[k] / data (cidx (k+1)-1);
1907  work[k] = tmp;
1908  for (octave_idx_type i = cidx (k); i < cidx (k+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[i] != 0.)
1933  {
1934  retval.xridx (ii) = i;
1935  retval.xdata (ii++) = work[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  if (work[k] != 0.)
1955  {
1956  Complex tmp = work[k] / data (cidx (k+1)-1);
1957  work[k] = tmp;
1958  for (octave_idx_type i = cidx (k);
1959  i < cidx (k+1)-1; i++)
1960  {
1961  octave_idx_type iidx = ridx (i);
1962  work[iidx] = work[iidx] - tmp * data (i);
1963  }
1964  }
1965  }
1966  double atmp = 0;
1967  for (octave_idx_type i = 0; i < j+1; i++)
1968  {
1969  atmp += std::abs (work[i]);
1970  work[i] = 0.;
1971  }
1972  if (atmp > ainvnorm)
1973  ainvnorm = atmp;
1974  }
1975  rcond = 1. / ainvnorm / anorm;
1976  }
1977  }
1978 
1979  triangular_error:
1980  if (err != 0)
1981  {
1982  if (sing_handler)
1983  {
1984  sing_handler (rcond);
1985  mattype.mark_as_rectangular ();
1986  }
1987  else
1989  }
1990 
1991  volatile double rcond_plus_one = rcond + 1.0;
1992 
1993  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
1994  {
1995  err = -2;
1996 
1997  if (sing_handler)
1998  {
1999  sing_handler (rcond);
2000  mattype.mark_as_rectangular ();
2001  }
2002  else
2004  }
2005  }
2006  return retval;
2007 }
2008 
2011  octave_idx_type& err, double& rcond,
2012  solve_singularity_handler sing_handler,
2013  bool calc_cond) const
2014 {
2016 
2017  octave_idx_type nr = rows ();
2018  octave_idx_type nc = cols ();
2019  octave_idx_type nm = (nc > nr ? nc : nr);
2020  err = 0;
2021 
2022  if (nr != b.rows ())
2024  ("matrix dimension mismatch solution of linear equations");
2025 
2026  if (nr == 0 || nc == 0 || b.cols () == 0)
2027  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2028  else
2029  {
2030  // Print spparms("spumoni") info if requested
2031  int typ = mattype.type ();
2032  mattype.info ();
2033 
2034  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2035  (*current_liboctave_error_handler) ("incorrect matrix type");
2036 
2037  double anorm = 0.;
2038  double ainvnorm = 0.;
2039  octave_idx_type b_nc = b.cols ();
2040  rcond = 1.;
2041 
2042  if (calc_cond)
2043  {
2044  // Calculate the 1-norm of matrix for rcond calculation
2045  for (octave_idx_type j = 0; j < nc; j++)
2046  {
2047  double atmp = 0.;
2048  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2049  atmp += std::abs (data (i));
2050  if (atmp > anorm)
2051  anorm = atmp;
2052  }
2053  }
2054 
2055  if (typ == MatrixType::Permuted_Upper)
2056  {
2057  retval.resize (nc, b_nc);
2058  octave_idx_type *perm = mattype.triangular_perm ();
2059  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2060 
2061  for (octave_idx_type j = 0; j < b_nc; j++)
2062  {
2063  for (octave_idx_type i = 0; i < nr; i++)
2064  work[i] = b(i,j);
2065  for (octave_idx_type i = nr; i < nc; i++)
2066  work[i] = 0.;
2067 
2068  for (octave_idx_type k = nc-1; k >= 0; k--)
2069  {
2070  octave_idx_type kidx = perm[k];
2071 
2072  if (work[k] != 0.)
2073  {
2074  if (ridx (cidx (kidx+1)-1) != k
2075  || data (cidx (kidx+1)-1) == 0.)
2076  {
2077  err = -2;
2078  goto triangular_error;
2079  }
2080 
2081  Complex tmp = work[k] / data (cidx (kidx+1)-1);
2082  work[k] = tmp;
2083  for (octave_idx_type i = cidx (kidx);
2084  i < cidx (kidx+1)-1; i++)
2085  {
2086  octave_idx_type iidx = ridx (i);
2087  work[iidx] = work[iidx] - tmp * data (i);
2088  }
2089  }
2090  }
2091 
2092  for (octave_idx_type i = 0; i < nc; i++)
2093  retval(perm[i], j) = work[i];
2094  }
2095 
2096  if (calc_cond)
2097  {
2098  // Calculation of 1-norm of inv(*this)
2099  for (octave_idx_type i = 0; i < nm; i++)
2100  work[i] = 0.;
2101 
2102  for (octave_idx_type j = 0; j < nr; j++)
2103  {
2104  work[j] = 1.;
2105 
2106  for (octave_idx_type k = j; k >= 0; k--)
2107  {
2108  octave_idx_type iidx = perm[k];
2109 
2110  if (work[k] != 0.)
2111  {
2112  Complex tmp = work[k] / data (cidx (iidx+1)-1);
2113  work[k] = tmp;
2114  for (octave_idx_type i = cidx (iidx);
2115  i < cidx (iidx+1)-1; i++)
2116  {
2117  octave_idx_type idx2 = ridx (i);
2118  work[idx2] = work[idx2] - tmp * data (i);
2119  }
2120  }
2121  }
2122  double atmp = 0;
2123  for (octave_idx_type i = 0; i < j+1; i++)
2124  {
2125  atmp += std::abs (work[i]);
2126  work[i] = 0.;
2127  }
2128  if (atmp > ainvnorm)
2129  ainvnorm = atmp;
2130  }
2131  rcond = 1. / ainvnorm / anorm;
2132  }
2133  }
2134  else
2135  {
2136  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2137  retval.resize (nc, b_nc);
2138 
2139  for (octave_idx_type j = 0; j < b_nc; j++)
2140  {
2141  for (octave_idx_type i = 0; i < nr; i++)
2142  work[i] = b(i,j);
2143  for (octave_idx_type i = nr; i < nc; i++)
2144  work[i] = 0.;
2145 
2146  for (octave_idx_type k = nc-1; k >= 0; k--)
2147  {
2148  if (work[k] != 0.)
2149  {
2150  if (ridx (cidx (k+1)-1) != k
2151  || data (cidx (k+1)-1) == 0.)
2152  {
2153  err = -2;
2154  goto triangular_error;
2155  }
2156 
2157  Complex tmp = work[k] / data (cidx (k+1)-1);
2158  work[k] = tmp;
2159  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2160  {
2161  octave_idx_type iidx = ridx (i);
2162  work[iidx] = work[iidx] - tmp * data (i);
2163  }
2164  }
2165  }
2166 
2167  for (octave_idx_type i = 0; i < nc; i++)
2168  retval.xelem (i, j) = work[i];
2169  }
2170 
2171  if (calc_cond)
2172  {
2173  // Calculation of 1-norm of inv(*this)
2174  for (octave_idx_type i = 0; i < nm; i++)
2175  work[i] = 0.;
2176 
2177  for (octave_idx_type j = 0; j < nr; j++)
2178  {
2179  work[j] = 1.;
2180 
2181  for (octave_idx_type k = j; k >= 0; k--)
2182  {
2183  if (work[k] != 0.)
2184  {
2185  Complex tmp = work[k] / data (cidx (k+1)-1);
2186  work[k] = tmp;
2187  for (octave_idx_type i = cidx (k);
2188  i < cidx (k+1)-1; i++)
2189  {
2190  octave_idx_type iidx = ridx (i);
2191  work[iidx] = work[iidx] - tmp * data (i);
2192  }
2193  }
2194  }
2195  double atmp = 0;
2196  for (octave_idx_type i = 0; i < j+1; i++)
2197  {
2198  atmp += std::abs (work[i]);
2199  work[i] = 0.;
2200  }
2201  if (atmp > ainvnorm)
2202  ainvnorm = atmp;
2203  }
2204  rcond = 1. / ainvnorm / anorm;
2205  }
2206  }
2207 
2208  triangular_error:
2209  if (err != 0)
2210  {
2211  if (sing_handler)
2212  {
2213  sing_handler (rcond);
2214  mattype.mark_as_rectangular ();
2215  }
2216  else
2218  }
2219 
2220  volatile double rcond_plus_one = rcond + 1.0;
2221 
2222  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2223  {
2224  err = -2;
2225 
2226  if (sing_handler)
2227  {
2228  sing_handler (rcond);
2229  mattype.mark_as_rectangular ();
2230  }
2231  else
2233  }
2234  }
2235 
2236  return retval;
2237 }
2238 
2241  octave_idx_type& err, double& rcond,
2242  solve_singularity_handler sing_handler,
2243  bool calc_cond) const
2244 {
2246 
2247  octave_idx_type nr = rows ();
2248  octave_idx_type nc = cols ();
2249  octave_idx_type nm = (nc > nr ? nc : nr);
2250  err = 0;
2251 
2252  if (nr != b.rows ())
2254  ("matrix dimension mismatch solution of linear equations");
2255 
2256  if (nr == 0 || nc == 0 || b.cols () == 0)
2257  retval = SparseComplexMatrix (nc, b.cols ());
2258  else
2259  {
2260  // Print spparms("spumoni") info if requested
2261  int typ = mattype.type ();
2262  mattype.info ();
2263 
2264  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2265  (*current_liboctave_error_handler) ("incorrect matrix type");
2266 
2267  double anorm = 0.;
2268  double ainvnorm = 0.;
2269  rcond = 1.;
2270 
2271  if (calc_cond)
2272  {
2273  // Calculate the 1-norm of matrix for rcond calculation
2274  for (octave_idx_type j = 0; j < nc; j++)
2275  {
2276  double atmp = 0.;
2277  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2278  atmp += std::abs (data (i));
2279  if (atmp > anorm)
2280  anorm = atmp;
2281  }
2282  }
2283 
2284  octave_idx_type b_nc = b.cols ();
2285  octave_idx_type b_nz = b.nnz ();
2286  retval = SparseComplexMatrix (nc, b_nc, b_nz);
2287  retval.xcidx (0) = 0;
2288  octave_idx_type ii = 0;
2289  octave_idx_type x_nz = b_nz;
2290 
2291  if (typ == MatrixType::Permuted_Upper)
2292  {
2293  octave_idx_type *perm = mattype.triangular_perm ();
2294  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2295 
2296  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
2297  for (octave_idx_type i = 0; i < nc; i++)
2298  rperm[perm[i]] = i;
2299 
2300  for (octave_idx_type j = 0; j < b_nc; j++)
2301  {
2302  for (octave_idx_type i = 0; i < nm; i++)
2303  work[i] = 0.;
2304  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2305  work[b.ridx (i)] = b.data (i);
2306 
2307  for (octave_idx_type k = nc-1; k >= 0; k--)
2308  {
2309  octave_idx_type kidx = perm[k];
2310 
2311  if (work[k] != 0.)
2312  {
2313  if (ridx (cidx (kidx+1)-1) != k
2314  || data (cidx (kidx+1)-1) == 0.)
2315  {
2316  err = -2;
2317  goto triangular_error;
2318  }
2319 
2320  Complex tmp = work[k] / data (cidx (kidx+1)-1);
2321  work[k] = tmp;
2322  for (octave_idx_type i = cidx (kidx);
2323  i < cidx (kidx+1)-1; i++)
2324  {
2325  octave_idx_type iidx = ridx (i);
2326  work[iidx] = work[iidx] - tmp * data (i);
2327  }
2328  }
2329  }
2330 
2331  // Count nonzeros in work vector and adjust space in
2332  // retval if needed
2333  octave_idx_type new_nnz = 0;
2334  for (octave_idx_type i = 0; i < nc; i++)
2335  if (work[i] != 0.)
2336  new_nnz++;
2337 
2338  if (ii + new_nnz > x_nz)
2339  {
2340  // Resize the sparse matrix
2341  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2342  retval.change_capacity (sz);
2343  x_nz = sz;
2344  }
2345 
2346  for (octave_idx_type i = 0; i < nc; i++)
2347  if (work[rperm[i]] != 0.)
2348  {
2349  retval.xridx (ii) = i;
2350  retval.xdata (ii++) = work[rperm[i]];
2351  }
2352  retval.xcidx (j+1) = ii;
2353  }
2354 
2355  retval.maybe_compress ();
2356 
2357  if (calc_cond)
2358  {
2359  // Calculation of 1-norm of inv(*this)
2360  for (octave_idx_type i = 0; i < nm; i++)
2361  work[i] = 0.;
2362 
2363  for (octave_idx_type j = 0; j < nr; j++)
2364  {
2365  work[j] = 1.;
2366 
2367  for (octave_idx_type k = j; k >= 0; k--)
2368  {
2369  octave_idx_type iidx = perm[k];
2370 
2371  if (work[k] != 0.)
2372  {
2373  Complex tmp = work[k] / data (cidx (iidx+1)-1);
2374  work[k] = tmp;
2375  for (octave_idx_type i = cidx (iidx);
2376  i < cidx (iidx+1)-1; i++)
2377  {
2378  octave_idx_type idx2 = ridx (i);
2379  work[idx2] = work[idx2] - tmp * data (i);
2380  }
2381  }
2382  }
2383  double atmp = 0;
2384  for (octave_idx_type i = 0; i < j+1; i++)
2385  {
2386  atmp += std::abs (work[i]);
2387  work[i] = 0.;
2388  }
2389  if (atmp > ainvnorm)
2390  ainvnorm = atmp;
2391  }
2392  rcond = 1. / ainvnorm / anorm;
2393  }
2394  }
2395  else
2396  {
2397  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2398 
2399  for (octave_idx_type j = 0; j < b_nc; j++)
2400  {
2401  for (octave_idx_type i = 0; i < nm; i++)
2402  work[i] = 0.;
2403  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2404  work[b.ridx (i)] = b.data (i);
2405 
2406  for (octave_idx_type k = nr-1; k >= 0; k--)
2407  {
2408  if (work[k] != 0.)
2409  {
2410  if (ridx (cidx (k+1)-1) != k
2411  || data (cidx (k+1)-1) == 0.)
2412  {
2413  err = -2;
2414  goto triangular_error;
2415  }
2416 
2417  Complex tmp = work[k] / data (cidx (k+1)-1);
2418  work[k] = tmp;
2419  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2420  {
2421  octave_idx_type iidx = ridx (i);
2422  work[iidx] = work[iidx] - tmp * data (i);
2423  }
2424  }
2425  }
2426 
2427  // Count nonzeros in work vector and adjust space in
2428  // retval if needed
2429  octave_idx_type new_nnz = 0;
2430  for (octave_idx_type i = 0; i < nc; i++)
2431  if (work[i] != 0.)
2432  new_nnz++;
2433 
2434  if (ii + new_nnz > x_nz)
2435  {
2436  // Resize the sparse matrix
2437  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2438  retval.change_capacity (sz);
2439  x_nz = sz;
2440  }
2441 
2442  for (octave_idx_type i = 0; i < nc; i++)
2443  if (work[i] != 0.)
2444  {
2445  retval.xridx (ii) = i;
2446  retval.xdata (ii++) = work[i];
2447  }
2448  retval.xcidx (j+1) = ii;
2449  }
2450 
2451  retval.maybe_compress ();
2452 
2453  if (calc_cond)
2454  {
2455  // Calculation of 1-norm of inv(*this)
2456  for (octave_idx_type i = 0; i < nm; i++)
2457  work[i] = 0.;
2458 
2459  for (octave_idx_type j = 0; j < nr; j++)
2460  {
2461  work[j] = 1.;
2462 
2463  for (octave_idx_type k = j; k >= 0; k--)
2464  {
2465  if (work[k] != 0.)
2466  {
2467  Complex tmp = work[k] / data (cidx (k+1)-1);
2468  work[k] = tmp;
2469  for (octave_idx_type i = cidx (k);
2470  i < cidx (k+1)-1; i++)
2471  {
2472  octave_idx_type iidx = ridx (i);
2473  work[iidx] = work[iidx] - tmp * data (i);
2474  }
2475  }
2476  }
2477  double atmp = 0;
2478  for (octave_idx_type i = 0; i < j+1; i++)
2479  {
2480  atmp += std::abs (work[i]);
2481  work[i] = 0.;
2482  }
2483  if (atmp > ainvnorm)
2484  ainvnorm = atmp;
2485  }
2486  rcond = 1. / ainvnorm / anorm;
2487  }
2488  }
2489 
2490  triangular_error:
2491  if (err != 0)
2492  {
2493  if (sing_handler)
2494  {
2495  sing_handler (rcond);
2496  mattype.mark_as_rectangular ();
2497  }
2498  else
2500  }
2501 
2502  volatile double rcond_plus_one = rcond + 1.0;
2503 
2504  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2505  {
2506  err = -2;
2507 
2508  if (sing_handler)
2509  {
2510  sing_handler (rcond);
2511  mattype.mark_as_rectangular ();
2512  }
2513  else
2515  }
2516  }
2517 
2518  return retval;
2519 }
2520 
2523  octave_idx_type& err, double& rcond,
2524  solve_singularity_handler sing_handler,
2525  bool calc_cond) const
2526 {
2528 
2529  octave_idx_type nr = rows ();
2530  octave_idx_type nc = cols ();
2531  octave_idx_type nm = (nc > nr ? nc : nr);
2532  err = 0;
2533 
2534  if (nr != b.rows ())
2536  ("matrix dimension mismatch solution of linear equations");
2537 
2538  if (nr == 0 || nc == 0 || b.cols () == 0)
2539  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2540  else
2541  {
2542  // Print spparms("spumoni") info if requested
2543  int typ = mattype.type ();
2544  mattype.info ();
2545 
2546  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2547  (*current_liboctave_error_handler) ("incorrect matrix type");
2548 
2549  double anorm = 0.;
2550  double ainvnorm = 0.;
2551  octave_idx_type b_nc = b.cols ();
2552  rcond = 1.;
2553 
2554  if (calc_cond)
2555  {
2556  // Calculate the 1-norm of matrix for rcond calculation
2557  for (octave_idx_type j = 0; j < nc; j++)
2558  {
2559  double atmp = 0.;
2560  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2561  atmp += std::abs (data (i));
2562  if (atmp > anorm)
2563  anorm = atmp;
2564  }
2565  }
2566 
2567  if (typ == MatrixType::Permuted_Lower)
2568  {
2569  retval.resize (nc, b_nc);
2570  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2571  octave_idx_type *perm = mattype.triangular_perm ();
2572 
2573  for (octave_idx_type j = 0; j < b_nc; j++)
2574  {
2575  for (octave_idx_type i = 0; i < nm; i++)
2576  work[i] = 0.;
2577  for (octave_idx_type i = 0; i < nr; i++)
2578  work[perm[i]] = b(i,j);
2579 
2580  for (octave_idx_type k = 0; k < nc; k++)
2581  {
2582  if (work[k] != 0.)
2583  {
2584  octave_idx_type minr = nr;
2585  octave_idx_type mini = 0;
2586 
2587  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2588  if (perm[ridx (i)] < minr)
2589  {
2590  minr = perm[ridx (i)];
2591  mini = i;
2592  }
2593 
2594  if (minr != k || data (mini) == 0.)
2595  {
2596  err = -2;
2597  goto triangular_error;
2598  }
2599 
2600  Complex tmp = work[k] / data (mini);
2601  work[k] = tmp;
2602  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2603  {
2604  if (i == mini)
2605  continue;
2606 
2607  octave_idx_type iidx = perm[ridx (i)];
2608  work[iidx] = work[iidx] - tmp * data (i);
2609  }
2610  }
2611  }
2612 
2613  for (octave_idx_type i = 0; i < nc; i++)
2614  retval(i, j) = work[i];
2615  }
2616 
2617  if (calc_cond)
2618  {
2619  // Calculation of 1-norm of inv(*this)
2620  for (octave_idx_type i = 0; i < nm; i++)
2621  work[i] = 0.;
2622 
2623  for (octave_idx_type j = 0; j < nr; j++)
2624  {
2625  work[j] = 1.;
2626 
2627  for (octave_idx_type k = 0; k < nc; k++)
2628  {
2629  if (work[k] != 0.)
2630  {
2631  octave_idx_type minr = nr;
2632  octave_idx_type mini = 0;
2633 
2634  for (octave_idx_type i = cidx (k);
2635  i < cidx (k+1); i++)
2636  if (perm[ridx (i)] < minr)
2637  {
2638  minr = perm[ridx (i)];
2639  mini = i;
2640  }
2641 
2642  Complex tmp = work[k] / data (mini);
2643  work[k] = tmp;
2644  for (octave_idx_type i = cidx (k);
2645  i < cidx (k+1); i++)
2646  {
2647  if (i == mini)
2648  continue;
2649 
2650  octave_idx_type iidx = perm[ridx (i)];
2651  work[iidx] = work[iidx] - tmp * data (i);
2652  }
2653  }
2654  }
2655 
2656  double atmp = 0;
2657  for (octave_idx_type i = j; i < nc; i++)
2658  {
2659  atmp += std::abs (work[i]);
2660  work[i] = 0.;
2661  }
2662  if (atmp > ainvnorm)
2663  ainvnorm = atmp;
2664  }
2665  rcond = 1. / ainvnorm / anorm;
2666  }
2667  }
2668  else
2669  {
2670  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2671  retval.resize (nc, b_nc, 0.);
2672 
2673  for (octave_idx_type j = 0; j < b_nc; j++)
2674  {
2675  for (octave_idx_type i = 0; i < nr; i++)
2676  work[i] = b(i,j);
2677  for (octave_idx_type i = nr; i < nc; i++)
2678  work[i] = 0.;
2679  for (octave_idx_type k = 0; k < nc; k++)
2680  {
2681  if (work[k] != 0.)
2682  {
2683  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2684  {
2685  err = -2;
2686  goto triangular_error;
2687  }
2688 
2689  Complex tmp = work[k] / data (cidx (k));
2690  work[k] = tmp;
2691  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2692  {
2693  octave_idx_type iidx = ridx (i);
2694  work[iidx] = work[iidx] - tmp * data (i);
2695  }
2696  }
2697  }
2698  for (octave_idx_type i = 0; i < nc; i++)
2699  retval.xelem (i, j) = work[i];
2700  }
2701 
2702  if (calc_cond)
2703  {
2704  // Calculation of 1-norm of inv(*this)
2705  for (octave_idx_type i = 0; i < nm; i++)
2706  work[i] = 0.;
2707 
2708  for (octave_idx_type j = 0; j < nr; j++)
2709  {
2710  work[j] = 1.;
2711 
2712  for (octave_idx_type k = j; k < nc; k++)
2713  {
2714 
2715  if (work[k] != 0.)
2716  {
2717  Complex tmp = work[k] / data (cidx (k));
2718  work[k] = tmp;
2719  for (octave_idx_type i = cidx (k)+1;
2720  i < cidx (k+1); i++)
2721  {
2722  octave_idx_type iidx = ridx (i);
2723  work[iidx] = work[iidx] - tmp * data (i);
2724  }
2725  }
2726  }
2727  double atmp = 0;
2728  for (octave_idx_type i = j; i < nc; i++)
2729  {
2730  atmp += std::abs (work[i]);
2731  work[i] = 0.;
2732  }
2733  if (atmp > ainvnorm)
2734  ainvnorm = atmp;
2735  }
2736  rcond = 1. / ainvnorm / anorm;
2737  }
2738  }
2739  triangular_error:
2740  if (err != 0)
2741  {
2742  if (sing_handler)
2743  {
2744  sing_handler (rcond);
2745  mattype.mark_as_rectangular ();
2746  }
2747  else
2749  }
2750 
2751  volatile double rcond_plus_one = rcond + 1.0;
2752 
2753  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2754  {
2755  err = -2;
2756 
2757  if (sing_handler)
2758  {
2759  sing_handler (rcond);
2760  mattype.mark_as_rectangular ();
2761  }
2762  else
2764  }
2765  }
2766 
2767  return retval;
2768 }
2769 
2772  octave_idx_type& err, double& rcond,
2773  solve_singularity_handler sing_handler,
2774  bool calc_cond) const
2775 {
2777 
2778  octave_idx_type nr = rows ();
2779  octave_idx_type nc = cols ();
2780  octave_idx_type nm = (nc > nr ? nc : nr);
2781 
2782  err = 0;
2783 
2784  if (nr != b.rows ())
2786  ("matrix dimension mismatch solution of linear equations");
2787 
2788  if (nr == 0 || nc == 0 || b.cols () == 0)
2789  retval = SparseComplexMatrix (nc, b.cols ());
2790  else
2791  {
2792  // Print spparms("spumoni") info if requested
2793  int typ = mattype.type ();
2794  mattype.info ();
2795 
2796  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2797  (*current_liboctave_error_handler) ("incorrect matrix type");
2798 
2799  double anorm = 0.;
2800  double ainvnorm = 0.;
2801  rcond = 1.;
2802 
2803  if (calc_cond)
2804  {
2805  // Calculate the 1-norm of matrix for rcond calculation
2806  for (octave_idx_type j = 0; j < nc; j++)
2807  {
2808  double atmp = 0.;
2809  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2810  atmp += std::abs (data (i));
2811  if (atmp > anorm)
2812  anorm = atmp;
2813  }
2814  }
2815 
2816  octave_idx_type b_nc = b.cols ();
2817  octave_idx_type b_nz = b.nnz ();
2818  retval = SparseComplexMatrix (nc, b_nc, b_nz);
2819  retval.xcidx (0) = 0;
2820  octave_idx_type ii = 0;
2821  octave_idx_type x_nz = b_nz;
2822 
2823  if (typ == MatrixType::Permuted_Lower)
2824  {
2825  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2826  octave_idx_type *perm = mattype.triangular_perm ();
2827 
2828  for (octave_idx_type j = 0; j < b_nc; j++)
2829  {
2830  for (octave_idx_type i = 0; i < nm; i++)
2831  work[i] = 0.;
2832  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2833  work[perm[b.ridx (i)]] = b.data (i);
2834 
2835  for (octave_idx_type k = 0; k < nc; k++)
2836  {
2837  if (work[k] != 0.)
2838  {
2839  octave_idx_type minr = nr;
2840  octave_idx_type mini = 0;
2841 
2842  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2843  if (perm[ridx (i)] < minr)
2844  {
2845  minr = perm[ridx (i)];
2846  mini = i;
2847  }
2848 
2849  if (minr != k || data (mini) == 0.)
2850  {
2851  err = -2;
2852  goto triangular_error;
2853  }
2854 
2855  Complex tmp = work[k] / data (mini);
2856  work[k] = tmp;
2857  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2858  {
2859  if (i == mini)
2860  continue;
2861 
2862  octave_idx_type iidx = perm[ridx (i)];
2863  work[iidx] = work[iidx] - tmp * data (i);
2864  }
2865  }
2866  }
2867 
2868  // Count nonzeros in work vector and adjust space in
2869  // retval if needed
2870  octave_idx_type new_nnz = 0;
2871  for (octave_idx_type i = 0; i < nc; i++)
2872  if (work[i] != 0.)
2873  new_nnz++;
2874 
2875  if (ii + new_nnz > x_nz)
2876  {
2877  // Resize the sparse matrix
2878  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2879  retval.change_capacity (sz);
2880  x_nz = sz;
2881  }
2882 
2883  for (octave_idx_type i = 0; i < nc; i++)
2884  if (work[i] != 0.)
2885  {
2886  retval.xridx (ii) = i;
2887  retval.xdata (ii++) = work[i];
2888  }
2889  retval.xcidx (j+1) = ii;
2890  }
2891 
2892  retval.maybe_compress ();
2893 
2894  if (calc_cond)
2895  {
2896  // Calculation of 1-norm of inv(*this)
2897  for (octave_idx_type i = 0; i < nm; i++)
2898  work[i] = 0.;
2899 
2900  for (octave_idx_type j = 0; j < nr; j++)
2901  {
2902  work[j] = 1.;
2903 
2904  for (octave_idx_type k = 0; k < nc; k++)
2905  {
2906  if (work[k] != 0.)
2907  {
2908  octave_idx_type minr = nr;
2909  octave_idx_type mini = 0;
2910 
2911  for (octave_idx_type i = cidx (k);
2912  i < cidx (k+1); i++)
2913  if (perm[ridx (i)] < minr)
2914  {
2915  minr = perm[ridx (i)];
2916  mini = i;
2917  }
2918 
2919  Complex tmp = work[k] / data (mini);
2920  work[k] = tmp;
2921  for (octave_idx_type i = cidx (k);
2922  i < cidx (k+1); i++)
2923  {
2924  if (i == mini)
2925  continue;
2926 
2927  octave_idx_type iidx = perm[ridx (i)];
2928  work[iidx] = work[iidx] - tmp * data (i);
2929  }
2930  }
2931  }
2932 
2933  double atmp = 0;
2934  for (octave_idx_type i = j; i < nc; i++)
2935  {
2936  atmp += std::abs (work[i]);
2937  work[i] = 0.;
2938  }
2939  if (atmp > ainvnorm)
2940  ainvnorm = atmp;
2941  }
2942  rcond = 1. / ainvnorm / anorm;
2943  }
2944  }
2945  else
2946  {
2947  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2948 
2949  for (octave_idx_type j = 0; j < b_nc; j++)
2950  {
2951  for (octave_idx_type i = 0; i < nm; i++)
2952  work[i] = 0.;
2953  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2954  work[b.ridx (i)] = b.data (i);
2955 
2956  for (octave_idx_type k = 0; k < nc; k++)
2957  {
2958  if (work[k] != 0.)
2959  {
2960  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2961  {
2962  err = -2;
2963  goto triangular_error;
2964  }
2965 
2966  Complex tmp = work[k] / data (cidx (k));
2967  work[k] = tmp;
2968  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2969  {
2970  octave_idx_type iidx = ridx (i);
2971  work[iidx] = work[iidx] - tmp * data (i);
2972  }
2973  }
2974  }
2975 
2976  // Count nonzeros in work vector and adjust space in
2977  // retval if needed
2978  octave_idx_type new_nnz = 0;
2979  for (octave_idx_type i = 0; i < nc; i++)
2980  if (work[i] != 0.)
2981  new_nnz++;
2982 
2983  if (ii + new_nnz > x_nz)
2984  {
2985  // Resize the sparse matrix
2986  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2987  retval.change_capacity (sz);
2988  x_nz = sz;
2989  }
2990 
2991  for (octave_idx_type i = 0; i < nc; i++)
2992  if (work[i] != 0.)
2993  {
2994  retval.xridx (ii) = i;
2995  retval.xdata (ii++) = work[i];
2996  }
2997  retval.xcidx (j+1) = ii;
2998  }
2999 
3000  retval.maybe_compress ();
3001 
3002  if (calc_cond)
3003  {
3004  // Calculation of 1-norm of inv(*this)
3005  for (octave_idx_type i = 0; i < nm; i++)
3006  work[i] = 0.;
3007 
3008  for (octave_idx_type j = 0; j < nr; j++)
3009  {
3010  work[j] = 1.;
3011 
3012  for (octave_idx_type k = j; k < nc; k++)
3013  {
3014 
3015  if (work[k] != 0.)
3016  {
3017  Complex tmp = work[k] / data (cidx (k));
3018  work[k] = tmp;
3019  for (octave_idx_type i = cidx (k)+1;
3020  i < cidx (k+1); i++)
3021  {
3022  octave_idx_type iidx = ridx (i);
3023  work[iidx] = work[iidx] - tmp * data (i);
3024  }
3025  }
3026  }
3027  double atmp = 0;
3028  for (octave_idx_type i = j; i < nc; i++)
3029  {
3030  atmp += std::abs (work[i]);
3031  work[i] = 0.;
3032  }
3033  if (atmp > ainvnorm)
3034  ainvnorm = atmp;
3035  }
3036  rcond = 1. / ainvnorm / anorm;
3037  }
3038  }
3039 
3040  triangular_error:
3041  if (err != 0)
3042  {
3043  if (sing_handler)
3044  {
3045  sing_handler (rcond);
3046  mattype.mark_as_rectangular ();
3047  }
3048  else
3050  }
3051 
3052  volatile double rcond_plus_one = rcond + 1.0;
3053 
3054  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
3055  {
3056  err = -2;
3057 
3058  if (sing_handler)
3059  {
3060  sing_handler (rcond);
3061  mattype.mark_as_rectangular ();
3062  }
3063  else
3065  }
3066  }
3067 
3068  return retval;
3069 }
3070 
3073  octave_idx_type& err, double& rcond,
3074  solve_singularity_handler sing_handler,
3075  bool calc_cond) const
3076 {
3078 
3079  octave_idx_type nr = rows ();
3080  octave_idx_type nc = cols ();
3081  octave_idx_type nm = (nc > nr ? nc : nr);
3082  err = 0;
3083 
3084  if (nr != b.rows ())
3086  ("matrix dimension mismatch solution of linear equations");
3087 
3088  if (nr == 0 || nc == 0 || b.cols () == 0)
3089  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3090  else
3091  {
3092  // Print spparms("spumoni") info if requested
3093  int typ = mattype.type ();
3094  mattype.info ();
3095 
3096  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3097  (*current_liboctave_error_handler) ("incorrect matrix type");
3098 
3099  double anorm = 0.;
3100  double ainvnorm = 0.;
3101  octave_idx_type b_nc = b.cols ();
3102  rcond = 1.;
3103 
3104  if (calc_cond)
3105  {
3106  // Calculate the 1-norm of matrix for rcond calculation
3107  for (octave_idx_type j = 0; j < nc; j++)
3108  {
3109  double atmp = 0.;
3110  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3111  atmp += std::abs (data (i));
3112  if (atmp > anorm)
3113  anorm = atmp;
3114  }
3115  }
3116 
3117  if (typ == MatrixType::Permuted_Lower)
3118  {
3119  retval.resize (nc, b_nc);
3120  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3121  octave_idx_type *perm = mattype.triangular_perm ();
3122 
3123  for (octave_idx_type j = 0; j < b_nc; j++)
3124  {
3125  for (octave_idx_type i = 0; i < nm; i++)
3126  work[i] = 0.;
3127  for (octave_idx_type i = 0; i < nr; i++)
3128  work[perm[i]] = b(i,j);
3129 
3130  for (octave_idx_type k = 0; k < nc; k++)
3131  {
3132  if (work[k] != 0.)
3133  {
3134  octave_idx_type minr = nr;
3135  octave_idx_type mini = 0;
3136 
3137  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3138  if (perm[ridx (i)] < minr)
3139  {
3140  minr = perm[ridx (i)];
3141  mini = i;
3142  }
3143 
3144  if (minr != k || data (mini) == 0.)
3145  {
3146  err = -2;
3147  goto triangular_error;
3148  }
3149 
3150  Complex tmp = work[k] / data (mini);
3151  work[k] = tmp;
3152  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3153  {
3154  if (i == mini)
3155  continue;
3156 
3157  octave_idx_type iidx = perm[ridx (i)];
3158  work[iidx] = work[iidx] - tmp * data (i);
3159  }
3160  }
3161  }
3162 
3163  for (octave_idx_type i = 0; i < nc; i++)
3164  retval(i, j) = work[i];
3165  }
3166 
3167  if (calc_cond)
3168  {
3169  // Calculation of 1-norm of inv(*this)
3170  for (octave_idx_type i = 0; i < nm; i++)
3171  work[i] = 0.;
3172 
3173  for (octave_idx_type j = 0; j < nr; j++)
3174  {
3175  work[j] = 1.;
3176 
3177  for (octave_idx_type k = 0; k < nc; k++)
3178  {
3179  if (work[k] != 0.)
3180  {
3181  octave_idx_type minr = nr;
3182  octave_idx_type mini = 0;
3183 
3184  for (octave_idx_type i = cidx (k);
3185  i < cidx (k+1); i++)
3186  if (perm[ridx (i)] < minr)
3187  {
3188  minr = perm[ridx (i)];
3189  mini = i;
3190  }
3191 
3192  Complex tmp = work[k] / data (mini);
3193  work[k] = tmp;
3194  for (octave_idx_type i = cidx (k);
3195  i < cidx (k+1); i++)
3196  {
3197  if (i == mini)
3198  continue;
3199 
3200  octave_idx_type iidx = perm[ridx (i)];
3201  work[iidx] = work[iidx] - tmp * data (i);
3202  }
3203  }
3204  }
3205 
3206  double atmp = 0;
3207  for (octave_idx_type i = j; i < nc; i++)
3208  {
3209  atmp += std::abs (work[i]);
3210  work[i] = 0.;
3211  }
3212  if (atmp > ainvnorm)
3213  ainvnorm = atmp;
3214  }
3215  rcond = 1. / ainvnorm / anorm;
3216  }
3217  }
3218  else
3219  {
3220  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3221  retval.resize (nc, b_nc, 0.);
3222 
3223  for (octave_idx_type j = 0; j < b_nc; j++)
3224  {
3225  for (octave_idx_type i = 0; i < nr; i++)
3226  work[i] = b(i,j);
3227  for (octave_idx_type i = nr; i < nc; i++)
3228  work[i] = 0.;
3229 
3230  for (octave_idx_type k = 0; k < nc; k++)
3231  {
3232  if (work[k] != 0.)
3233  {
3234  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3235  {
3236  err = -2;
3237  goto triangular_error;
3238  }
3239 
3240  Complex tmp = work[k] / data (cidx (k));
3241  work[k] = tmp;
3242  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3243  {
3244  octave_idx_type iidx = ridx (i);
3245  work[iidx] = work[iidx] - tmp * data (i);
3246  }
3247  }
3248  }
3249 
3250  for (octave_idx_type i = 0; i < nc; i++)
3251  retval.xelem (i, j) = work[i];
3252  }
3253 
3254  if (calc_cond)
3255  {
3256  // Calculation of 1-norm of inv(*this)
3257  for (octave_idx_type i = 0; i < nm; i++)
3258  work[i] = 0.;
3259 
3260  for (octave_idx_type j = 0; j < nr; j++)
3261  {
3262  work[j] = 1.;
3263 
3264  for (octave_idx_type k = j; k < nc; k++)
3265  {
3266 
3267  if (work[k] != 0.)
3268  {
3269  Complex tmp = work[k] / data (cidx (k));
3270  work[k] = tmp;
3271  for (octave_idx_type i = cidx (k)+1;
3272  i < cidx (k+1); i++)
3273  {
3274  octave_idx_type iidx = ridx (i);
3275  work[iidx] = work[iidx] - tmp * data (i);
3276  }
3277  }
3278  }
3279  double atmp = 0;
3280  for (octave_idx_type i = j; i < nc; i++)
3281  {
3282  atmp += std::abs (work[i]);
3283  work[i] = 0.;
3284  }
3285  if (atmp > ainvnorm)
3286  ainvnorm = atmp;
3287  }
3288  rcond = 1. / ainvnorm / anorm;
3289  }
3290  }
3291 
3292  triangular_error:
3293  if (err != 0)
3294  {
3295  if (sing_handler)
3296  {
3297  sing_handler (rcond);
3298  mattype.mark_as_rectangular ();
3299  }
3300  else
3302  }
3303 
3304  volatile double rcond_plus_one = rcond + 1.0;
3305 
3306  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
3307  {
3308  err = -2;
3309 
3310  if (sing_handler)
3311  {
3312  sing_handler (rcond);
3313  mattype.mark_as_rectangular ();
3314  }
3315  else
3317  }
3318  }
3319 
3320  return retval;
3321 }
3322 
3325  octave_idx_type& err, double& rcond,
3326  solve_singularity_handler sing_handler,
3327  bool calc_cond) const
3328 {
3330 
3331  octave_idx_type nr = rows ();
3332  octave_idx_type nc = cols ();
3333  octave_idx_type nm = (nc > nr ? nc : nr);
3334  err = 0;
3335 
3336  if (nr != b.rows ())
3338  ("matrix dimension mismatch solution of linear equations");
3339 
3340  if (nr == 0 || nc == 0 || b.cols () == 0)
3341  retval = SparseComplexMatrix (nc, b.cols ());
3342  else
3343  {
3344  // Print spparms("spumoni") info if requested
3345  int typ = mattype.type ();
3346  mattype.info ();
3347 
3348  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3349  (*current_liboctave_error_handler) ("incorrect matrix type");
3350 
3351  double anorm = 0.;
3352  double ainvnorm = 0.;
3353  rcond = 1.;
3354 
3355  if (calc_cond)
3356  {
3357  // Calculate the 1-norm of matrix for rcond calculation
3358  for (octave_idx_type j = 0; j < nc; j++)
3359  {
3360  double atmp = 0.;
3361  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3362  atmp += std::abs (data (i));
3363  if (atmp > anorm)
3364  anorm = atmp;
3365  }
3366  }
3367 
3368  octave_idx_type b_nc = b.cols ();
3369  octave_idx_type b_nz = b.nnz ();
3370  retval = SparseComplexMatrix (nc, b_nc, b_nz);
3371  retval.xcidx (0) = 0;
3372  octave_idx_type ii = 0;
3373  octave_idx_type x_nz = b_nz;
3374 
3375  if (typ == MatrixType::Permuted_Lower)
3376  {
3377  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3378  octave_idx_type *perm = mattype.triangular_perm ();
3379 
3380  for (octave_idx_type j = 0; j < b_nc; j++)
3381  {
3382  for (octave_idx_type i = 0; i < nm; i++)
3383  work[i] = 0.;
3384  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3385  work[perm[b.ridx (i)]] = b.data (i);
3386 
3387  for (octave_idx_type k = 0; k < nc; k++)
3388  {
3389  if (work[k] != 0.)
3390  {
3391  octave_idx_type minr = nr;
3392  octave_idx_type mini = 0;
3393 
3394  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3395  if (perm[ridx (i)] < minr)
3396  {
3397  minr = perm[ridx (i)];
3398  mini = i;
3399  }
3400 
3401  if (minr != k || data (mini) == 0.)
3402  {
3403  err = -2;
3404  goto triangular_error;
3405  }
3406 
3407  Complex tmp = work[k] / data (mini);
3408  work[k] = tmp;
3409  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3410  {
3411  if (i == mini)
3412  continue;
3413 
3414  octave_idx_type iidx = perm[ridx (i)];
3415  work[iidx] = work[iidx] - tmp * data (i);
3416  }
3417  }
3418  }
3419 
3420  // Count nonzeros in work vector and adjust space in
3421  // retval if needed
3422  octave_idx_type new_nnz = 0;
3423  for (octave_idx_type i = 0; i < nc; i++)
3424  if (work[i] != 0.)
3425  new_nnz++;
3426 
3427  if (ii + new_nnz > x_nz)
3428  {
3429  // Resize the sparse matrix
3430  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3431  retval.change_capacity (sz);
3432  x_nz = sz;
3433  }
3434 
3435  for (octave_idx_type i = 0; i < nc; i++)
3436  if (work[i] != 0.)
3437  {
3438  retval.xridx (ii) = i;
3439  retval.xdata (ii++) = work[i];
3440  }
3441  retval.xcidx (j+1) = ii;
3442  }
3443 
3444  retval.maybe_compress ();
3445 
3446  if (calc_cond)
3447  {
3448  // Calculation of 1-norm of inv(*this)
3449  for (octave_idx_type i = 0; i < nm; i++)
3450  work[i] = 0.;
3451 
3452  for (octave_idx_type j = 0; j < nr; j++)
3453  {
3454  work[j] = 1.;
3455 
3456  for (octave_idx_type k = 0; k < nc; k++)
3457  {
3458  if (work[k] != 0.)
3459  {
3460  octave_idx_type minr = nr;
3461  octave_idx_type mini = 0;
3462 
3463  for (octave_idx_type i = cidx (k);
3464  i < cidx (k+1); i++)
3465  if (perm[ridx (i)] < minr)
3466  {
3467  minr = perm[ridx (i)];
3468  mini = i;
3469  }
3470 
3471  Complex tmp = work[k] / data (mini);
3472  work[k] = tmp;
3473  for (octave_idx_type i = cidx (k);
3474  i < cidx (k+1); i++)
3475  {
3476  if (i == mini)
3477  continue;
3478 
3479  octave_idx_type iidx = perm[ridx (i)];
3480  work[iidx] = work[iidx] - tmp * data (i);
3481  }
3482  }
3483  }
3484 
3485  double atmp = 0;
3486  for (octave_idx_type i = j; i < nc; i++)
3487  {
3488  atmp += std::abs (work[i]);
3489  work[i] = 0.;
3490  }
3491  if (atmp > ainvnorm)
3492  ainvnorm = atmp;
3493  }
3494  rcond = 1. / ainvnorm / anorm;
3495  }
3496  }
3497  else
3498  {
3499  OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3500 
3501  for (octave_idx_type j = 0; j < b_nc; j++)
3502  {
3503  for (octave_idx_type i = 0; i < nm; i++)
3504  work[i] = 0.;
3505  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3506  work[b.ridx (i)] = b.data (i);
3507 
3508  for (octave_idx_type k = 0; k < nc; k++)
3509  {
3510  if (work[k] != 0.)
3511  {
3512  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3513  {
3514  err = -2;
3515  goto triangular_error;
3516  }
3517 
3518  Complex tmp = work[k] / data (cidx (k));
3519  work[k] = tmp;
3520  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3521  {
3522  octave_idx_type iidx = ridx (i);
3523  work[iidx] = work[iidx] - tmp * data (i);
3524  }
3525  }
3526  }
3527 
3528  // Count nonzeros in work vector and adjust space in
3529  // retval if needed
3530  octave_idx_type new_nnz = 0;
3531  for (octave_idx_type i = 0; i < nc; i++)
3532  if (work[i] != 0.)
3533  new_nnz++;
3534 
3535  if (ii + new_nnz > x_nz)
3536  {
3537  // Resize the sparse matrix
3538  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3539  retval.change_capacity (sz);
3540  x_nz = sz;
3541  }
3542 
3543  for (octave_idx_type i = 0; i < nc; i++)
3544  if (work[i] != 0.)
3545  {
3546  retval.xridx (ii) = i;
3547  retval.xdata (ii++) = work[i];
3548  }
3549  retval.xcidx (j+1) = ii;
3550  }
3551 
3552  retval.maybe_compress ();
3553 
3554  if (calc_cond)
3555  {
3556  // Calculation of 1-norm of inv(*this)
3557  for (octave_idx_type i = 0; i < nm; i++)
3558  work[i] = 0.;
3559 
3560  for (octave_idx_type j = 0; j < nr; j++)
3561  {
3562  work[j] = 1.;
3563 
3564  for (octave_idx_type k = j; k < nc; k++)
3565  {
3566 
3567  if (work[k] != 0.)
3568  {
3569  Complex tmp = work[k] / data (cidx (k));
3570  work[k] = tmp;
3571  for (octave_idx_type i = cidx (k)+1;
3572  i < cidx (k+1); i++)
3573  {
3574  octave_idx_type iidx = ridx (i);
3575  work[iidx] = work[iidx] - tmp * data (i);
3576  }
3577  }
3578  }
3579  double atmp = 0;
3580  for (octave_idx_type i = j; i < nc; i++)
3581  {
3582  atmp += std::abs (work[i]);
3583  work[i] = 0.;
3584  }
3585  if (atmp > ainvnorm)
3586  ainvnorm = atmp;
3587  }
3588  rcond = 1. / ainvnorm / anorm;
3589  }
3590  }
3591 
3592  triangular_error:
3593  if (err != 0)
3594  {
3595  if (sing_handler)
3596  {
3597  sing_handler (rcond);
3598  mattype.mark_as_rectangular ();
3599  }
3600  else
3602  }
3603 
3604  volatile double rcond_plus_one = rcond + 1.0;
3605 
3606  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
3607  {
3608  err = -2;
3609 
3610  if (sing_handler)
3611  {
3612  sing_handler (rcond);
3613  mattype.mark_as_rectangular ();
3614  }
3615  else
3617  }
3618  }
3619 
3620  return retval;
3621 }
3622 
3625  octave_idx_type& err, double& rcond,
3626  solve_singularity_handler sing_handler,
3627  bool calc_cond) const
3628 {
3630 
3631  octave_idx_type nr = rows ();
3632  octave_idx_type nc = cols ();
3633  err = 0;
3634 
3635  if (nr != nc || nr != b.rows ())
3637  ("matrix dimension mismatch solution of linear equations");
3638 
3639  if (nr == 0 || b.cols () == 0)
3640  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3641  else if (calc_cond)
3642  (*current_liboctave_error_handler)
3643  ("calculation of condition number not implemented");
3644  else
3645  {
3646  // Print spparms("spumoni") info if requested
3647  volatile int typ = mattype.type ();
3648  mattype.info ();
3649 
3651  {
3652  OCTAVE_LOCAL_BUFFER (double, D, nr);
3653  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3654 
3655  if (mattype.is_dense ())
3656  {
3657  octave_idx_type ii = 0;
3658 
3659  for (octave_idx_type j = 0; j < nc-1; j++)
3660  {
3661  D[j] = octave::math::real (data (ii++));
3662  DL[j] = data (ii);
3663  ii += 2;
3664  }
3665  D[nc-1] = octave::math::real (data (ii));
3666  }
3667  else
3668  {
3669  D[0] = 0.;
3670  for (octave_idx_type i = 0; i < nr - 1; i++)
3671  {
3672  D[i+1] = 0.;
3673  DL[i] = 0.;
3674  }
3675 
3676  for (octave_idx_type j = 0; j < nc; j++)
3677  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3678  {
3679  if (ridx (i) == j)
3680  D[j] = octave::math::real (data (i));
3681  else if (ridx (i) == j + 1)
3682  DL[j] = data (i);
3683  }
3684  }
3685 
3686  octave_idx_type b_nc = b.cols ();
3687  retval = ComplexMatrix (b);
3688  Complex *result = retval.fortran_vec ();
3689 
3690  F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL),
3691  F77_DBLE_CMPLX_ARG (result),
3692  b.rows (), err));
3693 
3694  if (err != 0)
3695  {
3696  err = 0;
3697  mattype.mark_as_unsymmetric ();
3699  }
3700  else
3701  rcond = 1.;
3702  }
3703 
3704  if (typ == MatrixType::Tridiagonal)
3705  {
3706  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3707  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3708  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3709 
3710  if (mattype.is_dense ())
3711  {
3712  octave_idx_type ii = 0;
3713 
3714  for (octave_idx_type j = 0; j < nc-1; j++)
3715  {
3716  D[j] = data (ii++);
3717  DL[j] = data (ii++);
3718  DU[j] = data (ii++);
3719  }
3720  D[nc-1] = data (ii);
3721  }
3722  else
3723  {
3724  D[0] = 0.;
3725  for (octave_idx_type i = 0; i < nr - 1; i++)
3726  {
3727  D[i+1] = 0.;
3728  DL[i] = 0.;
3729  DU[i] = 0.;
3730  }
3731 
3732  for (octave_idx_type j = 0; j < nc; j++)
3733  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3734  {
3735  if (ridx (i) == j)
3736  D[j] = data (i);
3737  else if (ridx (i) == j + 1)
3738  DL[j] = data (i);
3739  else if (ridx (i) == j - 1)
3740  DU[j-1] = data (i);
3741  }
3742  }
3743 
3744  octave_idx_type b_nc = b.cols ();
3745  retval = ComplexMatrix (b);
3746  Complex *result = retval.fortran_vec ();
3747 
3748  F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, F77_DBLE_CMPLX_ARG (DL),
3750  b.rows (), err));
3751 
3752  if (err != 0)
3753  {
3754  rcond = 0.;
3755  err = -2;
3756 
3757  if (sing_handler)
3758  {
3759  sing_handler (rcond);
3760  mattype.mark_as_rectangular ();
3761  }
3762  else
3764 
3765  }
3766  else
3767  rcond = 1.;
3768  }
3769  else if (typ != MatrixType::Tridiagonal_Hermitian)
3770  (*current_liboctave_error_handler) ("incorrect matrix type");
3771  }
3772 
3773  return retval;
3774 }
3775 
3778  octave_idx_type& err, double& rcond,
3779  solve_singularity_handler sing_handler,
3780  bool calc_cond) const
3781 {
3783 
3784  octave_idx_type nr = rows ();
3785  octave_idx_type nc = cols ();
3786  err = 0;
3787 
3788  if (nr != nc || nr != b.rows ())
3790  ("matrix dimension mismatch solution of linear equations");
3791 
3792  if (nr == 0 || b.cols () == 0)
3793  retval = SparseComplexMatrix (nc, b.cols ());
3794  else if (calc_cond)
3795  (*current_liboctave_error_handler)
3796  ("calculation of condition number not implemented");
3797  else
3798  {
3799  // Print spparms("spumoni") info if requested
3800  int typ = mattype.type ();
3801  mattype.info ();
3802 
3803  // Note can't treat symmetric case as there is no dpttrf function
3804  if (typ == MatrixType::Tridiagonal
3806  {
3807  OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
3808  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3809  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3810  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3811  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
3812  octave_idx_type *pipvt = ipvt.fortran_vec ();
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  F77_XFCN (zgttrf, ZGTTRF, (nr, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D),
3849  F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, err));
3850 
3851  if (err != 0)
3852  {
3853  err = -2;
3854  rcond = 0.0;
3855 
3856  if (sing_handler)
3857  {
3858  sing_handler (rcond);
3859  mattype.mark_as_rectangular ();
3860  }
3861  else
3863  }
3864  else
3865  {
3866  char job = 'N';
3867  volatile octave_idx_type x_nz = b.nnz ();
3868  octave_idx_type b_nc = b.cols ();
3869  retval = SparseComplexMatrix (nr, b_nc, x_nz);
3870  retval.xcidx (0) = 0;
3871  volatile octave_idx_type ii = 0;
3872  rcond = 1.0;
3873 
3874  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
3875 
3876  for (volatile octave_idx_type j = 0; j < b_nc; j++)
3877  {
3878  for (octave_idx_type i = 0; i < nr; i++)
3879  work[i] = 0.;
3880  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3881  work[b.ridx (i)] = b.data (i);
3882 
3883  F77_XFCN (zgttrs, ZGTTRS,
3884  (F77_CONST_CHAR_ARG2 (&job, 1),
3886  F77_DBLE_CMPLX_ARG (DU2), pipvt,
3887  F77_DBLE_CMPLX_ARG (work), b.rows (), err
3888  F77_CHAR_ARG_LEN (1)));
3889 
3890  // Count nonzeros in work vector and adjust
3891  // space in retval if needed
3892  octave_idx_type new_nnz = 0;
3893  for (octave_idx_type i = 0; i < nr; i++)
3894  if (work[i] != 0.)
3895  new_nnz++;
3896 
3897  if (ii + new_nnz > x_nz)
3898  {
3899  // Resize the sparse matrix
3900  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3901  retval.change_capacity (sz);
3902  x_nz = sz;
3903  }
3904 
3905  for (octave_idx_type i = 0; i < nr; i++)
3906  if (work[i] != 0.)
3907  {
3908  retval.xridx (ii) = i;
3909  retval.xdata (ii++) = work[i];
3910  }
3911  retval.xcidx (j+1) = ii;
3912  }
3913 
3914  retval.maybe_compress ();
3915  }
3916  }
3917  else if (typ != MatrixType::Tridiagonal_Hermitian)
3918  (*current_liboctave_error_handler) ("incorrect matrix type");
3919  }
3920 
3921  return retval;
3922 }
3923 
3926  octave_idx_type& err, double& rcond,
3927  solve_singularity_handler sing_handler,
3928  bool calc_cond) const
3929 {
3931 
3932  octave_idx_type nr = rows ();
3933  octave_idx_type nc = cols ();
3934  err = 0;
3935 
3936  if (nr != nc || nr != b.rows ())
3938  ("matrix dimension mismatch solution of linear equations");
3939 
3940  if (nr == 0 || b.cols () == 0)
3941  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3942  else if (calc_cond)
3943  (*current_liboctave_error_handler)
3944  ("calculation of condition number not implemented");
3945  else
3946  {
3947  // Print spparms("spumoni") info if requested
3948  volatile int typ = mattype.type ();
3949  mattype.info ();
3950 
3952  {
3953  OCTAVE_LOCAL_BUFFER (double, D, nr);
3954  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3955 
3956  if (mattype.is_dense ())
3957  {
3958  octave_idx_type ii = 0;
3959 
3960  for (octave_idx_type j = 0; j < nc-1; j++)
3961  {
3962  D[j] = octave::math::real (data (ii++));
3963  DL[j] = data (ii);
3964  ii += 2;
3965  }
3966  D[nc-1] = octave::math::real (data (ii));
3967  }
3968  else
3969  {
3970  D[0] = 0.;
3971  for (octave_idx_type i = 0; i < nr - 1; i++)
3972  {
3973  D[i+1] = 0.;
3974  DL[i] = 0.;
3975  }
3976 
3977  for (octave_idx_type j = 0; j < nc; j++)
3978  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3979  {
3980  if (ridx (i) == j)
3981  D[j] = octave::math::real (data (i));
3982  else if (ridx (i) == j + 1)
3983  DL[j] = data (i);
3984  }
3985  }
3986 
3987  octave_idx_type b_nr = b.rows ();
3988  octave_idx_type b_nc = b.cols ();
3989  rcond = 1.;
3990 
3991  retval = ComplexMatrix (b);
3992  Complex *result = retval.fortran_vec ();
3993 
3994  F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL),
3995  F77_DBLE_CMPLX_ARG (result),
3996  b_nr, err));
3997 
3998  if (err != 0)
3999  {
4000  err = 0;
4001  mattype.mark_as_unsymmetric ();
4003  }
4004  }
4005 
4006  if (typ == MatrixType::Tridiagonal)
4007  {
4008  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4009  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4010  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4011 
4012  if (mattype.is_dense ())
4013  {
4014  octave_idx_type ii = 0;
4015 
4016  for (octave_idx_type j = 0; j < nc-1; j++)
4017  {
4018  D[j] = data (ii++);
4019  DL[j] = data (ii++);
4020  DU[j] = data (ii++);
4021  }
4022  D[nc-1] = data (ii);
4023  }
4024  else
4025  {
4026  D[0] = 0.;
4027  for (octave_idx_type i = 0; i < nr - 1; i++)
4028  {
4029  D[i+1] = 0.;
4030  DL[i] = 0.;
4031  DU[i] = 0.;
4032  }
4033 
4034  for (octave_idx_type j = 0; j < nc; j++)
4035  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4036  {
4037  if (ridx (i) == j)
4038  D[j] = data (i);
4039  else if (ridx (i) == j + 1)
4040  DL[j] = data (i);
4041  else if (ridx (i) == j - 1)
4042  DU[j-1] = data (i);
4043  }
4044  }
4045 
4046  octave_idx_type b_nr = b.rows ();
4047  octave_idx_type b_nc = b.cols ();
4048  rcond = 1.;
4049 
4050  retval = ComplexMatrix (b);
4051  Complex *result = retval.fortran_vec ();
4052 
4053  F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, F77_DBLE_CMPLX_ARG (DL),
4055  b_nr, err));
4056 
4057  if (err != 0)
4058  {
4059  rcond = 0.;
4060  err = -2;
4061 
4062  if (sing_handler)
4063  {
4064  sing_handler (rcond);
4065  mattype.mark_as_rectangular ();
4066  }
4067  else
4069  }
4070  }
4071  else if (typ != MatrixType::Tridiagonal_Hermitian)
4072  (*current_liboctave_error_handler) ("incorrect matrix type");
4073  }
4074 
4075  return retval;
4076 }
4077 
4080  const SparseComplexMatrix& b,
4081  octave_idx_type& err, double& rcond,
4082  solve_singularity_handler sing_handler,
4083  bool calc_cond) const
4084 {
4086 
4087  octave_idx_type nr = rows ();
4088  octave_idx_type nc = cols ();
4089  err = 0;
4090 
4091  if (nr != nc || nr != b.rows ())
4093  ("matrix dimension mismatch solution of linear equations");
4094 
4095  if (nr == 0 || b.cols () == 0)
4096  retval = SparseComplexMatrix (nc, b.cols ());
4097  else if (calc_cond)
4098  (*current_liboctave_error_handler)
4099  ("calculation of condition number not implemented");
4100  else
4101  {
4102  // Print spparms("spumoni") info if requested
4103  int typ = mattype.type ();
4104  mattype.info ();
4105 
4106  // Note can't treat symmetric case as there is no dpttrf function
4107  if (typ == MatrixType::Tridiagonal
4109  {
4110  OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
4111  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4112  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4113  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4114  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4115  octave_idx_type *pipvt = ipvt.fortran_vec ();
4116 
4117  if (mattype.is_dense ())
4118  {
4119  octave_idx_type ii = 0;
4120 
4121  for (octave_idx_type j = 0; j < nc-1; j++)
4122  {
4123  D[j] = data (ii++);
4124  DL[j] = data (ii++);
4125  DU[j] = data (ii++);
4126  }
4127  D[nc-1] = data (ii);
4128  }
4129  else
4130  {
4131  D[0] = 0.;
4132  for (octave_idx_type i = 0; i < nr - 1; i++)
4133  {
4134  D[i+1] = 0.;
4135  DL[i] = 0.;
4136  DU[i] = 0.;
4137  }
4138 
4139  for (octave_idx_type j = 0; j < nc; j++)
4140  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4141  {
4142  if (ridx (i) == j)
4143  D[j] = data (i);
4144  else if (ridx (i) == j + 1)
4145  DL[j] = data (i);
4146  else if (ridx (i) == j - 1)
4147  DU[j-1] = data (i);
4148  }
4149  }
4150 
4151  F77_XFCN (zgttrf, ZGTTRF, (nr, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D),
4152  F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, err));
4153 
4154  if (err != 0)
4155  {
4156  rcond = 0.0;
4157  err = -2;
4158 
4159  if (sing_handler)
4160  {
4161  sing_handler (rcond);
4162  mattype.mark_as_rectangular ();
4163  }
4164  else
4166  }
4167  else
4168  {
4169  rcond = 1.;
4170  char job = 'N';
4171  octave_idx_type b_nr = b.rows ();
4172  octave_idx_type b_nc = b.cols ();
4173  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4174 
4175  // Take a first guess that the number of nonzero terms
4176  // will be as many as in b
4177  volatile octave_idx_type x_nz = b.nnz ();
4178  volatile octave_idx_type ii = 0;
4179  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4180 
4181  retval.xcidx (0) = 0;
4182  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4183  {
4184 
4185  for (octave_idx_type i = 0; i < b_nr; i++)
4186  Bx[i] = b(i,j);
4187 
4188  F77_XFCN (zgttrs, ZGTTRS,
4189  (F77_CONST_CHAR_ARG2 (&job, 1),
4191  F77_DBLE_CMPLX_ARG (DU2), pipvt,
4192  F77_DBLE_CMPLX_ARG (Bx), b_nr, err
4193  F77_CHAR_ARG_LEN (1)));
4194 
4195  if (err != 0)
4196  {
4197  // FIXME: This should probably be a warning so that
4198  // error value can be passed back.
4199  (*current_liboctave_error_handler)
4200  ("SparseComplexMatrix::solve solve failed");
4201 
4202  err = -1;
4203  break;
4204  }
4205 
4206  // Count nonzeros in work vector and adjust
4207  // space in retval if needed
4208  octave_idx_type new_nnz = 0;
4209  for (octave_idx_type i = 0; i < nr; i++)
4210  if (Bx[i] != 0.)
4211  new_nnz++;
4212 
4213  if (ii + new_nnz > x_nz)
4214  {
4215  // Resize the sparse matrix
4216  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4217  retval.change_capacity (sz);
4218  x_nz = sz;
4219  }
4220 
4221  for (octave_idx_type i = 0; i < nr; i++)
4222  if (Bx[i] != 0.)
4223  {
4224  retval.xridx (ii) = i;
4225  retval.xdata (ii++) = Bx[i];
4226  }
4227 
4228  retval.xcidx (j+1) = ii;
4229  }
4230 
4231  retval.maybe_compress ();
4232  }
4233  }
4234  else if (typ != MatrixType::Tridiagonal_Hermitian)
4235  (*current_liboctave_error_handler) ("incorrect matrix type");
4236  }
4237 
4238  return retval;
4239 }
4240 
4243  octave_idx_type& err, double& rcond,
4244  solve_singularity_handler sing_handler,
4245  bool calc_cond) const
4246 {
4248 
4249  octave_idx_type nr = rows ();
4250  octave_idx_type nc = cols ();
4251  err = 0;
4252 
4253  if (nr != nc || nr != b.rows ())
4255  ("matrix dimension mismatch solution of linear equations");
4256 
4257  if (nr == 0 || b.cols () == 0)
4258  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4259  else
4260  {
4261  // Print spparms("spumoni") info if requested
4262  volatile int typ = mattype.type ();
4263  mattype.info ();
4264 
4265  if (typ == MatrixType::Banded_Hermitian)
4266  {
4267  octave_idx_type n_lower = mattype.nlower ();
4268  octave_idx_type ldm = n_lower + 1;
4269  ComplexMatrix m_band (ldm, nc);
4270  Complex *tmp_data = m_band.fortran_vec ();
4271 
4272  if (! mattype.is_dense ())
4273  {
4274  octave_idx_type ii = 0;
4275 
4276  for (octave_idx_type j = 0; j < ldm; j++)
4277  for (octave_idx_type i = 0; i < nc; i++)
4278  tmp_data[ii++] = 0.;
4279  }
4280 
4281  for (octave_idx_type j = 0; j < nc; j++)
4282  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4283  {
4284  octave_idx_type ri = ridx (i);
4285  if (ri >= j)
4286  m_band(ri - j, j) = data (i);
4287  }
4288 
4289  // Calculate the norm of the matrix, for later use.
4290  double anorm;
4291  if (calc_cond)
4292  anorm = m_band.abs ().sum ().row (0).max ();
4293 
4294  char job = 'L';
4295  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4296  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err
4297  F77_CHAR_ARG_LEN (1)));
4298 
4299  if (err != 0)
4300  {
4301  rcond = 0.0;
4302  // Matrix is not positive definite!! Fall through to
4303  // unsymmetric banded solver.
4304  mattype.mark_as_unsymmetric ();
4305  typ = MatrixType::Banded;
4306  err = 0;
4307  }
4308  else
4309  {
4310  if (calc_cond)
4311  {
4312  Array<Complex> z (dim_vector (2 * nr, 1));
4313  Complex *pz = z.fortran_vec ();
4314  Array<double> iz (dim_vector (nr, 1));
4315  double *piz = iz.fortran_vec ();
4316 
4317  F77_XFCN (zpbcon, ZPBCON,
4318  (F77_CONST_CHAR_ARG2 (&job, 1),
4319  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm,
4320  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
4321  F77_CHAR_ARG_LEN (1)));
4322 
4323  if (err != 0)
4324  err = -2;
4325 
4326  volatile double rcond_plus_one = rcond + 1.0;
4327 
4328  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4329  {
4330  err = -2;
4331 
4332  if (sing_handler)
4333  {
4334  sing_handler (rcond);
4335  mattype.mark_as_rectangular ();
4336  }
4337  else
4339  }
4340  }
4341  else
4342  rcond = 1.0;
4343 
4344  if (err == 0)
4345  {
4346  retval = ComplexMatrix (b);
4347  Complex *result = retval.fortran_vec ();
4348 
4349  octave_idx_type b_nc = b.cols ();
4350 
4351  F77_XFCN (zpbtrs, ZPBTRS,
4352  (F77_CONST_CHAR_ARG2 (&job, 1),
4353  nr, n_lower, b_nc, F77_DBLE_CMPLX_ARG (tmp_data),
4354  ldm, F77_DBLE_CMPLX_ARG (result), b.rows (), err
4355  F77_CHAR_ARG_LEN (1)));
4356 
4357  if (err != 0)
4358  {
4359  // FIXME: Probably should be a warning.
4360  (*current_liboctave_error_handler)
4361  ("SparseMatrix::solve solve failed");
4362  err = -1;
4363  }
4364  }
4365  }
4366  }
4367 
4368  if (typ == MatrixType::Banded)
4369  {
4370  // Create the storage for the banded form of the sparse matrix
4371  octave_idx_type n_upper = mattype.nupper ();
4372  octave_idx_type n_lower = mattype.nlower ();
4373  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4374 
4375  ComplexMatrix m_band (ldm, nc);
4376  Complex *tmp_data = m_band.fortran_vec ();
4377 
4378  if (! mattype.is_dense ())
4379  {
4380  octave_idx_type ii = 0;
4381 
4382  for (octave_idx_type j = 0; j < ldm; j++)
4383  for (octave_idx_type i = 0; i < nc; i++)
4384  tmp_data[ii++] = 0.;
4385  }
4386 
4387  for (octave_idx_type j = 0; j < nc; j++)
4388  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4389  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4390 
4391  // Calculate the norm of the matrix, for later use.
4392  double anorm = 0.0;
4393  if (calc_cond)
4394  {
4395  for (octave_idx_type j = 0; j < nr; j++)
4396  {
4397  double atmp = 0.;
4398  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4399  atmp += std::abs (data (i));
4400  if (atmp > anorm)
4401  anorm = atmp;
4402  }
4403  }
4404 
4405  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4406  octave_idx_type *pipvt = ipvt.fortran_vec ();
4407 
4408  F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper,
4409  F77_DBLE_CMPLX_ARG (tmp_data),
4410  ldm, pipvt, err));
4411 
4412  // Throw-away extra info LAPACK gives so as to not
4413  // change output.
4414  if (err != 0)
4415  {
4416  rcond = 0.0;
4417  err = -2;
4418 
4419  if (sing_handler)
4420  {
4421  sing_handler (rcond);
4422  mattype.mark_as_rectangular ();
4423  }
4424  else
4426  }
4427  else
4428  {
4429  if (calc_cond)
4430  {
4431  char job = '1';
4432  Array<Complex> z (dim_vector (2 * nr, 1));
4433  Complex *pz = z.fortran_vec ();
4434  Array<double> iz (dim_vector (nr, 1));
4435  double *piz = iz.fortran_vec ();
4436 
4437  F77_XFCN (zgbcon, ZGBCON,
4438  (F77_CONST_CHAR_ARG2 (&job, 1),
4439  nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt,
4440  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
4441  F77_CHAR_ARG_LEN (1)));
4442 
4443  if (err != 0)
4444  err = -2;
4445 
4446  volatile double rcond_plus_one = rcond + 1.0;
4447 
4448  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4449  {
4450  err = -2;
4451 
4452  if (sing_handler)
4453  {
4454  sing_handler (rcond);
4455  mattype.mark_as_rectangular ();
4456  }
4457  else
4459  }
4460  }
4461  else
4462  rcond = 1.;
4463 
4464  if (err == 0)
4465  {
4466  retval = ComplexMatrix (b);
4467  Complex *result = retval.fortran_vec ();
4468 
4469  octave_idx_type b_nc = b.cols ();
4470 
4471  char job = 'N';
4472  F77_XFCN (zgbtrs, ZGBTRS,
4473  (F77_CONST_CHAR_ARG2 (&job, 1),
4474  nr, n_lower, n_upper, b_nc, F77_DBLE_CMPLX_ARG (tmp_data),
4475  ldm, pipvt, F77_DBLE_CMPLX_ARG (result), b.rows (), err
4476  F77_CHAR_ARG_LEN (1)));
4477  }
4478  }
4479  }
4480  else if (typ != MatrixType::Banded_Hermitian)
4481  (*current_liboctave_error_handler) ("incorrect matrix type");
4482  }
4483 
4484  return retval;
4485 }
4486 
4489  octave_idx_type& err, double& rcond,
4490  solve_singularity_handler sing_handler,
4491  bool calc_cond) const
4492 {
4494 
4495  octave_idx_type nr = rows ();
4496  octave_idx_type nc = cols ();
4497  err = 0;
4498 
4499  if (nr != nc || nr != b.rows ())
4501  ("matrix dimension mismatch solution of linear equations");
4502 
4503  if (nr == 0 || b.cols () == 0)
4504  retval = SparseComplexMatrix (nc, b.cols ());
4505  else
4506  {
4507  // Print spparms("spumoni") info if requested
4508  volatile int typ = mattype.type ();
4509  mattype.info ();
4510 
4511  if (typ == MatrixType::Banded_Hermitian)
4512  {
4513  octave_idx_type n_lower = mattype.nlower ();
4514  octave_idx_type ldm = n_lower + 1;
4515 
4516  ComplexMatrix m_band (ldm, nc);
4517  Complex *tmp_data = m_band.fortran_vec ();
4518 
4519  if (! mattype.is_dense ())
4520  {
4521  octave_idx_type ii = 0;
4522 
4523  for (octave_idx_type j = 0; j < ldm; j++)
4524  for (octave_idx_type i = 0; i < nc; i++)
4525  tmp_data[ii++] = 0.;
4526  }
4527 
4528  for (octave_idx_type j = 0; j < nc; j++)
4529  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4530  {
4531  octave_idx_type ri = ridx (i);
4532  if (ri >= j)
4533  m_band(ri - j, j) = data (i);
4534  }
4535 
4536  // Calculate the norm of the matrix, for later use.
4537  double anorm;
4538  if (calc_cond)
4539  anorm = m_band.abs ().sum ().row (0).max ();
4540 
4541  char job = 'L';
4542  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4543  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err
4544  F77_CHAR_ARG_LEN (1)));
4545 
4546  if (err != 0)
4547  {
4548  rcond = 0.0;
4549  mattype.mark_as_unsymmetric ();
4550  typ = MatrixType::Banded;
4551  err = 0;
4552  }
4553  else
4554  {
4555  if (calc_cond)
4556  {
4557  Array<Complex> z (dim_vector (2 * nr, 1));
4558  Complex *pz = z.fortran_vec ();
4559  Array<double> iz (dim_vector (nr, 1));
4560  double *piz = iz.fortran_vec ();
4561 
4562  F77_XFCN (zpbcon, ZPBCON,
4563  (F77_CONST_CHAR_ARG2 (&job, 1),
4564  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm,
4565  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
4566  F77_CHAR_ARG_LEN (1)));
4567 
4568  if (err != 0)
4569  err = -2;
4570 
4571  volatile double rcond_plus_one = rcond + 1.0;
4572 
4573  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4574  {
4575  err = -2;
4576 
4577  if (sing_handler)
4578  {
4579  sing_handler (rcond);
4580  mattype.mark_as_rectangular ();
4581  }
4582  else
4584  }
4585  }
4586  else
4587  rcond = 1.0;
4588 
4589  if (err == 0)
4590  {
4591  octave_idx_type b_nr = b.rows ();
4592  octave_idx_type b_nc = b.cols ();
4593  OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4594 
4595  // Take a first guess that the number of nonzero terms
4596  // will be as many as in b
4597  volatile octave_idx_type x_nz = b.nnz ();
4598  volatile octave_idx_type ii = 0;
4599  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4600 
4601  retval.xcidx (0) = 0;
4602  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4603  {
4604  for (octave_idx_type i = 0; i < b_nr; i++)
4605  Bx[i] = b.elem (i, j);
4606 
4607  F77_XFCN (zpbtrs, ZPBTRS,
4608  (F77_CONST_CHAR_ARG2 (&job, 1),
4609  nr, n_lower, 1, F77_DBLE_CMPLX_ARG (tmp_data),
4610  ldm, F77_DBLE_CMPLX_ARG (Bx), b_nr, err
4611  F77_CHAR_ARG_LEN (1)));
4612 
4613  if (err != 0)
4614  {
4615  // FIXME: Probably should be a warning.
4616  (*current_liboctave_error_handler)
4617  ("SparseComplexMatrix::solve solve failed");
4618  err = -1;
4619  break;
4620  }
4621 
4622  for (octave_idx_type i = 0; i < b_nr; i++)
4623  {
4624  Complex tmp = Bx[i];
4625  if (tmp != 0.0)
4626  {
4627  if (ii == x_nz)
4628  {
4629  // Resize the sparse matrix
4630  octave_idx_type sz = x_nz *
4631  (b_nc - j) / b_nc;
4632  sz = (sz > 10 ? sz : 10) + x_nz;
4633  retval.change_capacity (sz);
4634  x_nz = sz;
4635  }
4636  retval.xdata (ii) = tmp;
4637  retval.xridx (ii++) = i;
4638  }
4639  }
4640  retval.xcidx (j+1) = ii;
4641  }
4642 
4643  retval.maybe_compress ();
4644  }
4645  }
4646  }
4647 
4648  if (typ == MatrixType::Banded)
4649  {
4650  // Create the storage for the banded form of the sparse matrix
4651  octave_idx_type n_upper = mattype.nupper ();
4652  octave_idx_type n_lower = mattype.nlower ();
4653  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4654 
4655  ComplexMatrix m_band (ldm, nc);
4656  Complex *tmp_data = m_band.fortran_vec ();
4657 
4658  if (! mattype.is_dense ())
4659  {
4660  octave_idx_type ii = 0;
4661 
4662  for (octave_idx_type j = 0; j < ldm; j++)
4663  for (octave_idx_type i = 0; i < nc; i++)
4664  tmp_data[ii++] = 0.;
4665  }
4666 
4667  for (octave_idx_type j = 0; j < nc; j++)
4668  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4669  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4670 
4671  // Calculate the norm of the matrix, for later use.
4672  double anorm = 0.0;
4673  if (calc_cond)
4674  {
4675  for (octave_idx_type j = 0; j < nr; j++)
4676  {
4677  double atmp = 0.;
4678  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4679  atmp += std::abs (data (i));
4680  if (atmp > anorm)
4681  anorm = atmp;
4682  }
4683  }
4684 
4685  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4686  octave_idx_type *pipvt = ipvt.fortran_vec ();
4687 
4688  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper,
4689  F77_DBLE_CMPLX_ARG (tmp_data),
4690  ldm, pipvt, err));
4691 
4692  if (err != 0)
4693  {
4694  rcond = 0.0;
4695  err = -2;
4696 
4697  if (sing_handler)
4698  {
4699  sing_handler (rcond);
4700  mattype.mark_as_rectangular ();
4701  }
4702  else
4704  }
4705  else
4706  {
4707  if (calc_cond)
4708  {
4709  char job = '1';
4710  Array<Complex> z (dim_vector (2 * nr, 1));
4711  Complex *pz = z.fortran_vec ();
4712  Array<double> iz (dim_vector (nr, 1));
4713  double *piz = iz.fortran_vec ();
4714 
4715  F77_XFCN (zgbcon, ZGBCON,
4716  (F77_CONST_CHAR_ARG2 (&job, 1),
4717  nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt,
4718  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
4719  F77_CHAR_ARG_LEN (1)));
4720 
4721  if (err != 0)
4722  err = -2;
4723 
4724  volatile double rcond_plus_one = rcond + 1.0;
4725 
4726  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4727  {
4728  err = -2;
4729 
4730  if (sing_handler)
4731  {
4732  sing_handler (rcond);
4733  mattype.mark_as_rectangular ();
4734  }
4735  else
4737  }
4738  }
4739  else
4740  rcond = 1.;
4741 
4742  if (err == 0)
4743  {
4744  char job = 'N';
4745  volatile octave_idx_type x_nz = b.nnz ();
4746  octave_idx_type b_nc = b.cols ();
4747  retval = SparseComplexMatrix (nr, b_nc, x_nz);
4748  retval.xcidx (0) = 0;
4749  volatile octave_idx_type ii = 0;
4750 
4751  OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4752 
4753  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4754  {
4755  for (octave_idx_type i = 0; i < nr; i++)
4756  work[i] = 0.;
4757  for (octave_idx_type i = b.cidx (j);
4758  i < b.cidx (j+1); i++)
4759  work[b.ridx (i)] = b.data (i);
4760 
4761  F77_XFCN (zgbtrs, ZGBTRS,
4762  (F77_CONST_CHAR_ARG2 (&job, 1),
4763  nr, n_lower, n_upper, 1, F77_DBLE_CMPLX_ARG (tmp_data),
4764  ldm, pipvt, F77_DBLE_CMPLX_ARG (work), b.rows (), err
4765  F77_CHAR_ARG_LEN (1)));
4766 
4767  // Count nonzeros in work vector and adjust
4768  // space in retval if needed
4769  octave_idx_type new_nnz = 0;
4770  for (octave_idx_type i = 0; i < nr; i++)
4771  if (work[i] != 0.)
4772  new_nnz++;
4773 
4774  if (ii + new_nnz > x_nz)
4775  {
4776  // Resize the sparse matrix
4777  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4778  retval.change_capacity (sz);
4779  x_nz = sz;
4780  }
4781 
4782  for (octave_idx_type i = 0; i < nr; i++)
4783  if (work[i] != 0.)
4784  {
4785  retval.xridx (ii) = i;
4786  retval.xdata (ii++) = work[i];
4787  }
4788  retval.xcidx (j+1) = ii;
4789  }
4790 
4791  retval.maybe_compress ();
4792  }
4793  }
4794  }
4795  else if (typ != MatrixType::Banded_Hermitian)
4796  (*current_liboctave_error_handler) ("incorrect matrix type");
4797  }
4798 
4799  return retval;
4800 }
4801 
4804  octave_idx_type& err, double& rcond,
4805  solve_singularity_handler sing_handler,
4806  bool calc_cond) const
4807 {
4809 
4810  octave_idx_type nr = rows ();
4811  octave_idx_type nc = cols ();
4812  err = 0;
4813 
4814  if (nr != nc || nr != b.rows ())
4816  ("matrix dimension mismatch solution of linear equations");
4817 
4818  if (nr == 0 || b.cols () == 0)
4819  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4820  else
4821  {
4822  // Print spparms("spumoni") info if requested
4823  volatile int typ = mattype.type ();
4824  mattype.info ();
4825 
4826  if (typ == MatrixType::Banded_Hermitian)
4827  {
4828  octave_idx_type n_lower = mattype.nlower ();
4829  octave_idx_type ldm = n_lower + 1;
4830 
4831  ComplexMatrix m_band (ldm, nc);
4832  Complex *tmp_data = m_band.fortran_vec ();
4833 
4834  if (! mattype.is_dense ())
4835  {
4836  octave_idx_type ii = 0;
4837 
4838  for (octave_idx_type j = 0; j < ldm; j++)
4839  for (octave_idx_type i = 0; i < nc; i++)
4840  tmp_data[ii++] = 0.;
4841  }
4842 
4843  for (octave_idx_type j = 0; j < nc; j++)
4844  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4845  {
4846  octave_idx_type ri = ridx (i);
4847  if (ri >= j)
4848  m_band(ri - j, j) = data (i);
4849  }
4850 
4851  // Calculate the norm of the matrix, for later use.
4852  double anorm;
4853  if (calc_cond)
4854  anorm = m_band.abs ().sum ().row (0).max ();
4855 
4856  char job = 'L';
4857  F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4858  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err
4859  F77_CHAR_ARG_LEN (1)));
4860 
4861  if (err != 0)
4862  {
4863  // Matrix is not positive definite!! Fall through to
4864  // unsymmetric banded solver.
4865  rcond = 0.0;
4866  mattype.mark_as_unsymmetric ();
4867  typ = MatrixType::Banded;
4868  err = 0;
4869  }
4870  else
4871  {
4872  if (calc_cond)
4873  {
4874  Array<Complex> z (dim_vector (2 * nr, 1));
4875  Complex *pz = z.fortran_vec ();
4876  Array<double> iz (dim_vector (nr, 1));
4877  double *piz = iz.fortran_vec ();
4878 
4879  F77_XFCN (zpbcon, ZPBCON,
4880  (F77_CONST_CHAR_ARG2 (&job, 1),
4881  nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm,
4882  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
4883  F77_CHAR_ARG_LEN (1)));
4884 
4885  if (err != 0)
4886  err = -2;
4887 
4888  volatile double rcond_plus_one = rcond + 1.0;
4889 
4890  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4891  {
4892  err = -2;
4893 
4894  if (sing_handler)
4895  {
4896  sing_handler (rcond);
4897  mattype.mark_as_rectangular ();
4898  }
4899  else
4901  }
4902  }
4903  else
4904  rcond = 1.0;
4905 
4906  if (err == 0)
4907  {
4908  octave_idx_type b_nr = b.rows ();
4909  octave_idx_type b_nc = b.cols ();
4910  retval = ComplexMatrix (b);
4911  Complex *result = retval.fortran_vec ();
4912 
4913  F77_XFCN (zpbtrs, ZPBTRS,
4914  (F77_CONST_CHAR_ARG2 (&job, 1),
4915  nr, n_lower, b_nc, F77_DBLE_CMPLX_ARG (tmp_data),
4916  ldm, F77_DBLE_CMPLX_ARG (result), b_nr, err
4917  F77_CHAR_ARG_LEN (1)));
4918 
4919  if (err != 0)
4920  {
4921  // FIXME: Probably should be a warning.
4922  (*current_liboctave_error_handler)
4923  ("SparseComplexMatrix::solve solve failed");
4924  err = -1;
4925  }
4926  }
4927  }
4928  }
4929 
4930  if (typ == MatrixType::Banded)
4931  {
4932  // Create the storage for the banded form of the sparse matrix
4933  octave_idx_type n_upper = mattype.nupper ();
4934  octave_idx_type n_lower = mattype.nlower ();
4935  octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4936 
4937  ComplexMatrix m_band (ldm, nc);
4938  Complex *tmp_data = m_band.fortran_vec ();
4939 
4940  if (! mattype.is_dense ())
4941  {
4942  octave_idx_type ii = 0;
4943 
4944  for (octave_idx_type j = 0; j < ldm; j++)
4945  for (octave_idx_type i = 0; i < nc; i++)
4946  tmp_data[ii++] = 0.;
4947  }
4948 
4949  for (octave_idx_type j = 0; j < nc; j++)
4950  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4951  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4952 
4953  // Calculate the norm of the matrix, for later use.
4954  double anorm = 0.0;
4955  if (calc_cond)
4956  {
4957  for (octave_idx_type j = 0; j < nr; j++)
4958  {
4959  double atmp = 0.;
4960  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4961  atmp += std::abs (data (i));
4962  if (atmp > anorm)
4963  anorm = atmp;
4964  }
4965  }
4966 
4967  Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4968  octave_idx_type *pipvt = ipvt.fortran_vec ();
4969 
4970  F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper,
4971  F77_DBLE_CMPLX_ARG (tmp_data),
4972  ldm, pipvt, err));
4973 
4974  if (err != 0)
4975  {
4976  err = -2;
4977  rcond = 0.0;
4978 
4979  if (sing_handler)
4980  {
4981  sing_handler (rcond);
4982  mattype.mark_as_rectangular ();
4983  }
4984  else
4986  }
4987  else
4988  {
4989  if (calc_cond)
4990  {
4991  char job = '1';
4992  Array<Complex> z (dim_vector (2 * nr, 1));
4993  Complex *pz = z.fortran_vec ();
4994  Array<double> iz (dim_vector (nr, 1));
4995  double *piz = iz.fortran_vec ();
4996 
4997  F77_XFCN (zgbcon, ZGBCON,
4998  (F77_CONST_CHAR_ARG2 (&job, 1),
4999  nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt,
5000  anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err
5001  F77_CHAR_ARG_LEN (1)));
5002 
5003  if (err != 0)
5004  err = -2;
5005 
5006  volatile double rcond_plus_one = rcond + 1.0;