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