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_ca