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