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