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