GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
CSparse.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 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