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