GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dSparse.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2018 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
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13 
14 Octave is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License 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 <https://www.gnu.org/licenses/>.
22 
23 */
24 
25 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28 
29 #include <iostream>
30 
31 #include "quit.h"
32 #include "lo-ieee.h"
33 #include "lo-lapack-proto.h"
34 #include "lo-mappers.h"
35 #include "dRowVector.h"
36 #include "oct-locbuf.h"
37 
38 #include "dDiagMatrix.h"
39 #include "CSparse.h"
40 #include "boolSparse.h"
41 #include "dSparse.h"
42 #include "functor.h"
43 #include "oct-spparms.h"
44 #include "sparse-lu.h"
45 #include "MatrixType.h"
46 #include "oct-sparse.h"
47 #include "sparse-util.h"
48 #include "sparse-chol.h"
49 #include "sparse-qr.h"
50 
51 #include "Sparse-op-defs.h"
52 
53 #include "Sparse-diag-op-defs.h"
54 
55 #include "Sparse-perm-op-defs.h"
56 
57 // Define whether to use a basic QR solver or one that uses a Dulmange
58 // Mendelsohn factorization to seperate the problem into under-determined,
59 // well-determined and over-determined parts and solves them seperately
60 #if ! defined (USE_QRSOLVE)
61 # include "sparse-dmsolve.h"
62 #endif
63 
65  : MSparse<double> (a.rows (), a.cols (), a.nnz ())
66 {
67  octave_idx_type nc = cols ();
68  octave_idx_type nz = a.nnz ();
69 
70  for (octave_idx_type i = 0; i < nc + 1; i++)
71  cidx (i) = a.cidx (i);
72 
73  for (octave_idx_type i = 0; i < nz; i++)
74  {
75  data (i) = a.data (i);
76  ridx (i) = a.ridx (i);
77  }
78 }
79 
81  : MSparse<double> (a.rows (), a.cols (), a.length ())
82 {
83  octave_idx_type j = 0;
84  octave_idx_type l = a.length ();
85  for (octave_idx_type i = 0; i < l; i++)
86  {
87  cidx (i) = j;
88  if (a(i, i) != 0.0)
89  {
90  data (j) = a(i, i);
91  ridx (j) = i;
92  j++;
93  }
94  }
95  for (octave_idx_type i = l; i <= a.cols (); i++)
96  cidx (i) = j;
97 }
98 
99 bool
101 {
102  octave_idx_type nr = rows ();
103  octave_idx_type nc = cols ();
104  octave_idx_type nz = nnz ();
105  octave_idx_type nr_a = a.rows ();
106  octave_idx_type nc_a = a.cols ();
107  octave_idx_type nz_a = a.nnz ();
108 
109  if (nr != nr_a || nc != nc_a || nz != nz_a)
110  return false;
111 
112  for (octave_idx_type i = 0; i < nc + 1; i++)
113  if (cidx (i) != a.cidx (i))
114  return false;
115 
116  for (octave_idx_type i = 0; i < nz; i++)
117  if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
118  return false;
119 
120  return true;
121 }
122 
123 bool
125 {
126  return !(*this == a);
127 }
128 
129 bool
131 {
132  octave_idx_type nr = rows ();
133  octave_idx_type nc = cols ();
134 
135  if (nr == nc && nr > 0)
136  {
137  for (octave_idx_type j = 0; j < nc; j++)
138  {
139  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
140  {
141  octave_idx_type ri = ridx (i);
142 
143  if (ri != j)
144  {
145  bool found = false;
146 
147  for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
148  {
149  if (ridx (k) == j)
150  {
151  if (data (i) == data (k))
152  found = true;
153  break;
154  }
155  }
156 
157  if (! found)
158  return false;
159  }
160  }
161  }
162 
163  return true;
164  }
165 
166  return false;
167 }
168 
172 {
174  return *this;
175 }
176 
179 {
180  MSparse<double>::insert (a, indx);
181  return *this;
182 }
183 
185 SparseMatrix::max (int dim) const
186 {
187  Array<octave_idx_type> dummy_idx;
188  return max (dummy_idx, dim);
189 }
190 
192 SparseMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
193 {
195  dim_vector dv = dims ();
196  octave_idx_type nr = dv(0);
197  octave_idx_type nc = dv(1);
198 
199  if (dim >= dv.ndims ())
200  {
201  idx_arg.resize (dim_vector (nr, nc), 0);
202  return *this;
203  }
204 
205  if (dim < 0)
206  dim = dv.first_non_singleton ();
207 
208  if (dim == 0)
209  {
210  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
211 
212  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
213  return SparseMatrix (nr == 0 ? 0 : 1, nc);
214 
215  octave_idx_type nel = 0;
216  for (octave_idx_type j = 0; j < nc; j++)
217  {
218  double tmp_max = octave::numeric_limits<double>::NaN ();
219  octave_idx_type idx_j = 0;
220  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
221  {
222  if (ridx (i) != idx_j)
223  break;
224  else
225  idx_j++;
226  }
227 
228  if (idx_j != nr)
229  tmp_max = 0.;
230 
231  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
232  {
233  double tmp = data (i);
234 
235  if (octave::math::isnan (tmp))
236  continue;
237  else if (octave::math::isnan (tmp_max) || tmp > tmp_max)
238  {
239  idx_j = ridx (i);
240  tmp_max = tmp;
241  }
242 
243  }
244 
245  idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
246  if (tmp_max != 0.)
247  nel++;
248  }
249 
250  result = SparseMatrix (1, nc, nel);
251 
252  octave_idx_type ii = 0;
253  result.xcidx (0) = 0;
254  for (octave_idx_type j = 0; j < nc; j++)
255  {
256  double tmp = elem (idx_arg(j), j);
257  if (tmp != 0.)
258  {
259  result.xdata (ii) = tmp;
260  result.xridx (ii++) = 0;
261  }
262  result.xcidx (j+1) = ii;
263 
264  }
265  }
266  else
267  {
268  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
269 
270  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
271  return SparseMatrix (nr, nc == 0 ? 0 : 1);
272 
274 
275  for (octave_idx_type i = 0; i < nr; i++)
276  found[i] = 0;
277 
278  for (octave_idx_type j = 0; j < nc; j++)
279  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
280  if (found[ridx (i)] == -j)
281  found[ridx (i)] = -j - 1;
282 
283  for (octave_idx_type i = 0; i < nr; i++)
284  if (found[i] > -nc && found[i] < 0)
285  idx_arg.elem (i) = -found[i];
286 
287  for (octave_idx_type j = 0; j < nc; j++)
288  {
289  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
290  {
291  octave_idx_type ir = ridx (i);
292  octave_idx_type ix = idx_arg.elem (ir);
293  double tmp = data (i);
294 
295  if (octave::math::isnan (tmp))
296  continue;
297  else if (ix == -1 || tmp > elem (ir, ix))
298  idx_arg.elem (ir) = j;
299  }
300  }
301 
302  octave_idx_type nel = 0;
303  for (octave_idx_type j = 0; j < nr; j++)
304  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
305  nel++;
306 
307  result = SparseMatrix (nr, 1, nel);
308 
309  octave_idx_type ii = 0;
310  result.xcidx (0) = 0;
311  result.xcidx (1) = nel;
312  for (octave_idx_type j = 0; j < nr; j++)
313  {
314  if (idx_arg(j) == -1)
315  {
316  idx_arg(j) = 0;
318  result.xridx (ii++) = j;
319  }
320  else
321  {
322  double tmp = elem (j, idx_arg(j));
323  if (tmp != 0.)
324  {
325  result.xdata (ii) = tmp;
326  result.xridx (ii++) = j;
327  }
328  }
329  }
330  }
331 
332  return result;
333 }
334 
336 SparseMatrix::min (int dim) const
337 {
338  Array<octave_idx_type> dummy_idx;
339  return min (dummy_idx, dim);
340 }
341 
343 SparseMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
344 {
346  dim_vector dv = dims ();
347  octave_idx_type nr = dv(0);
348  octave_idx_type nc = dv(1);
349 
350  if (dim >= dv.ndims ())
351  {
352  idx_arg.resize (dim_vector (nr, nc), 0);
353  return *this;
354  }
355 
356  if (dim < 0)
357  dim = dv.first_non_singleton ();
358 
359  if (dim == 0)
360  {
361  idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
362 
363  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
364  return SparseMatrix (nr == 0 ? 0 : 1, nc);
365 
366  octave_idx_type nel = 0;
367  for (octave_idx_type j = 0; j < nc; j++)
368  {
369  double tmp_min = octave::numeric_limits<double>::NaN ();
370  octave_idx_type idx_j = 0;
371  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
372  {
373  if (ridx (i) != idx_j)
374  break;
375  else
376  idx_j++;
377  }
378 
379  if (idx_j != nr)
380  tmp_min = 0.;
381 
382  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
383  {
384  double tmp = data (i);
385 
386  if (octave::math::isnan (tmp))
387  continue;
388  else if (octave::math::isnan (tmp_min) || tmp < tmp_min)
389  {
390  idx_j = ridx (i);
391  tmp_min = tmp;
392  }
393 
394  }
395 
396  idx_arg.elem (j) = (octave::math::isnan (tmp_min) ? 0 : idx_j);
397  if (tmp_min != 0.)
398  nel++;
399  }
400 
401  result = SparseMatrix (1, nc, nel);
402 
403  octave_idx_type ii = 0;
404  result.xcidx (0) = 0;
405  for (octave_idx_type j = 0; j < nc; j++)
406  {
407  double tmp = elem (idx_arg(j), j);
408  if (tmp != 0.)
409  {
410  result.xdata (ii) = tmp;
411  result.xridx (ii++) = 0;
412  }
413  result.xcidx (j+1) = ii;
414 
415  }
416  }
417  else
418  {
419  idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
420 
421  if (nr == 0 || nc == 0 || dim >= dv.ndims ())
422  return SparseMatrix (nr, nc == 0 ? 0 : 1);
423 
425 
426  for (octave_idx_type i = 0; i < nr; i++)
427  found[i] = 0;
428 
429  for (octave_idx_type j = 0; j < nc; j++)
430  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
431  if (found[ridx (i)] == -j)
432  found[ridx (i)] = -j - 1;
433 
434  for (octave_idx_type i = 0; i < nr; i++)
435  if (found[i] > -nc && found[i] < 0)
436  idx_arg.elem (i) = -found[i];
437 
438  for (octave_idx_type j = 0; j < nc; j++)
439  {
440  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
441  {
442  octave_idx_type ir = ridx (i);
443  octave_idx_type ix = idx_arg.elem (ir);
444  double tmp = data (i);
445 
446  if (octave::math::isnan (tmp))
447  continue;
448  else if (ix == -1 || tmp < elem (ir, ix))
449  idx_arg.elem (ir) = j;
450  }
451  }
452 
453  octave_idx_type nel = 0;
454  for (octave_idx_type j = 0; j < nr; j++)
455  if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
456  nel++;
457 
458  result = SparseMatrix (nr, 1, nel);
459 
460  octave_idx_type ii = 0;
461  result.xcidx (0) = 0;
462  result.xcidx (1) = nel;
463  for (octave_idx_type j = 0; j < nr; j++)
464  {
465  if (idx_arg(j) == -1)
466  {
467  idx_arg(j) = 0;
469  result.xridx (ii++) = j;
470  }
471  else
472  {
473  double tmp = elem (j, idx_arg(j));
474  if (tmp != 0.)
475  {
476  result.xdata (ii) = tmp;
477  result.xridx (ii++) = j;
478  }
479  }
480  }
481  }
482 
483  return result;
484 }
485 
486 /*
487 
488 %!assert (max (max (speye (65536))), sparse (1))
489 %!assert (min (min (speye (65536))), sparse (0))
490 %!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
491 %!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
492 %!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
493 %!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
494 %!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
495 %!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
496 %!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
497 %!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
498 
499 */
500 
501 RowVector
503 {
504  octave_idx_type nc = columns ();
505  RowVector retval (nc, 0);
506 
507  for (octave_idx_type j = 0; j < nc; j++)
508  for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
509  {
510  if (ridx (k) == i)
511  {
512  retval(j) = data (k);
513  break;
514  }
515  }
516 
517  return retval;
518 }
519 
522 {
523  octave_idx_type nr = rows ();
524  ColumnVector retval (nr, 0);
525 
526  for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
527  retval(ridx (k)) = data (k);
528 
529  return retval;
530 }
531 
535 {
536  // Don't use numel to avoid all possiblity of an overflow
537  if (rb.rows () > 0 && rb.cols () > 0)
538  insert (rb, ra_idx(0), ra_idx(1));
539  return *this;
540 }
541 
545 {
546  SparseComplexMatrix retval (*this);
547  if (rb.rows () > 0 && rb.cols () > 0)
548  retval.insert (rb, ra_idx(0), ra_idx(1));
549  return retval;
550 }
551 
554 {
555  octave_idx_type nr = a.rows ();
556  octave_idx_type nc = a.cols ();
557  octave_idx_type nz = a.nnz ();
558  SparseMatrix r (nr, nc, nz);
559 
560  for (octave_idx_type i = 0; i < nc +1; i++)
561  r.cidx (i) = a.cidx (i);
562 
563  for (octave_idx_type i = 0; i < nz; i++)
564  {
565  r.data (i) = std::real (a.data (i));
566  r.ridx (i) = a.ridx (i);
567  }
568 
569  r.maybe_compress (true);
570  return r;
571 }
572 
575 {
576  octave_idx_type nr = a.rows ();
577  octave_idx_type nc = a.cols ();
578  octave_idx_type nz = a.nnz ();
579  SparseMatrix r (nr, nc, nz);
580 
581  for (octave_idx_type i = 0; i < nc +1; i++)
582  r.cidx (i) = a.cidx (i);
583 
584  for (octave_idx_type i = 0; i < nz; i++)
585  {
586  r.data (i) = std::imag (a.data (i));
587  r.ridx (i) = a.ridx (i);
588  }
589 
590  r.maybe_compress (true);
591  return r;
592 }
593 
594 /*
595 
596 %!assert (nnz (real (sparse ([1i,1]))), 1)
597 %!assert (nnz (real (sparse ([1i,1]))), 1)
598 
599 */
600 
603 {
604  octave_idx_type info;
605  double rcond;
606  MatrixType mattype (*this);
607  return inverse (mattype, info, rcond, 0, 0);
608 }
609 
612 {
613  octave_idx_type info;
614  double rcond;
615  return inverse (mattype, info, rcond, 0, 0);
616 }
617 
620 {
621  double rcond;
622  return inverse (mattype, info, rcond, 0, 0);
623 }
624 
627  double& rcond, const bool,
628  const bool calccond) const
629 {
631 
632  octave_idx_type nr = rows ();
633  octave_idx_type nc = cols ();
634  info = 0;
635 
636  if (nr == 0 || nc == 0 || nr != nc)
637  (*current_liboctave_error_handler) ("inverse requires square matrix");
638 
639  // Print spparms("spumoni") info if requested
640  int typ = mattype.type ();
641  mattype.info ();
642 
644  (*current_liboctave_error_handler) ("incorrect matrix type");
645 
647  retval = transpose ();
648  else
649  retval = *this;
650 
651  // Force make_unique to be called
652  double *v = retval.data ();
653 
654  if (calccond)
655  {
656  double dmax = 0.;
657  double dmin = octave::numeric_limits<double>::Inf ();
658  for (octave_idx_type i = 0; i < nr; i++)
659  {
660  double tmp = fabs (v[i]);
661  if (tmp > dmax)
662  dmax = tmp;
663  if (tmp < dmin)
664  dmin = tmp;
665  }
666  rcond = dmin / dmax;
667  }
668 
669  for (octave_idx_type i = 0; i < nr; i++)
670  v[i] = 1.0 / v[i];
671 
672  return retval;
673 }
674 
677  double& rcond, const bool,
678  const bool calccond) const
679 {
681 
682  octave_idx_type nr = rows ();
683  octave_idx_type nc = cols ();
684  info = 0;
685 
686  if (nr == 0 || nc == 0 || nr != nc)
687  (*current_liboctave_error_handler) ("inverse requires square matrix");
688 
689  // Print spparms("spumoni") info if requested
690  int typ = mattype.type ();
691  mattype.info ();
692 
693  if (typ != MatrixType::Upper && typ != MatrixType::Permuted_Upper
694  && typ != MatrixType::Lower && typ != MatrixType::Permuted_Lower)
695  (*current_liboctave_error_handler) ("incorrect matrix type");
696 
697  double anorm = 0.;
698  double ainvnorm = 0.;
699 
700  if (calccond)
701  {
702  // Calculate the 1-norm of matrix for rcond calculation
703  for (octave_idx_type j = 0; j < nr; j++)
704  {
705  double atmp = 0.;
706  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
707  atmp += fabs (data (i));
708  if (atmp > anorm)
709  anorm = atmp;
710  }
711  }
712 
713  if (typ == MatrixType::Upper || typ == MatrixType::Lower)
714  {
715  octave_idx_type nz = nnz ();
716  octave_idx_type cx = 0;
717  octave_idx_type nz2 = nz;
718  retval = SparseMatrix (nr, nc, nz2);
719 
720  for (octave_idx_type i = 0; i < nr; i++)
721  {
722  octave_quit ();
723  // place the 1 in the identity position
724  octave_idx_type cx_colstart = cx;
725 
726  if (cx == nz2)
727  {
728  nz2 *= 2;
729  retval.change_capacity (nz2);
730  }
731 
732  retval.xcidx (i) = cx;
733  retval.xridx (cx) = i;
734  retval.xdata (cx) = 1.0;
735  cx++;
736 
737  // iterate accross columns of input matrix
738  for (octave_idx_type j = i+1; j < nr; j++)
739  {
740  double v = 0.;
741  // iterate to calculate sum
742  octave_idx_type colXp = retval.xcidx (i);
743  octave_idx_type colUp = cidx (j);
744  octave_idx_type rpX, rpU;
745 
746  if (cidx (j) == cidx (j+1))
747  (*current_liboctave_error_handler) ("division by zero");
748 
749  do
750  {
751  octave_quit ();
752  rpX = retval.xridx (colXp);
753  rpU = ridx (colUp);
754 
755  if (rpX < rpU)
756  colXp++;
757  else if (rpX > rpU)
758  colUp++;
759  else
760  {
761  v -= retval.xdata (colXp) * data (colUp);
762  colXp++;
763  colUp++;
764  }
765  }
766  while (rpX < j && rpU < j && colXp < cx && colUp < nz);
767 
768  // get A(m,m)
769  if (typ == MatrixType::Upper)
770  colUp = cidx (j+1) - 1;
771  else
772  colUp = cidx (j);
773  double pivot = data (colUp);
774  if (pivot == 0. || ridx (colUp) != j)
775  (*current_liboctave_error_handler) ("division by zero");
776 
777  if (v != 0.)
778  {
779  if (cx == nz2)
780  {
781  nz2 *= 2;
782  retval.change_capacity (nz2);
783  }
784 
785  retval.xridx (cx) = j;
786  retval.xdata (cx) = v / pivot;
787  cx++;
788  }
789  }
790 
791  // get A(m,m)
792  octave_idx_type colUp;
793  if (typ == MatrixType::Upper)
794  colUp = cidx (i+1) - 1;
795  else
796  colUp = cidx (i);
797  double pivot = data (colUp);
798  if (pivot == 0. || ridx (colUp) != i)
799  (*current_liboctave_error_handler) ("division by zero");
800 
801  if (pivot != 1.0)
802  for (octave_idx_type j = cx_colstart; j < cx; j++)
803  retval.xdata (j) /= pivot;
804  }
805  retval.xcidx (nr) = cx;
806  retval.maybe_compress ();
807  }
808  else
809  {
810  octave_idx_type nz = nnz ();
811  octave_idx_type cx = 0;
812  octave_idx_type nz2 = nz;
813  retval = SparseMatrix (nr, nc, nz2);
814 
815  OCTAVE_LOCAL_BUFFER (double, work, nr);
817 
818  octave_idx_type *perm = mattype.triangular_perm ();
819  if (typ == MatrixType::Permuted_Upper)
820  {
821  for (octave_idx_type i = 0; i < nr; i++)
822  rperm[perm[i]] = i;
823  }
824  else
825  {
826  for (octave_idx_type i = 0; i < nr; i++)
827  rperm[i] = perm[i];
828  for (octave_idx_type i = 0; i < nr; i++)
829  perm[rperm[i]] = i;
830  }
831 
832  for (octave_idx_type i = 0; i < nr; i++)
833  {
834  octave_quit ();
835  octave_idx_type iidx = rperm[i];
836 
837  for (octave_idx_type j = 0; j < nr; j++)
838  work[j] = 0.;
839 
840  // place the 1 in the identity position
841  work[iidx] = 1.0;
842 
843  // iterate accross columns of input matrix
844  for (octave_idx_type j = iidx+1; j < nr; j++)
845  {
846  double v = 0.;
847  octave_idx_type jidx = perm[j];
848  // iterate to calculate sum
849  for (octave_idx_type k = cidx (jidx);
850  k < cidx (jidx+1); k++)
851  {
852  octave_quit ();
853  v -= work[ridx (k)] * data (k);
854  }
855 
856  // get A(m,m)
857  double pivot;
858  if (typ == MatrixType::Permuted_Upper)
859  pivot = data (cidx (jidx+1) - 1);
860  else
861  pivot = data (cidx (jidx));
862  if (pivot == 0.)
863  (*current_liboctave_error_handler) ("division by zero");
864 
865  work[j] = v / pivot;
866  }
867 
868  // get A(m,m)
869  octave_idx_type colUp;
870  if (typ == MatrixType::Permuted_Upper)
871  colUp = cidx (perm[iidx]+1) - 1;
872  else
873  colUp = cidx (perm[iidx]);
874 
875  double pivot = data (colUp);
876  if (pivot == 0.)
877  (*current_liboctave_error_handler) ("division by zero");
878 
879  octave_idx_type new_cx = cx;
880  for (octave_idx_type j = iidx; j < nr; j++)
881  if (work[j] != 0.0)
882  {
883  new_cx++;
884  if (pivot != 1.0)
885  work[j] /= pivot;
886  }
887 
888  if (cx < new_cx)
889  {
890  nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
891  retval.change_capacity (nz2);
892  }
893 
894  retval.xcidx (i) = cx;
895  for (octave_idx_type j = iidx; j < nr; j++)
896  if (work[j] != 0.)
897  {
898  retval.xridx (cx) = j;
899  retval.xdata (cx++) = work[j];
900  }
901  }
902 
903  retval.xcidx (nr) = cx;
904  retval.maybe_compress ();
905  }
906 
907  if (calccond)
908  {
909  // Calculate the 1-norm of inverse matrix for rcond calculation
910  for (octave_idx_type j = 0; j < nr; j++)
911  {
912  double atmp = 0.;
913  for (octave_idx_type i = retval.cidx (j);
914  i < retval.cidx (j+1); i++)
915  atmp += fabs (retval.data (i));
916  if (atmp > ainvnorm)
917  ainvnorm = atmp;
918  }
919 
920  rcond = 1. / ainvnorm / anorm;
921  }
922 
923  return retval;
924 }
925 
928  double& rcond, bool, bool calc_cond) const
929 {
930  int typ = mattype.type (false);
931  SparseMatrix ret;
932 
933  if (typ == MatrixType::Unknown)
934  typ = mattype.type (*this);
935 
937  ret = dinverse (mattype, info, rcond, true, calc_cond);
938  else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
939  ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
940  else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
941  {
942  MatrixType newtype = mattype.transpose ();
943  ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
944  }
945  else
946  {
947  if (mattype.ishermitian ())
948  {
949  MatrixType tmp_typ (MatrixType::Upper);
950  octave::math::sparse_chol<SparseMatrix> fact (*this, info, false);
951  rcond = fact.rcond ();
952  if (info == 0)
953  {
954  double rcond2;
955  SparseMatrix Q = fact.Q ();
956  SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
957  info, rcond2, true, false);
958  ret = Q * InvL.transpose () * InvL * Q.transpose ();
959  }
960  else
961  {
962  // Matrix is either singular or not positive definite
963  mattype.mark_as_unsymmetric ();
964  }
965  }
966 
967  if (! mattype.ishermitian ())
968  {
969  octave_idx_type n = rows ();
970  ColumnVector Qinit(n);
971  for (octave_idx_type i = 0; i < n; i++)
972  Qinit(i) = i;
973 
974  MatrixType tmp_typ (MatrixType::Upper);
976  Qinit, Matrix (),
977  false, false);
978  rcond = fact.rcond ();
979  double rcond2;
980  SparseMatrix InvL = fact.L ().transpose ().tinverse (tmp_typ,
981  info, rcond2, true, false);
982  SparseMatrix InvU = fact.U ().tinverse (tmp_typ, info, rcond2,
983  true, false).transpose ();
984  ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
985  }
986  }
987 
988  return ret;
989 }
990 
991 DET
993 {
994  octave_idx_type info;
995  double rcond;
996  return determinant (info, rcond, 0);
997 }
998 
999 DET
1001 {
1002  double rcond;
1003  return determinant (info, rcond, 0);
1004 }
1005 
1006 DET
1007 SparseMatrix::determinant (octave_idx_type& err, double& rcond, bool) const
1008 {
1009  DET retval;
1010 
1011 #if defined (HAVE_UMFPACK)
1012 
1013  octave_idx_type nr = rows ();
1014  octave_idx_type nc = cols ();
1015 
1016  if (nr == 0 || nc == 0 || nr != nc)
1017  {
1018  retval = DET (1.0);
1019  }
1020  else
1021  {
1022  err = 0;
1023 
1024  // Setup the control parameters
1025  Matrix Control (UMFPACK_CONTROL, 1);
1026  double *control = Control.fortran_vec ();
1027  UMFPACK_DNAME (defaults) (control);
1028 
1029  double tmp = octave_sparse_params::get_key ("spumoni");
1030  if (! octave::math::isnan (tmp))
1031  Control (UMFPACK_PRL) = tmp;
1032 
1033  tmp = octave_sparse_params::get_key ("piv_tol");
1034  if (! octave::math::isnan (tmp))
1035  {
1036  Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1037  Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1038  }
1039 
1040  // Set whether we are allowed to modify Q or not
1041  tmp = octave_sparse_params::get_key ("autoamd");
1042  if (! octave::math::isnan (tmp))
1043  Control (UMFPACK_FIXQ) = tmp;
1044 
1045  // Turn-off UMFPACK scaling for LU
1046  Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1047 
1048  UMFPACK_DNAME (report_control) (control);
1049 
1050  const octave_idx_type *Ap = cidx ();
1051  const octave_idx_type *Ai = ridx ();
1052  const double *Ax = data ();
1053 
1054  UMFPACK_DNAME (report_matrix) (nr, nc,
1057  Ax, 1, control);
1058 
1059  void *Symbolic;
1060  Matrix Info (1, UMFPACK_INFO);
1061  double *info = Info.fortran_vec ();
1062  int status = UMFPACK_DNAME (qsymbolic) (nr, nc,
1065  Ax, nullptr, &Symbolic, control, info);
1066 
1067  if (status < 0)
1068  {
1069  UMFPACK_DNAME (report_status) (control, status);
1070  UMFPACK_DNAME (report_info) (control, info);
1071 
1072  UMFPACK_DNAME (free_symbolic) (&Symbolic);
1073 
1074  (*current_liboctave_error_handler)
1075  ("SparseMatrix::determinant symbolic factorization failed");
1076  }
1077  else
1078  {
1079  UMFPACK_DNAME (report_symbolic) (Symbolic, control);
1080 
1081  void *Numeric;
1082  status = UMFPACK_DNAME (numeric) (octave::to_suitesparse_intptr (Ap),
1084  Ax, Symbolic,
1085  &Numeric, control, info);
1086  UMFPACK_DNAME (free_symbolic) (&Symbolic);
1087 
1088  rcond = Info (UMFPACK_RCOND);
1089 
1090  if (status < 0)
1091  {
1092  UMFPACK_DNAME (report_status) (control, status);
1093  UMFPACK_DNAME (report_info) (control, info);
1094 
1095  UMFPACK_DNAME (free_numeric) (&Numeric);
1096  (*current_liboctave_error_handler)
1097  ("SparseMatrix::determinant numeric factorization failed");
1098  }
1099  else
1100  {
1101  UMFPACK_DNAME (report_numeric) (Numeric, control);
1102 
1103  double c10, e10;
1104 
1105  status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1106  info);
1107 
1108  if (status < 0)
1109  {
1110  UMFPACK_DNAME (report_status) (control, status);
1111  UMFPACK_DNAME (report_info) (control, info);
1112 
1113  (*current_liboctave_error_handler)
1114  ("SparseMatrix::determinant error calculating determinant");
1115  }
1116  else
1117  retval = DET (c10, e10, 10);
1118 
1119  UMFPACK_DNAME (free_numeric) (&Numeric);
1120  }
1121  }
1122  }
1123 
1124 #else
1125 
1126  octave_unused_parameter (err);
1127  octave_unused_parameter (rcond);
1128 
1129  (*current_liboctave_error_handler)
1130  ("support for UMFPACK was unavailable or disabled "
1131  "when liboctave was built");
1132 
1133 #endif
1134 
1135  return retval;
1136 }
1137 
1138 Matrix
1141  double& rcond, solve_singularity_handler,
1142  bool calc_cond) const
1143 {
1144  Matrix retval;
1145 
1146  octave_idx_type nr = rows ();
1147  octave_idx_type nc = cols ();
1148  octave_idx_type nm = (nc < nr ? nc : nr);
1149  err = 0;
1150 
1151  if (nr != b.rows ())
1153  ("matrix dimension mismatch solution of linear equations");
1154 
1155  if (nr == 0 || nc == 0 || b.cols () == 0)
1156  retval = Matrix (nc, b.cols (), 0.0);
1157  else
1158  {
1159  // Print spparms("spumoni") info if requested
1160  int typ = mattype.type ();
1161  mattype.info ();
1162 
1164  (*current_liboctave_error_handler) ("incorrect matrix type");
1165 
1166  retval.resize (nc, b.cols (), 0.);
1167  if (typ == MatrixType::Diagonal)
1168  for (octave_idx_type j = 0; j < b.cols (); j++)
1169  for (octave_idx_type i = 0; i < nm; i++)
1170  retval(i,j) = b(i,j) / data (i);
1171  else
1172  for (octave_idx_type j = 0; j < b.cols (); j++)
1173  for (octave_idx_type k = 0; k < nc; k++)
1174  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1175  retval(k,j) = b(ridx (i),j) / data (i);
1176 
1177  if (calc_cond)
1178  {
1179  double dmax = 0.;
1180  double dmin = octave::numeric_limits<double>::Inf ();
1181  for (octave_idx_type i = 0; i < nm; i++)
1182  {
1183  double tmp = fabs (data (i));
1184  if (tmp > dmax)
1185  dmax = tmp;
1186  if (tmp < dmin)
1187  dmin = tmp;
1188  }
1189  rcond = dmin / dmax;
1190  }
1191  else
1192  rcond = 1.;
1193  }
1194 
1195  return retval;
1196 }
1197 
1200  octave_idx_type& err, double& rcond,
1201  solve_singularity_handler, bool calc_cond) const
1202 {
1204 
1205  octave_idx_type nr = rows ();
1206  octave_idx_type nc = cols ();
1207  octave_idx_type nm = (nc < nr ? nc : nr);
1208  err = 0;
1209 
1210  if (nr != b.rows ())
1212  ("matrix dimension mismatch solution of linear equations");
1213 
1214  if (nr == 0 || nc == 0 || b.cols () == 0)
1215  retval = SparseMatrix (nc, b.cols ());
1216  else
1217  {
1218  // Print spparms("spumoni") info if requested
1219  int typ = mattype.type ();
1220  mattype.info ();
1221 
1223  (*current_liboctave_error_handler) ("incorrect matrix type");
1224 
1225  octave_idx_type b_nc = b.cols ();
1226  octave_idx_type b_nz = b.nnz ();
1227  retval = SparseMatrix (nc, b_nc, b_nz);
1228 
1229  retval.xcidx (0) = 0;
1230  octave_idx_type ii = 0;
1231  if (typ == MatrixType::Diagonal)
1232  for (octave_idx_type j = 0; j < b_nc; j++)
1233  {
1234  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1235  {
1236  if (b.ridx (i) >= nm)
1237  break;
1238  retval.xridx (ii) = b.ridx (i);
1239  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1240  }
1241  retval.xcidx (j+1) = ii;
1242  }
1243  else
1244  for (octave_idx_type j = 0; j < b_nc; j++)
1245  {
1246  for (octave_idx_type l = 0; l < nc; l++)
1247  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1248  {
1249  bool found = false;
1251  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1252  if (ridx (i) == b.ridx (k))
1253  {
1254  found = true;
1255  break;
1256  }
1257  if (found)
1258  {
1259  retval.xridx (ii) = l;
1260  retval.xdata (ii++) = b.data (k) / data (i);
1261  }
1262  }
1263  retval.xcidx (j+1) = ii;
1264  }
1265 
1266  if (calc_cond)
1267  {
1268  double dmax = 0.;
1269  double dmin = octave::numeric_limits<double>::Inf ();
1270  for (octave_idx_type i = 0; i < nm; i++)
1271  {
1272  double tmp = fabs (data (i));
1273  if (tmp > dmax)
1274  dmax = tmp;
1275  if (tmp < dmin)
1276  dmin = tmp;
1277  }
1278  rcond = dmin / dmax;
1279  }
1280  else
1281  rcond = 1.;
1282  }
1283 
1284  return retval;
1285 }
1286 
1289  octave_idx_type& err, double& rcond,
1290  solve_singularity_handler, bool calc_cond) const
1291 {
1293 
1294  octave_idx_type nr = rows ();
1295  octave_idx_type nc = cols ();
1296  octave_idx_type nm = (nc < nr ? nc : nr);
1297  err = 0;
1298 
1299  if (nr != b.rows ())
1301  ("matrix dimension mismatch solution of linear equations");
1302 
1303  if (nr == 0 || nc == 0 || b.cols () == 0)
1304  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1305  else
1306  {
1307  // Print spparms("spumoni") info if requested
1308  int typ = mattype.type ();
1309  mattype.info ();
1310 
1312  (*current_liboctave_error_handler) ("incorrect matrix type");
1313 
1314  retval.resize (nc, b.cols (), 0);
1315  if (typ == MatrixType::Diagonal)
1316  for (octave_idx_type j = 0; j < b.cols (); j++)
1317  for (octave_idx_type i = 0; i < nm; i++)
1318  retval(i,j) = b(i,j) / data (i);
1319  else
1320  for (octave_idx_type j = 0; j < b.cols (); j++)
1321  for (octave_idx_type k = 0; k < nc; k++)
1322  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1323  retval(k,j) = b(ridx (i),j) / data (i);
1324 
1325  if (calc_cond)
1326  {
1327  double dmax = 0.;
1328  double dmin = octave::numeric_limits<double>::Inf ();
1329  for (octave_idx_type i = 0; i < nm; i++)
1330  {
1331  double tmp = fabs (data (i));
1332  if (tmp > dmax)
1333  dmax = tmp;
1334  if (tmp < dmin)
1335  dmin = tmp;
1336  }
1337  rcond = dmin / dmax;
1338  }
1339  else
1340  rcond = 1.;
1341  }
1342 
1343  return retval;
1344 }
1345 
1348  octave_idx_type& err, double& rcond,
1349  solve_singularity_handler, bool calc_cond) const
1350 {
1352 
1353  octave_idx_type nr = rows ();
1354  octave_idx_type nc = cols ();
1355  octave_idx_type nm = (nc < nr ? nc : nr);
1356  err = 0;
1357 
1358  if (nr != b.rows ())
1360  ("matrix dimension mismatch solution of linear equations");
1361 
1362  if (nr == 0 || nc == 0 || b.cols () == 0)
1363  retval = SparseComplexMatrix (nc, b.cols ());
1364  else
1365  {
1366  // Print spparms("spumoni") info if requested
1367  int typ = mattype.type ();
1368  mattype.info ();
1369 
1371  (*current_liboctave_error_handler) ("incorrect matrix type");
1372 
1373  octave_idx_type b_nc = b.cols ();
1374  octave_idx_type b_nz = b.nnz ();
1375  retval = SparseComplexMatrix (nc, b_nc, b_nz);
1376 
1377  retval.xcidx (0) = 0;
1378  octave_idx_type ii = 0;
1379  if (typ == MatrixType::Diagonal)
1380  for (octave_idx_type j = 0; j < b.cols (); j++)
1381  {
1382  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1383  {
1384  if (b.ridx (i) >= nm)
1385  break;
1386  retval.xridx (ii) = b.ridx (i);
1387  retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1388  }
1389  retval.xcidx (j+1) = ii;
1390  }
1391  else
1392  for (octave_idx_type j = 0; j < b.cols (); j++)
1393  {
1394  for (octave_idx_type l = 0; l < nc; l++)
1395  for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1396  {
1397  bool found = false;
1399  for (k = b.cidx (j); k < b.cidx (j+1); k++)
1400  if (ridx (i) == b.ridx (k))
1401  {
1402  found = true;
1403  break;
1404  }
1405  if (found)
1406  {
1407  retval.xridx (ii) = l;
1408  retval.xdata (ii++) = b.data (k) / data (i);
1409  }
1410  }
1411  retval.xcidx (j+1) = ii;
1412  }
1413 
1414  if (calc_cond)
1415  {
1416  double dmax = 0.;
1417  double dmin = octave::numeric_limits<double>::Inf ();
1418  for (octave_idx_type i = 0; i < nm; i++)
1419  {
1420  double tmp = fabs (data (i));
1421  if (tmp > dmax)
1422  dmax = tmp;
1423  if (tmp < dmin)
1424  dmin = tmp;
1425  }
1426  rcond = dmin / dmax;
1427  }
1428  else
1429  rcond = 1.;
1430  }
1431 
1432  return retval;
1433 }
1434 
1435 Matrix
1437  octave_idx_type& err, double& rcond,
1438  solve_singularity_handler sing_handler,
1439  bool calc_cond) const
1440 {
1441  Matrix 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 ())
1450  ("matrix dimension mismatch solution of linear equations");
1451 
1452  if (nr == 0 || nc == 0 || b.cols () == 0)
1453  retval = Matrix (nc, b.cols (), 0.0);
1454  else
1455  {
1456  // Print spparms("spumoni") info if requested
1457  int typ = mattype.type ();
1458  mattype.info ();
1459 
1460  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1461  (*current_liboctave_error_handler) ("incorrect matrix type");
1462 
1463  double anorm = 0.;
1464  double ainvnorm = 0.;
1465  octave_idx_type b_nc = b.cols ();
1466  rcond = 1.;
1467 
1468  if (calc_cond)
1469  {
1470  // Calculate the 1-norm of matrix for rcond calculation
1471  for (octave_idx_type j = 0; j < nc; j++)
1472  {
1473  double atmp = 0.;
1474  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1475  atmp += fabs (data (i));
1476  if (atmp > anorm)
1477  anorm = atmp;
1478  }
1479  }
1480 
1481  if (typ == MatrixType::Permuted_Upper)
1482  {
1483  retval.resize (nc, b_nc);
1484  octave_idx_type *perm = mattype.triangular_perm ();
1485  OCTAVE_LOCAL_BUFFER (double, work, nm);
1486 
1487  for (octave_idx_type j = 0; j < b_nc; j++)
1488  {
1489  for (octave_idx_type i = 0; i < nr; i++)
1490  work[i] = b(i,j);
1491  for (octave_idx_type i = nr; i < nc; i++)
1492  work[i] = 0.;
1493 
1494  for (octave_idx_type k = nc-1; k >= 0; k--)
1495  {
1496  octave_idx_type kidx = perm[k];
1497 
1498  if (work[k] != 0.)
1499  {
1500  if (ridx (cidx (kidx+1)-1) != k
1501  || data (cidx (kidx+1)-1) == 0.)
1502  {
1503  err = -2;
1504  goto triangular_error;
1505  }
1506 
1507  double tmp = work[k] / data (cidx (kidx+1)-1);
1508  work[k] = tmp;
1509  for (octave_idx_type i = cidx (kidx);
1510  i < cidx (kidx+1)-1; i++)
1511  {
1512  octave_idx_type iidx = ridx (i);
1513  work[iidx] = work[iidx] - tmp * data (i);
1514  }
1515  }
1516  }
1517 
1518  for (octave_idx_type i = 0; i < nc; i++)
1519  retval.xelem (perm[i], j) = work[i];
1520  }
1521 
1522  if (calc_cond)
1523  {
1524  // Calculation of 1-norm of inv(*this)
1525  for (octave_idx_type i = 0; i < nm; i++)
1526  work[i] = 0.;
1527 
1528  for (octave_idx_type j = 0; j < nr; j++)
1529  {
1530  work[j] = 1.;
1531 
1532  for (octave_idx_type k = j; k >= 0; k--)
1533  {
1534  octave_idx_type iidx = perm[k];
1535 
1536  if (work[k] != 0.)
1537  {
1538  double tmp = work[k] / data (cidx (iidx+1)-1);
1539  work[k] = tmp;
1540  for (octave_idx_type i = cidx (iidx);
1541  i < cidx (iidx+1)-1; i++)
1542  {
1543  octave_idx_type idx2 = ridx (i);
1544  work[idx2] = work[idx2] - tmp * data (i);
1545  }
1546  }
1547  }
1548  double atmp = 0;
1549  for (octave_idx_type i = 0; i < j+1; i++)
1550  {
1551  atmp += fabs (work[i]);
1552  work[i] = 0.;
1553  }
1554  if (atmp > ainvnorm)
1555  ainvnorm = atmp;
1556  }
1557  rcond = 1. / ainvnorm / anorm;
1558  }
1559  }
1560  else
1561  {
1562  OCTAVE_LOCAL_BUFFER (double, work, nm);
1563  retval.resize (nc, b_nc);
1564 
1565  for (octave_idx_type j = 0; j < b_nc; j++)
1566  {
1567  for (octave_idx_type i = 0; i < nr; i++)
1568  work[i] = b(i,j);
1569  for (octave_idx_type i = nr; i < nc; i++)
1570  work[i] = 0.;
1571 
1572  for (octave_idx_type k = nc-1; k >= 0; k--)
1573  {
1574  if (work[k] != 0.)
1575  {
1576  if (ridx (cidx (k+1)-1) != k
1577  || data (cidx (k+1)-1) == 0.)
1578  {
1579  err = -2;
1580  goto triangular_error;
1581  }
1582 
1583  double tmp = work[k] / data (cidx (k+1)-1);
1584  work[k] = tmp;
1585  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1586  {
1587  octave_idx_type iidx = ridx (i);
1588  work[iidx] = work[iidx] - tmp * data (i);
1589  }
1590  }
1591  }
1592 
1593  for (octave_idx_type i = 0; i < nc; i++)
1594  retval.xelem (i, j) = work[i];
1595  }
1596 
1597  if (calc_cond)
1598  {
1599  // Calculation of 1-norm of inv(*this)
1600  for (octave_idx_type i = 0; i < nm; i++)
1601  work[i] = 0.;
1602 
1603  for (octave_idx_type j = 0; j < nr; j++)
1604  {
1605  work[j] = 1.;
1606 
1607  for (octave_idx_type k = j; k >= 0; k--)
1608  {
1609  if (work[k] != 0.)
1610  {
1611  double tmp = work[k] / data (cidx (k+1)-1);
1612  work[k] = tmp;
1613  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1614  {
1615  octave_idx_type iidx = ridx (i);
1616  work[iidx] = work[iidx] - tmp * data (i);
1617  }
1618  }
1619  }
1620  double atmp = 0;
1621  for (octave_idx_type i = 0; i < j+1; i++)
1622  {
1623  atmp += fabs (work[i]);
1624  work[i] = 0.;
1625  }
1626  if (atmp > ainvnorm)
1627  ainvnorm = atmp;
1628  }
1629  rcond = 1. / ainvnorm / anorm;
1630  }
1631  }
1632 
1633  triangular_error:
1634  if (err != 0)
1635  {
1636  if (sing_handler)
1637  {
1638  sing_handler (rcond);
1639  mattype.mark_as_rectangular ();
1640  }
1641  else
1643  }
1644 
1645  volatile double rcond_plus_one = rcond + 1.0;
1646 
1647  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
1648  {
1649  err = -2;
1650 
1651  if (sing_handler)
1652  {
1653  sing_handler (rcond);
1654  mattype.mark_as_rectangular ();
1655  }
1656  else
1658  }
1659  }
1660 
1661  return retval;
1662 }
1663 
1666  octave_idx_type& err, double& rcond,
1667  solve_singularity_handler sing_handler,
1668  bool calc_cond) const
1669 {
1671 
1672  octave_idx_type nr = rows ();
1673  octave_idx_type nc = cols ();
1674  octave_idx_type nm = (nc > nr ? nc : nr);
1675  err = 0;
1676 
1677  if (nr != b.rows ())
1679  ("matrix dimension mismatch solution of linear equations");
1680 
1681  if (nr == 0 || nc == 0 || b.cols () == 0)
1682  retval = SparseMatrix (nc, b.cols ());
1683  else
1684  {
1685  // Print spparms("spumoni") info if requested
1686  int typ = mattype.type ();
1687  mattype.info ();
1688 
1689  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1690  (*current_liboctave_error_handler) ("incorrect matrix type");
1691 
1692  double anorm = 0.;
1693  double ainvnorm = 0.;
1694  rcond = 1.;
1695 
1696  if (calc_cond)
1697  {
1698  // Calculate the 1-norm of matrix for rcond calculation
1699  for (octave_idx_type j = 0; j < nc; j++)
1700  {
1701  double atmp = 0.;
1702  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1703  atmp += fabs (data (i));
1704  if (atmp > anorm)
1705  anorm = atmp;
1706  }
1707  }
1708 
1709  octave_idx_type b_nc = b.cols ();
1710  octave_idx_type b_nz = b.nnz ();
1711  retval = SparseMatrix (nc, b_nc, b_nz);
1712  retval.xcidx (0) = 0;
1713  octave_idx_type ii = 0;
1714  octave_idx_type x_nz = b_nz;
1715 
1716  if (typ == MatrixType::Permuted_Upper)
1717  {
1718  octave_idx_type *perm = mattype.triangular_perm ();
1719  OCTAVE_LOCAL_BUFFER (double, work, nm);
1720 
1721  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1722  for (octave_idx_type i = 0; i < nc; i++)
1723  rperm[perm[i]] = i;
1724 
1725  for (octave_idx_type j = 0; j < b_nc; j++)
1726  {
1727  for (octave_idx_type i = 0; i < nm; i++)
1728  work[i] = 0.;
1729  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1730  work[b.ridx (i)] = b.data (i);
1731 
1732  for (octave_idx_type k = nc-1; k >= 0; k--)
1733  {
1734  octave_idx_type kidx = perm[k];
1735 
1736  if (work[k] != 0.)
1737  {
1738  if (ridx (cidx (kidx+1)-1) != k
1739  || data (cidx (kidx+1)-1) == 0.)
1740  {
1741  err = -2;
1742  goto triangular_error;
1743  }
1744 
1745  double tmp = work[k] / data (cidx (kidx+1)-1);
1746  work[k] = tmp;
1747  for (octave_idx_type i = cidx (kidx);
1748  i < cidx (kidx+1)-1; i++)
1749  {
1750  octave_idx_type iidx = ridx (i);
1751  work[iidx] = work[iidx] - tmp * data (i);
1752  }
1753  }
1754  }
1755 
1756  // Count nonzeros in work vector and adjust space in
1757  // retval if needed
1758  octave_idx_type new_nnz = 0;
1759  for (octave_idx_type i = 0; i < nc; i++)
1760  if (work[i] != 0.)
1761  new_nnz++;
1762 
1763  if (ii + new_nnz > x_nz)
1764  {
1765  // Resize the sparse matrix
1766  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1767  retval.change_capacity (sz);
1768  x_nz = sz;
1769  }
1770 
1771  for (octave_idx_type i = 0; i < nc; i++)
1772  if (work[rperm[i]] != 0.)
1773  {
1774  retval.xridx (ii) = i;
1775  retval.xdata (ii++) = work[rperm[i]];
1776  }
1777  retval.xcidx (j+1) = ii;
1778  }
1779 
1780  retval.maybe_compress ();
1781 
1782  if (calc_cond)
1783  {
1784  // Calculation of 1-norm of inv(*this)
1785  for (octave_idx_type i = 0; i < nm; i++)
1786  work[i] = 0.;
1787 
1788  for (octave_idx_type j = 0; j < nr; j++)
1789  {
1790  work[j] = 1.;
1791 
1792  for (octave_idx_type k = j; k >= 0; k--)
1793  {
1794  octave_idx_type iidx = perm[k];
1795 
1796  if (work[k] != 0.)
1797  {
1798  double tmp = work[k] / data (cidx (iidx+1)-1);
1799  work[k] = tmp;
1800  for (octave_idx_type i = cidx (iidx);
1801  i < cidx (iidx+1)-1; i++)
1802  {
1803  octave_idx_type idx2 = ridx (i);
1804  work[idx2] = work[idx2] - tmp * data (i);
1805  }
1806  }
1807  }
1808  double atmp = 0;
1809  for (octave_idx_type i = 0; i < j+1; i++)
1810  {
1811  atmp += fabs (work[i]);
1812  work[i] = 0.;
1813  }
1814  if (atmp > ainvnorm)
1815  ainvnorm = atmp;
1816  }
1817  rcond = 1. / ainvnorm / anorm;
1818  }
1819  }
1820  else
1821  {
1822  OCTAVE_LOCAL_BUFFER (double, work, nm);
1823 
1824  for (octave_idx_type j = 0; j < b_nc; j++)
1825  {
1826  for (octave_idx_type i = 0; i < nm; i++)
1827  work[i] = 0.;
1828  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1829  work[b.ridx (i)] = b.data (i);
1830 
1831  for (octave_idx_type k = nc-1; k >= 0; k--)
1832  {
1833  if (work[k] != 0.)
1834  {
1835  if (ridx (cidx (k+1)-1) != k
1836  || data (cidx (k+1)-1) == 0.)
1837  {
1838  err = -2;
1839  goto triangular_error;
1840  }
1841 
1842  double tmp = work[k] / data (cidx (k+1)-1);
1843  work[k] = tmp;
1844  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1845  {
1846  octave_idx_type iidx = ridx (i);
1847  work[iidx] = work[iidx] - tmp * data (i);
1848  }
1849  }
1850  }
1851 
1852  // Count nonzeros in work vector and adjust space in
1853  // retval if needed
1854  octave_idx_type new_nnz = 0;
1855  for (octave_idx_type i = 0; i < nc; i++)
1856  if (work[i] != 0.)
1857  new_nnz++;
1858 
1859  if (ii + new_nnz > x_nz)
1860  {
1861  // Resize the sparse matrix
1862  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1863  retval.change_capacity (sz);
1864  x_nz = sz;
1865  }
1866 
1867  for (octave_idx_type i = 0; i < nc; i++)
1868  if (work[i] != 0.)
1869  {
1870  retval.xridx (ii) = i;
1871  retval.xdata (ii++) = work[i];
1872  }
1873  retval.xcidx (j+1) = ii;
1874  }
1875 
1876  retval.maybe_compress ();
1877 
1878  if (calc_cond)
1879  {
1880  // Calculation of 1-norm of inv(*this)
1881  for (octave_idx_type i = 0; i < nm; i++)
1882  work[i] = 0.;
1883 
1884  for (octave_idx_type j = 0; j < nr; j++)
1885  {
1886  work[j] = 1.;
1887 
1888  for (octave_idx_type k = j; k >= 0; k--)
1889  {
1890  if (work[k] != 0.)
1891  {
1892  double tmp = work[k] / data (cidx (k+1)-1);
1893  work[k] = tmp;
1894  for (octave_idx_type i = cidx (k);
1895  i < cidx (k+1)-1; i++)
1896  {
1897  octave_idx_type iidx = ridx (i);
1898  work[iidx] = work[iidx] - tmp * data (i);
1899  }
1900  }
1901  }
1902  double atmp = 0;
1903  for (octave_idx_type i = 0; i < j+1; i++)
1904  {
1905  atmp += fabs (work[i]);
1906  work[i] = 0.;
1907  }
1908  if (atmp > ainvnorm)
1909  ainvnorm = atmp;
1910  }
1911  rcond = 1. / ainvnorm / anorm;
1912  }
1913  }
1914 
1915  triangular_error:
1916  if (err != 0)
1917  {
1918  if (sing_handler)
1919  {
1920  sing_handler (rcond);
1921  mattype.mark_as_rectangular ();
1922  }
1923  else
1925  }
1926 
1927  volatile double rcond_plus_one = rcond + 1.0;
1928 
1929  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
1930  {
1931  err = -2;
1932 
1933  if (sing_handler)
1934  {
1935  sing_handler (rcond);
1936  mattype.mark_as_rectangular ();
1937  }
1938  else
1940  }
1941  }
1942  return retval;
1943 }
1944 
1947  octave_idx_type& err, double& rcond,
1948  solve_singularity_handler sing_handler,
1949  bool calc_cond) const
1950 {
1952 
1953  octave_idx_type nr = rows ();
1954  octave_idx_type nc = cols ();
1955  octave_idx_type nm = (nc > nr ? nc : nr);
1956  err = 0;
1957 
1958  if (nr != b.rows ())
1960  ("matrix dimension mismatch solution of linear equations");
1961 
1962  if (nr == 0 || nc == 0 || b.cols () == 0)
1963  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1964  else
1965  {
1966  // Print spparms("spumoni") info if requested
1967  int typ = mattype.type ();
1968  mattype.info ();
1969 
1970  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1971  (*current_liboctave_error_handler) ("incorrect matrix type");
1972 
1973  double anorm = 0.;
1974  double ainvnorm = 0.;
1975  octave_idx_type b_nc = b.cols ();
1976  rcond = 1.;
1977 
1978  if (calc_cond)
1979  {
1980  // Calculate the 1-norm of matrix for rcond calculation
1981  for (octave_idx_type j = 0; j < nc; j++)
1982  {
1983  double atmp = 0.;
1984  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1985  atmp += fabs (data (i));
1986  if (atmp > anorm)
1987  anorm = atmp;
1988  }
1989  }
1990 
1991  if (typ == MatrixType::Permuted_Upper)
1992  {
1993  retval.resize (nc, b_nc);
1994  octave_idx_type *perm = mattype.triangular_perm ();
1995  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
1996 
1997  for (octave_idx_type j = 0; j < b_nc; j++)
1998  {
1999  for (octave_idx_type i = 0; i < nr; i++)
2000  cwork[i] = b(i,j);
2001  for (octave_idx_type i = nr; i < nc; i++)
2002  cwork[i] = 0.;
2003 
2004  for (octave_idx_type k = nc-1; k >= 0; k--)
2005  {
2006  octave_idx_type kidx = perm[k];
2007 
2008  if (cwork[k] != 0.)
2009  {
2010  if (ridx (cidx (kidx+1)-1) != k
2011  || data (cidx (kidx+1)-1) == 0.)
2012  {
2013  err = -2;
2014  goto triangular_error;
2015  }
2016 
2017  Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2018  cwork[k] = tmp;
2019  for (octave_idx_type i = cidx (kidx);
2020  i < cidx (kidx+1)-1; i++)
2021  {
2022  octave_idx_type iidx = ridx (i);
2023  cwork[iidx] = cwork[iidx] - tmp * data (i);
2024  }
2025  }
2026  }
2027 
2028  for (octave_idx_type i = 0; i < nc; i++)
2029  retval.xelem (perm[i], j) = cwork[i];
2030  }
2031 
2032  if (calc_cond)
2033  {
2034  // Calculation of 1-norm of inv(*this)
2035  OCTAVE_LOCAL_BUFFER (double, work, nm);
2036  for (octave_idx_type i = 0; i < nm; i++)
2037  work[i] = 0.;
2038 
2039  for (octave_idx_type j = 0; j < nr; j++)
2040  {
2041  work[j] = 1.;
2042 
2043  for (octave_idx_type k = j; k >= 0; k--)
2044  {
2045  octave_idx_type iidx = perm[k];
2046 
2047  if (work[k] != 0.)
2048  {
2049  double tmp = work[k] / data (cidx (iidx+1)-1);
2050  work[k] = tmp;
2051  for (octave_idx_type i = cidx (iidx);
2052  i < cidx (iidx+1)-1; i++)
2053  {
2054  octave_idx_type idx2 = ridx (i);
2055  work[idx2] = work[idx2] - tmp * data (i);
2056  }
2057  }
2058  }
2059  double atmp = 0;
2060  for (octave_idx_type i = 0; i < j+1; i++)
2061  {
2062  atmp += fabs (work[i]);
2063  work[i] = 0.;
2064  }
2065  if (atmp > ainvnorm)
2066  ainvnorm = atmp;
2067  }
2068  rcond = 1. / ainvnorm / anorm;
2069  }
2070  }
2071  else
2072  {
2073  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2074  retval.resize (nc, b_nc);
2075 
2076  for (octave_idx_type j = 0; j < b_nc; j++)
2077  {
2078  for (octave_idx_type i = 0; i < nr; i++)
2079  cwork[i] = b(i,j);
2080  for (octave_idx_type i = nr; i < nc; i++)
2081  cwork[i] = 0.;
2082 
2083  for (octave_idx_type k = nc-1; k >= 0; k--)
2084  {
2085  if (cwork[k] != 0.)
2086  {
2087  if (ridx (cidx (k+1)-1) != k
2088  || data (cidx (k+1)-1) == 0.)
2089  {
2090  err = -2;
2091  goto triangular_error;
2092  }
2093 
2094  Complex tmp = cwork[k] / data (cidx (k+1)-1);
2095  cwork[k] = tmp;
2096  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2097  {
2098  octave_idx_type iidx = ridx (i);
2099  cwork[iidx] = cwork[iidx] - tmp * data (i);
2100  }
2101  }
2102  }
2103 
2104  for (octave_idx_type i = 0; i < nc; i++)
2105  retval.xelem (i, j) = cwork[i];
2106  }
2107 
2108  if (calc_cond)
2109  {
2110  // Calculation of 1-norm of inv(*this)
2111  OCTAVE_LOCAL_BUFFER (double, work, nm);
2112  for (octave_idx_type i = 0; i < nm; i++)
2113  work[i] = 0.;
2114 
2115  for (octave_idx_type j = 0; j < nr; j++)
2116  {
2117  work[j] = 1.;
2118 
2119  for (octave_idx_type k = j; k >= 0; k--)
2120  {
2121  if (work[k] != 0.)
2122  {
2123  double tmp = work[k] / data (cidx (k+1)-1);
2124  work[k] = tmp;
2125  for (octave_idx_type i = cidx (k);
2126  i < cidx (k+1)-1; i++)
2127  {
2128  octave_idx_type iidx = ridx (i);
2129  work[iidx] = work[iidx] - tmp * data (i);
2130  }
2131  }
2132  }
2133  double atmp = 0;
2134  for (octave_idx_type i = 0; i < j+1; i++)
2135  {
2136  atmp += fabs (work[i]);
2137  work[i] = 0.;
2138  }
2139  if (atmp > ainvnorm)
2140  ainvnorm = atmp;
2141  }
2142  rcond = 1. / ainvnorm / anorm;
2143  }
2144  }
2145 
2146  triangular_error:
2147  if (err != 0)
2148  {
2149  if (sing_handler)
2150  {
2151  sing_handler (rcond);
2152  mattype.mark_as_rectangular ();
2153  }
2154  else
2156  }
2157 
2158  volatile double rcond_plus_one = rcond + 1.0;
2159 
2160  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2161  {
2162  err = -2;
2163 
2164  if (sing_handler)
2165  {
2166  sing_handler (rcond);
2167  mattype.mark_as_rectangular ();
2168  }
2169  else
2171  }
2172  }
2173 
2174  return retval;
2175 }
2176 
2179  octave_idx_type& err, double& rcond,
2180  solve_singularity_handler sing_handler,
2181  bool calc_cond) const
2182 {
2184 
2185  octave_idx_type nr = rows ();
2186  octave_idx_type nc = cols ();
2187  octave_idx_type nm = (nc > nr ? nc : nr);
2188  err = 0;
2189 
2190  if (nr != b.rows ())
2192  ("matrix dimension mismatch solution of linear equations");
2193 
2194  if (nr == 0 || nc == 0 || b.cols () == 0)
2195  retval = SparseComplexMatrix (nc, b.cols ());
2196  else
2197  {
2198  // Print spparms("spumoni") info if requested
2199  int typ = mattype.type ();
2200  mattype.info ();
2201 
2202  if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2203  (*current_liboctave_error_handler) ("incorrect matrix type");
2204 
2205  double anorm = 0.;
2206  double ainvnorm = 0.;
2207  rcond = 1.;
2208 
2209  if (calc_cond)
2210  {
2211  // Calculate the 1-norm of matrix for rcond calculation
2212  for (octave_idx_type j = 0; j < nc; j++)
2213  {
2214  double atmp = 0.;
2215  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2216  atmp += fabs (data (i));
2217  if (atmp > anorm)
2218  anorm = atmp;
2219  }
2220  }
2221 
2222  octave_idx_type b_nc = b.cols ();
2223  octave_idx_type b_nz = b.nnz ();
2224  retval = SparseComplexMatrix (nc, b_nc, b_nz);
2225  retval.xcidx (0) = 0;
2226  octave_idx_type ii = 0;
2227  octave_idx_type x_nz = b_nz;
2228 
2229  if (typ == MatrixType::Permuted_Upper)
2230  {
2231  octave_idx_type *perm = mattype.triangular_perm ();
2232  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2233 
2234  OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
2235  for (octave_idx_type i = 0; i < nc; i++)
2236  rperm[perm[i]] = i;
2237 
2238  for (octave_idx_type j = 0; j < b_nc; j++)
2239  {
2240  for (octave_idx_type i = 0; i < nm; i++)
2241  cwork[i] = 0.;
2242  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2243  cwork[b.ridx (i)] = b.data (i);
2244 
2245  for (octave_idx_type k = nc-1; k >= 0; k--)
2246  {
2247  octave_idx_type kidx = perm[k];
2248 
2249  if (cwork[k] != 0.)
2250  {
2251  if (ridx (cidx (kidx+1)-1) != k
2252  || data (cidx (kidx+1)-1) == 0.)
2253  {
2254  err = -2;
2255  goto triangular_error;
2256  }
2257 
2258  Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2259  cwork[k] = tmp;
2260  for (octave_idx_type i = cidx (kidx);
2261  i < cidx (kidx+1)-1; i++)
2262  {
2263  octave_idx_type iidx = ridx (i);
2264  cwork[iidx] = cwork[iidx] - tmp * data (i);
2265  }
2266  }
2267  }
2268 
2269  // Count nonzeros in work vector and adjust space in
2270  // retval if needed
2271  octave_idx_type new_nnz = 0;
2272  for (octave_idx_type i = 0; i < nc; i++)
2273  if (cwork[i] != 0.)
2274  new_nnz++;
2275 
2276  if (ii + new_nnz > x_nz)
2277  {
2278  // Resize the sparse matrix
2279  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2280  retval.change_capacity (sz);
2281  x_nz = sz;
2282  }
2283 
2284  for (octave_idx_type i = 0; i < nc; i++)
2285  if (cwork[rperm[i]] != 0.)
2286  {
2287  retval.xridx (ii) = i;
2288  retval.xdata (ii++) = cwork[rperm[i]];
2289  }
2290  retval.xcidx (j+1) = ii;
2291  }
2292 
2293  retval.maybe_compress ();
2294 
2295  if (calc_cond)
2296  {
2297  // Calculation of 1-norm of inv(*this)
2298  OCTAVE_LOCAL_BUFFER (double, work, nm);
2299  for (octave_idx_type i = 0; i < nm; i++)
2300  work[i] = 0.;
2301 
2302  for (octave_idx_type j = 0; j < nr; j++)
2303  {
2304  work[j] = 1.;
2305 
2306  for (octave_idx_type k = j; k >= 0; k--)
2307  {
2308  octave_idx_type iidx = perm[k];
2309 
2310  if (work[k] != 0.)
2311  {
2312  double tmp = work[k] / data (cidx (iidx+1)-1);
2313  work[k] = tmp;
2314  for (octave_idx_type i = cidx (iidx);
2315  i < cidx (iidx+1)-1; i++)
2316  {
2317  octave_idx_type idx2 = ridx (i);
2318  work[idx2] = work[idx2] - tmp * data (i);
2319  }
2320  }
2321  }
2322  double atmp = 0;
2323  for (octave_idx_type i = 0; i < j+1; i++)
2324  {
2325  atmp += fabs (work[i]);
2326  work[i] = 0.;
2327  }
2328  if (atmp > ainvnorm)
2329  ainvnorm = atmp;
2330  }
2331  rcond = 1. / ainvnorm / anorm;
2332  }
2333  }
2334  else
2335  {
2336  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2337 
2338  for (octave_idx_type j = 0; j < b_nc; j++)
2339  {
2340  for (octave_idx_type i = 0; i < nm; i++)
2341  cwork[i] = 0.;
2342  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2343  cwork[b.ridx (i)] = b.data (i);
2344 
2345  for (octave_idx_type k = nc-1; k >= 0; k--)
2346  {
2347  if (cwork[k] != 0.)
2348  {
2349  if (ridx (cidx (k+1)-1) != k
2350  || data (cidx (k+1)-1) == 0.)
2351  {
2352  err = -2;
2353  goto triangular_error;
2354  }
2355 
2356  Complex tmp = cwork[k] / data (cidx (k+1)-1);
2357  cwork[k] = tmp;
2358  for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2359  {
2360  octave_idx_type iidx = ridx (i);
2361  cwork[iidx] = cwork[iidx] - tmp * data (i);
2362  }
2363  }
2364  }
2365 
2366  // Count nonzeros in work vector and adjust space in
2367  // retval if needed
2368  octave_idx_type new_nnz = 0;
2369  for (octave_idx_type i = 0; i < nc; i++)
2370  if (cwork[i] != 0.)
2371  new_nnz++;
2372 
2373  if (ii + new_nnz > x_nz)
2374  {
2375  // Resize the sparse matrix
2376  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2377  retval.change_capacity (sz);
2378  x_nz = sz;
2379  }
2380 
2381  for (octave_idx_type i = 0; i < nc; i++)
2382  if (cwork[i] != 0.)
2383  {
2384  retval.xridx (ii) = i;
2385  retval.xdata (ii++) = cwork[i];
2386  }
2387  retval.xcidx (j+1) = ii;
2388  }
2389 
2390  retval.maybe_compress ();
2391 
2392  if (calc_cond)
2393  {
2394  // Calculation of 1-norm of inv(*this)
2395  OCTAVE_LOCAL_BUFFER (double, work, nm);
2396  for (octave_idx_type i = 0; i < nm; i++)
2397  work[i] = 0.;
2398 
2399  for (octave_idx_type j = 0; j < nr; j++)
2400  {
2401  work[j] = 1.;
2402 
2403  for (octave_idx_type k = j; k >= 0; k--)
2404  {
2405  if (work[k] != 0.)
2406  {
2407  double tmp = work[k] / data (cidx (k+1)-1);
2408  work[k] = tmp;
2409  for (octave_idx_type i = cidx (k);
2410  i < cidx (k+1)-1; i++)
2411  {
2412  octave_idx_type iidx = ridx (i);
2413  work[iidx] = work[iidx] - tmp * data (i);
2414  }
2415  }
2416  }
2417  double atmp = 0;
2418  for (octave_idx_type i = 0; i < j+1; i++)
2419  {
2420  atmp += fabs (work[i]);
2421  work[i] = 0.;
2422  }
2423  if (atmp > ainvnorm)
2424  ainvnorm = atmp;
2425  }
2426  rcond = 1. / ainvnorm / anorm;
2427  }
2428  }
2429 
2430  triangular_error:
2431  if (err != 0)
2432  {
2433  if (sing_handler)
2434  {
2435  sing_handler (rcond);
2436  mattype.mark_as_rectangular ();
2437  }
2438  else
2440  }
2441 
2442  volatile double rcond_plus_one = rcond + 1.0;
2443 
2444  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2445  {
2446  err = -2;
2447 
2448  if (sing_handler)
2449  {
2450  sing_handler (rcond);
2451  mattype.mark_as_rectangular ();
2452  }
2453  else
2455  }
2456  }
2457 
2458  return retval;
2459 }
2460 
2461 Matrix
2463  octave_idx_type& err, double& rcond,
2464  solve_singularity_handler sing_handler,
2465  bool calc_cond) const
2466 {
2467  Matrix retval;
2468 
2469  octave_idx_type nr = rows ();
2470  octave_idx_type nc = cols ();
2471  octave_idx_type nm = (nc > nr ? nc : nr);
2472  err = 0;
2473 
2474  if (nr != b.rows ())
2476  ("matrix dimension mismatch solution of linear equations");
2477 
2478  if (nr == 0 || nc == 0 || b.cols () == 0)
2479  retval = Matrix (nc, b.cols (), 0.0);
2480  else
2481  {
2482  // Print spparms("spumoni") info if requested
2483  int typ = mattype.type ();
2484  mattype.info ();
2485 
2486  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2487  (*current_liboctave_error_handler) ("incorrect matrix type");
2488 
2489  double anorm = 0.;
2490  double ainvnorm = 0.;
2491  octave_idx_type b_nc = b.cols ();
2492  rcond = 1.;
2493 
2494  if (calc_cond)
2495  {
2496  // Calculate the 1-norm of matrix for rcond calculation
2497  for (octave_idx_type j = 0; j < nc; j++)
2498  {
2499  double atmp = 0.;
2500  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2501  atmp += fabs (data (i));
2502  if (atmp > anorm)
2503  anorm = atmp;
2504  }
2505  }
2506 
2507  if (typ == MatrixType::Permuted_Lower)
2508  {
2509  retval.resize (nc, b_nc);
2510  OCTAVE_LOCAL_BUFFER (double, work, nm);
2511  octave_idx_type *perm = mattype.triangular_perm ();
2512 
2513  for (octave_idx_type j = 0; j < b_nc; j++)
2514  {
2515  if (nc > nr)
2516  for (octave_idx_type i = 0; i < nm; i++)
2517  work[i] = 0.;
2518  for (octave_idx_type i = 0; i < nr; i++)
2519  work[perm[i]] = b(i,j);
2520 
2521  for (octave_idx_type k = 0; k < nc; k++)
2522  {
2523  if (work[k] != 0.)
2524  {
2525  octave_idx_type minr = nr;
2526  octave_idx_type mini = 0;
2527 
2528  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2529  if (perm[ridx (i)] < minr)
2530  {
2531  minr = perm[ridx (i)];
2532  mini = i;
2533  }
2534 
2535  if (minr != k || data (mini) == 0)
2536  {
2537  err = -2;
2538  goto triangular_error;
2539  }
2540 
2541  double tmp = work[k] / data (mini);
2542  work[k] = tmp;
2543  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2544  {
2545  if (i == mini)
2546  continue;
2547 
2548  octave_idx_type iidx = perm[ridx (i)];
2549  work[iidx] = work[iidx] - tmp * data (i);
2550  }
2551  }
2552  }
2553 
2554  for (octave_idx_type i = 0; i < nc; i++)
2555  retval(i, j) = work[i];
2556  }
2557 
2558  if (calc_cond)
2559  {
2560  // Calculation of 1-norm of inv(*this)
2561  for (octave_idx_type i = 0; i < nm; i++)
2562  work[i] = 0.;
2563 
2564  for (octave_idx_type j = 0; j < nr; j++)
2565  {
2566  work[j] = 1.;
2567 
2568  for (octave_idx_type k = 0; k < nc; k++)
2569  {
2570  if (work[k] != 0.)
2571  {
2572  octave_idx_type minr = nr;
2573  octave_idx_type mini = 0;
2574 
2575  for (octave_idx_type i = cidx (k);
2576  i < cidx (k+1); i++)
2577  if (perm[ridx (i)] < minr)
2578  {
2579  minr = perm[ridx (i)];
2580  mini = i;
2581  }
2582 
2583  double tmp = work[k] / data (mini);
2584  work[k] = tmp;
2585  for (octave_idx_type i = cidx (k);
2586  i < cidx (k+1); i++)
2587  {
2588  if (i == mini)
2589  continue;
2590 
2591  octave_idx_type iidx = perm[ridx (i)];
2592  work[iidx] = work[iidx] - tmp * data (i);
2593  }
2594  }
2595  }
2596 
2597  double atmp = 0;
2598  for (octave_idx_type i = j; i < nc; i++)
2599  {
2600  atmp += fabs (work[i]);
2601  work[i] = 0.;
2602  }
2603  if (atmp > ainvnorm)
2604  ainvnorm = atmp;
2605  }
2606  rcond = 1. / ainvnorm / anorm;
2607  }
2608  }
2609  else
2610  {
2611  OCTAVE_LOCAL_BUFFER (double, work, nm);
2612  retval.resize (nc, b_nc, 0.);
2613 
2614  for (octave_idx_type j = 0; j < b_nc; j++)
2615  {
2616  for (octave_idx_type i = 0; i < nr; i++)
2617  work[i] = b(i,j);
2618  for (octave_idx_type i = nr; i < nc; i++)
2619  work[i] = 0.;
2620  for (octave_idx_type k = 0; k < nc; k++)
2621  {
2622  if (work[k] != 0.)
2623  {
2624  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2625  {
2626  err = -2;
2627  goto triangular_error;
2628  }
2629 
2630  double tmp = work[k] / data (cidx (k));
2631  work[k] = tmp;
2632  for (octave_idx_type i = cidx (k)+1;
2633  i < cidx (k+1); i++)
2634  {
2635  octave_idx_type iidx = ridx (i);
2636  work[iidx] = work[iidx] - tmp * data (i);
2637  }
2638  }
2639  }
2640 
2641  for (octave_idx_type i = 0; i < nc; i++)
2642  retval.xelem (i, j) = work[i];
2643  }
2644 
2645  if (calc_cond)
2646  {
2647  // Calculation of 1-norm of inv(*this)
2648  for (octave_idx_type i = 0; i < nm; i++)
2649  work[i] = 0.;
2650 
2651  for (octave_idx_type j = 0; j < nr; j++)
2652  {
2653  work[j] = 1.;
2654 
2655  for (octave_idx_type k = j; k < nc; k++)
2656  {
2657 
2658  if (work[k] != 0.)
2659  {
2660  double tmp = work[k] / data (cidx (k));
2661  work[k] = tmp;
2662  for (octave_idx_type i = cidx (k)+1;
2663  i < cidx (k+1); i++)
2664  {
2665  octave_idx_type iidx = ridx (i);
2666  work[iidx] = work[iidx] - tmp * data (i);
2667  }
2668  }
2669  }
2670  double atmp = 0;
2671  for (octave_idx_type i = j; i < nc; i++)
2672  {
2673  atmp += fabs (work[i]);
2674  work[i] = 0.;
2675  }
2676  if (atmp > ainvnorm)
2677  ainvnorm = atmp;
2678  }
2679  rcond = 1. / ainvnorm / anorm;
2680  }
2681  }
2682 
2683  triangular_error:
2684  if (err != 0)
2685  {
2686  if (sing_handler)
2687  {
2688  sing_handler (rcond);
2689  mattype.mark_as_rectangular ();
2690  }
2691  else
2693  }
2694 
2695  volatile double rcond_plus_one = rcond + 1.0;
2696 
2697  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2698  {
2699  err = -2;
2700 
2701  if (sing_handler)
2702  {
2703  sing_handler (rcond);
2704  mattype.mark_as_rectangular ();
2705  }
2706  else
2708  }
2709  }
2710 
2711  return retval;
2712 }
2713 
2716  octave_idx_type& err, double& rcond,
2717  solve_singularity_handler sing_handler,
2718  bool calc_cond) const
2719 {
2721 
2722  octave_idx_type nr = rows ();
2723  octave_idx_type nc = cols ();
2724  octave_idx_type nm = (nc > nr ? nc : nr);
2725  err = 0;
2726 
2727  if (nr != b.rows ())
2729  ("matrix dimension mismatch solution of linear equations");
2730 
2731  if (nr == 0 || nc == 0 || b.cols () == 0)
2732  retval = SparseMatrix (nc, b.cols ());
2733  else
2734  {
2735  // Print spparms("spumoni") info if requested
2736  int typ = mattype.type ();
2737  mattype.info ();
2738 
2739  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2740  (*current_liboctave_error_handler) ("incorrect matrix type");
2741 
2742  double anorm = 0.;
2743  double ainvnorm = 0.;
2744  rcond = 1.;
2745 
2746  if (calc_cond)
2747  {
2748  // Calculate the 1-norm of matrix for rcond calculation
2749  for (octave_idx_type j = 0; j < nc; j++)
2750  {
2751  double atmp = 0.;
2752  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2753  atmp += fabs (data (i));
2754  if (atmp > anorm)
2755  anorm = atmp;
2756  }
2757  }
2758 
2759  octave_idx_type b_nc = b.cols ();
2760  octave_idx_type b_nz = b.nnz ();
2761  retval = SparseMatrix (nc, b_nc, b_nz);
2762  retval.xcidx (0) = 0;
2763  octave_idx_type ii = 0;
2764  octave_idx_type x_nz = b_nz;
2765 
2766  if (typ == MatrixType::Permuted_Lower)
2767  {
2768  OCTAVE_LOCAL_BUFFER (double, work, nm);
2769  octave_idx_type *perm = mattype.triangular_perm ();
2770 
2771  for (octave_idx_type j = 0; j < b_nc; j++)
2772  {
2773  for (octave_idx_type i = 0; i < nm; i++)
2774  work[i] = 0.;
2775  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2776  work[perm[b.ridx (i)]] = b.data (i);
2777 
2778  for (octave_idx_type k = 0; k < nc; k++)
2779  {
2780  if (work[k] != 0.)
2781  {
2782  octave_idx_type minr = nr;
2783  octave_idx_type mini = 0;
2784 
2785  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2786  if (perm[ridx (i)] < minr)
2787  {
2788  minr = perm[ridx (i)];
2789  mini = i;
2790  }
2791 
2792  if (minr != k || data (mini) == 0)
2793  {
2794  err = -2;
2795  goto triangular_error;
2796  }
2797 
2798  double tmp = work[k] / data (mini);
2799  work[k] = tmp;
2800  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2801  {
2802  if (i == mini)
2803  continue;
2804 
2805  octave_idx_type iidx = perm[ridx (i)];
2806  work[iidx] = work[iidx] - tmp * data (i);
2807  }
2808  }
2809  }
2810 
2811  // Count nonzeros in work vector and adjust space in
2812  // retval if needed
2813  octave_idx_type new_nnz = 0;
2814  for (octave_idx_type i = 0; i < nc; i++)
2815  if (work[i] != 0.)
2816  new_nnz++;
2817 
2818  if (ii + new_nnz > x_nz)
2819  {
2820  // Resize the sparse matrix
2821  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2822  retval.change_capacity (sz);
2823  x_nz = sz;
2824  }
2825 
2826  for (octave_idx_type i = 0; i < nc; i++)
2827  if (work[i] != 0.)
2828  {
2829  retval.xridx (ii) = i;
2830  retval.xdata (ii++) = work[i];
2831  }
2832  retval.xcidx (j+1) = ii;
2833  }
2834 
2835  retval.maybe_compress ();
2836 
2837  if (calc_cond)
2838  {
2839  // Calculation of 1-norm of inv(*this)
2840  for (octave_idx_type i = 0; i < nm; i++)
2841  work[i] = 0.;
2842 
2843  for (octave_idx_type j = 0; j < nr; j++)
2844  {
2845  work[j] = 1.;
2846 
2847  for (octave_idx_type k = 0; k < nc; k++)
2848  {
2849  if (work[k] != 0.)
2850  {
2851  octave_idx_type minr = nr;
2852  octave_idx_type mini = 0;
2853 
2854  for (octave_idx_type i = cidx (k);
2855  i < cidx (k+1); i++)
2856  if (perm[ridx (i)] < minr)
2857  {
2858  minr = perm[ridx (i)];
2859  mini = i;
2860  }
2861 
2862  double tmp = work[k] / data (mini);
2863  work[k] = tmp;
2864  for (octave_idx_type i = cidx (k);
2865  i < cidx (k+1); i++)
2866  {
2867  if (i == mini)
2868  continue;
2869 
2870  octave_idx_type iidx = perm[ridx (i)];
2871  work[iidx] = work[iidx] - tmp * data (i);
2872  }
2873  }
2874  }
2875 
2876  double atmp = 0;
2877  for (octave_idx_type i = j; i < nr; i++)
2878  {
2879  atmp += fabs (work[i]);
2880  work[i] = 0.;
2881  }
2882  if (atmp > ainvnorm)
2883  ainvnorm = atmp;
2884  }
2885  rcond = 1. / ainvnorm / anorm;
2886  }
2887  }
2888  else
2889  {
2890  OCTAVE_LOCAL_BUFFER (double, work, nm);
2891 
2892  for (octave_idx_type j = 0; j < b_nc; j++)
2893  {
2894  for (octave_idx_type i = 0; i < nm; i++)
2895  work[i] = 0.;
2896  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2897  work[b.ridx (i)] = b.data (i);
2898 
2899  for (octave_idx_type k = 0; k < nc; k++)
2900  {
2901  if (work[k] != 0.)
2902  {
2903  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2904  {
2905  err = -2;
2906  goto triangular_error;
2907  }
2908 
2909  double tmp = work[k] / data (cidx (k));
2910  work[k] = tmp;
2911  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2912  {
2913  octave_idx_type iidx = ridx (i);
2914  work[iidx] = work[iidx] - tmp * data (i);
2915  }
2916  }
2917  }
2918 
2919  // Count nonzeros in work vector and adjust space in
2920  // retval if needed
2921  octave_idx_type new_nnz = 0;
2922  for (octave_idx_type i = 0; i < nc; i++)
2923  if (work[i] != 0.)
2924  new_nnz++;
2925 
2926  if (ii + new_nnz > x_nz)
2927  {
2928  // Resize the sparse matrix
2929  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2930  retval.change_capacity (sz);
2931  x_nz = sz;
2932  }
2933 
2934  for (octave_idx_type i = 0; i < nc; i++)
2935  if (work[i] != 0.)
2936  {
2937  retval.xridx (ii) = i;
2938  retval.xdata (ii++) = work[i];
2939  }
2940  retval.xcidx (j+1) = ii;
2941  }
2942 
2943  retval.maybe_compress ();
2944 
2945  if (calc_cond)
2946  {
2947  // Calculation of 1-norm of inv(*this)
2948  for (octave_idx_type i = 0; i < nm; i++)
2949  work[i] = 0.;
2950 
2951  for (octave_idx_type j = 0; j < nr; j++)
2952  {
2953  work[j] = 1.;
2954 
2955  for (octave_idx_type k = j; k < nc; k++)
2956  {
2957 
2958  if (work[k] != 0.)
2959  {
2960  double tmp = work[k] / data (cidx (k));
2961  work[k] = tmp;
2962  for (octave_idx_type i = cidx (k)+1;
2963  i < cidx (k+1); i++)
2964  {
2965  octave_idx_type iidx = ridx (i);
2966  work[iidx] = work[iidx] - tmp * data (i);
2967  }
2968  }
2969  }
2970  double atmp = 0;
2971  for (octave_idx_type i = j; i < nc; i++)
2972  {
2973  atmp += fabs (work[i]);
2974  work[i] = 0.;
2975  }
2976  if (atmp > ainvnorm)
2977  ainvnorm = atmp;
2978  }
2979  rcond = 1. / ainvnorm / anorm;
2980  }
2981  }
2982 
2983  triangular_error:
2984  if (err != 0)
2985  {
2986  if (sing_handler)
2987  {
2988  sing_handler (rcond);
2989  mattype.mark_as_rectangular ();
2990  }
2991  else
2993  }
2994 
2995  volatile double rcond_plus_one = rcond + 1.0;
2996 
2997  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
2998  {
2999  err = -2;
3000 
3001  if (sing_handler)
3002  {
3003  sing_handler (rcond);
3004  mattype.mark_as_rectangular ();
3005  }
3006  else
3008  }
3009  }
3010 
3011  return retval;
3012 }
3013 
3016  octave_idx_type& err, double& rcond,
3017  solve_singularity_handler sing_handler,
3018  bool calc_cond) const
3019 {
3021 
3022  octave_idx_type nr = rows ();
3023  octave_idx_type nc = cols ();
3024  octave_idx_type nm = (nc > nr ? nc : nr);
3025  err = 0;
3026 
3027  if (nr != b.rows ())
3029  ("matrix dimension mismatch solution of linear equations");
3030 
3031  if (nr == 0 || nc == 0 || b.cols () == 0)
3032  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3033  else
3034  {
3035  // Print spparms("spumoni") info if requested
3036  int typ = mattype.type ();
3037  mattype.info ();
3038 
3039  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3040  (*current_liboctave_error_handler) ("incorrect matrix type");
3041 
3042  double anorm = 0.;
3043  double ainvnorm = 0.;
3044  octave_idx_type b_nc = b.cols ();
3045  rcond = 1.;
3046 
3047  if (calc_cond)
3048  {
3049  // Calculate the 1-norm of matrix for rcond calculation
3050  for (octave_idx_type j = 0; j < nc; j++)
3051  {
3052  double atmp = 0.;
3053  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3054  atmp += fabs (data (i));
3055  if (atmp > anorm)
3056  anorm = atmp;
3057  }
3058  }
3059 
3060  if (typ == MatrixType::Permuted_Lower)
3061  {
3062  retval.resize (nc, b_nc);
3063  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3064  octave_idx_type *perm = mattype.triangular_perm ();
3065 
3066  for (octave_idx_type j = 0; j < b_nc; j++)
3067  {
3068  for (octave_idx_type i = 0; i < nm; i++)
3069  cwork[i] = 0.;
3070  for (octave_idx_type i = 0; i < nr; i++)
3071  cwork[perm[i]] = b(i,j);
3072 
3073  for (octave_idx_type k = 0; k < nc; k++)
3074  {
3075  if (cwork[k] != 0.)
3076  {
3077  octave_idx_type minr = nr;
3078  octave_idx_type mini = 0;
3079 
3080  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3081  if (perm[ridx (i)] < minr)
3082  {
3083  minr = perm[ridx (i)];
3084  mini = i;
3085  }
3086 
3087  if (minr != k || data (mini) == 0)
3088  {
3089  err = -2;
3090  goto triangular_error;
3091  }
3092 
3093  Complex tmp = cwork[k] / data (mini);
3094  cwork[k] = tmp;
3095  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3096  {
3097  if (i == mini)
3098  continue;
3099 
3100  octave_idx_type iidx = perm[ridx (i)];
3101  cwork[iidx] = cwork[iidx] - tmp * data (i);
3102  }
3103  }
3104  }
3105 
3106  for (octave_idx_type i = 0; i < nc; i++)
3107  retval(i, j) = cwork[i];
3108  }
3109 
3110  if (calc_cond)
3111  {
3112  // Calculation of 1-norm of inv(*this)
3113  OCTAVE_LOCAL_BUFFER (double, work, nm);
3114  for (octave_idx_type i = 0; i < nm; i++)
3115  work[i] = 0.;
3116 
3117  for (octave_idx_type j = 0; j < nr; j++)
3118  {
3119  work[j] = 1.;
3120 
3121  for (octave_idx_type k = 0; k < nc; k++)
3122  {
3123  if (work[k] != 0.)
3124  {
3125  octave_idx_type minr = nr;
3126  octave_idx_type mini = 0;
3127 
3128  for (octave_idx_type i = cidx (k);
3129  i < cidx (k+1); i++)
3130  if (perm[ridx (i)] < minr)
3131  {
3132  minr = perm[ridx (i)];
3133  mini = i;
3134  }
3135 
3136  double tmp = work[k] / data (mini);
3137  work[k] = tmp;
3138  for (octave_idx_type i = cidx (k);
3139  i < cidx (k+1); i++)
3140  {
3141  if (i == mini)
3142  continue;
3143 
3144  octave_idx_type iidx = perm[ridx (i)];
3145  work[iidx] = work[iidx] - tmp * data (i);
3146  }
3147  }
3148  }
3149 
3150  double atmp = 0;
3151  for (octave_idx_type i = j; i < nc; i++)
3152  {
3153  atmp += fabs (work[i]);
3154  work[i] = 0.;
3155  }
3156  if (atmp > ainvnorm)
3157  ainvnorm = atmp;
3158  }
3159  rcond = 1. / ainvnorm / anorm;
3160  }
3161  }
3162  else
3163  {
3164  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3165  retval.resize (nc, b_nc, 0.);
3166 
3167  for (octave_idx_type j = 0; j < b_nc; j++)
3168  {
3169  for (octave_idx_type i = 0; i < nr; i++)
3170  cwork[i] = b(i,j);
3171  for (octave_idx_type i = nr; i < nc; i++)
3172  cwork[i] = 0.;
3173 
3174  for (octave_idx_type k = 0; k < nc; k++)
3175  {
3176  if (cwork[k] != 0.)
3177  {
3178  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3179  {
3180  err = -2;
3181  goto triangular_error;
3182  }
3183 
3184  Complex tmp = cwork[k] / data (cidx (k));
3185  cwork[k] = tmp;
3186  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3187  {
3188  octave_idx_type iidx = ridx (i);
3189  cwork[iidx] = cwork[iidx] - tmp * data (i);
3190  }
3191  }
3192  }
3193 
3194  for (octave_idx_type i = 0; i < nc; i++)
3195  retval.xelem (i, j) = cwork[i];
3196  }
3197 
3198  if (calc_cond)
3199  {
3200  // Calculation of 1-norm of inv(*this)
3201  OCTAVE_LOCAL_BUFFER (double, work, nm);
3202  for (octave_idx_type i = 0; i < nm; i++)
3203  work[i] = 0.;
3204 
3205  for (octave_idx_type j = 0; j < nr; j++)
3206  {
3207  work[j] = 1.;
3208 
3209  for (octave_idx_type k = j; k < nc; k++)
3210  {
3211 
3212  if (work[k] != 0.)
3213  {
3214  double tmp = work[k] / data (cidx (k));
3215  work[k] = tmp;
3216  for (octave_idx_type i = cidx (k)+1;
3217  i < cidx (k+1); i++)
3218  {
3219  octave_idx_type iidx = ridx (i);
3220  work[iidx] = work[iidx] - tmp * data (i);
3221  }
3222  }
3223  }
3224  double atmp = 0;
3225  for (octave_idx_type i = j; i < nc; i++)
3226  {
3227  atmp += fabs (work[i]);
3228  work[i] = 0.;
3229  }
3230  if (atmp > ainvnorm)
3231  ainvnorm = atmp;
3232  }
3233  rcond = 1. / ainvnorm / anorm;
3234  }
3235  }
3236 
3237  triangular_error:
3238  if (err != 0)
3239  {
3240  if (sing_handler)
3241  {
3242  sing_handler (rcond);
3243  mattype.mark_as_rectangular ();
3244  }
3245  else
3247  }
3248 
3249  volatile double rcond_plus_one = rcond + 1.0;
3250 
3251  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
3252  {
3253  err = -2;
3254 
3255  if (sing_handler)
3256  {
3257  sing_handler (rcond);
3258  mattype.mark_as_rectangular ();
3259  }
3260  else
3262  }
3263  }
3264 
3265  return retval;
3266 }
3267 
3270  octave_idx_type& err, double& rcond,
3271  solve_singularity_handler sing_handler,
3272  bool calc_cond) const
3273 {
3275 
3276  octave_idx_type nr = rows ();
3277  octave_idx_type nc = cols ();
3278  octave_idx_type nm = (nc > nr ? nc : nr);
3279  err = 0;
3280 
3281  if (nr != b.rows ())
3283  ("matrix dimension mismatch solution of linear equations");
3284 
3285  if (nr == 0 || nc == 0 || b.cols () == 0)
3286  retval = SparseComplexMatrix (nc, b.cols ());
3287  else
3288  {
3289  // Print spparms("spumoni") info if requested
3290  int typ = mattype.type ();
3291  mattype.info ();
3292 
3293  if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3294  (*current_liboctave_error_handler) ("incorrect matrix type");
3295 
3296  double anorm = 0.;
3297  double ainvnorm = 0.;
3298  rcond = 1.;
3299 
3300  if (calc_cond)
3301  {
3302  // Calculate the 1-norm of matrix for rcond calculation
3303  for (octave_idx_type j = 0; j < nc; j++)
3304  {
3305  double atmp = 0.;
3306  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3307  atmp += fabs (data (i));
3308  if (atmp > anorm)
3309  anorm = atmp;
3310  }
3311  }
3312 
3313  octave_idx_type b_nc = b.cols ();
3314  octave_idx_type b_nz = b.nnz ();
3315  retval = SparseComplexMatrix (nc, b_nc, b_nz);
3316  retval.xcidx (0) = 0;
3317  octave_idx_type ii = 0;
3318  octave_idx_type x_nz = b_nz;
3319 
3320  if (typ == MatrixType::Permuted_Lower)
3321  {
3322  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3323  octave_idx_type *perm = mattype.triangular_perm ();
3324 
3325  for (octave_idx_type j = 0; j < b_nc; j++)
3326  {
3327  for (octave_idx_type i = 0; i < nm; i++)
3328  cwork[i] = 0.;
3329  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3330  cwork[perm[b.ridx (i)]] = b.data (i);
3331 
3332  for (octave_idx_type k = 0; k < nc; k++)
3333  {
3334  if (cwork[k] != 0.)
3335  {
3336  octave_idx_type minr = nr;
3337  octave_idx_type mini = 0;
3338 
3339  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3340  if (perm[ridx (i)] < minr)
3341  {
3342  minr = perm[ridx (i)];
3343  mini = i;
3344  }
3345 
3346  if (minr != k || data (mini) == 0)
3347  {
3348  err = -2;
3349  goto triangular_error;
3350  }
3351 
3352  Complex tmp = cwork[k] / data (mini);
3353  cwork[k] = tmp;
3354  for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3355  {
3356  if (i == mini)
3357  continue;
3358 
3359  octave_idx_type iidx = perm[ridx (i)];
3360  cwork[iidx] = cwork[iidx] - tmp * data (i);
3361  }
3362  }
3363  }
3364 
3365  // Count nonzeros in work vector and adjust space in
3366  // retval if needed
3367  octave_idx_type new_nnz = 0;
3368  for (octave_idx_type i = 0; i < nc; i++)
3369  if (cwork[i] != 0.)
3370  new_nnz++;
3371 
3372  if (ii + new_nnz > x_nz)
3373  {
3374  // Resize the sparse matrix
3375  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3376  retval.change_capacity (sz);
3377  x_nz = sz;
3378  }
3379 
3380  for (octave_idx_type i = 0; i < nc; i++)
3381  if (cwork[i] != 0.)
3382  {
3383  retval.xridx (ii) = i;
3384  retval.xdata (ii++) = cwork[i];
3385  }
3386  retval.xcidx (j+1) = ii;
3387  }
3388 
3389  retval.maybe_compress ();
3390 
3391  if (calc_cond)
3392  {
3393  // Calculation of 1-norm of inv(*this)
3394  OCTAVE_LOCAL_BUFFER (double, work, nm);
3395  for (octave_idx_type i = 0; i < nm; i++)
3396  work[i] = 0.;
3397 
3398  for (octave_idx_type j = 0; j < nr; j++)
3399  {
3400  work[j] = 1.;
3401 
3402  for (octave_idx_type k = 0; k < nc; k++)
3403  {
3404  if (work[k] != 0.)
3405  {
3406  octave_idx_type minr = nr;
3407  octave_idx_type mini = 0;
3408 
3409  for (octave_idx_type i = cidx (k);
3410  i < cidx (k+1); i++)
3411  if (perm[ridx (i)] < minr)
3412  {
3413  minr = perm[ridx (i)];
3414  mini = i;
3415  }
3416 
3417  double tmp = work[k] / data (mini);
3418  work[k] = tmp;
3419  for (octave_idx_type i = cidx (k);
3420  i < cidx (k+1); i++)
3421  {
3422  if (i == mini)
3423  continue;
3424 
3425  octave_idx_type iidx = perm[ridx (i)];
3426  work[iidx] = work[iidx] - tmp * data (i);
3427  }
3428  }
3429  }
3430 
3431  double atmp = 0;
3432  for (octave_idx_type i = j; i < nc; i++)
3433  {
3434  atmp += fabs (work[i]);
3435  work[i] = 0.;
3436  }
3437  if (atmp > ainvnorm)
3438  ainvnorm = atmp;
3439  }
3440  rcond = 1. / ainvnorm / anorm;
3441  }
3442  }
3443  else
3444  {
3445  OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3446 
3447  for (octave_idx_type j = 0; j < b_nc; j++)
3448  {
3449  for (octave_idx_type i = 0; i < nm; i++)
3450  cwork[i] = 0.;
3451  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3452  cwork[b.ridx (i)] = b.data (i);
3453 
3454  for (octave_idx_type k = 0; k < nc; k++)
3455  {
3456  if (cwork[k] != 0.)
3457  {
3458  if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3459  {
3460  err = -2;
3461  goto triangular_error;
3462  }
3463 
3464  Complex tmp = cwork[k] / data (cidx (k));
3465  cwork[k] = tmp;
3466  for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3467  {
3468  octave_idx_type iidx = ridx (i);
3469  cwork[iidx] = cwork[iidx] - tmp * data (i);
3470  }
3471  }
3472  }
3473 
3474  // Count nonzeros in work vector and adjust space in
3475  // retval if needed
3476  octave_idx_type new_nnz = 0;
3477  for (octave_idx_type i = 0; i < nc; i++)
3478  if (cwork[i] != 0.)
3479  new_nnz++;
3480 
3481  if (ii + new_nnz > x_nz)
3482  {
3483  // Resize the sparse matrix
3484  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3485  retval.change_capacity (sz);
3486  x_nz = sz;
3487  }
3488 
3489  for (octave_idx_type i = 0; i < nc; i++)
3490  if (cwork[i] != 0.)
3491  {
3492  retval.xridx (ii) = i;
3493  retval.xdata (ii++) = cwork[i];
3494  }
3495  retval.xcidx (j+1) = ii;
3496  }
3497 
3498  retval.maybe_compress ();
3499 
3500  if (calc_cond)
3501  {
3502  // Calculation of 1-norm of inv(*this)
3503  OCTAVE_LOCAL_BUFFER (double, work, nm);
3504  for (octave_idx_type i = 0; i < nm; i++)
3505  work[i] = 0.;
3506 
3507  for (octave_idx_type j = 0; j < nr; j++)
3508  {
3509  work[j] = 1.;
3510 
3511  for (octave_idx_type k = j; k < nc; k++)
3512  {
3513 
3514  if (work[k] != 0.)
3515  {
3516  double tmp = work[k] / data (cidx (k));
3517  work[k] = tmp;
3518  for (octave_idx_type i = cidx (k)+1;
3519  i < cidx (k+1); i++)
3520  {
3521  octave_idx_type iidx = ridx (i);
3522  work[iidx] = work[iidx] - tmp * data (i);
3523  }
3524  }
3525  }
3526  double atmp = 0;
3527  for (octave_idx_type i = j; i < nc; i++)
3528  {
3529  atmp += fabs (work[i]);
3530  work[i] = 0.;
3531  }
3532  if (atmp > ainvnorm)
3533  ainvnorm = atmp;
3534  }
3535  rcond = 1. / ainvnorm / anorm;
3536  }
3537  }
3538 
3539  triangular_error:
3540  if (err != 0)
3541  {
3542  if (sing_handler)
3543  {
3544  sing_handler (rcond);
3545  mattype.mark_as_rectangular ();
3546  }
3547  else
3549  }
3550 
3551  volatile double rcond_plus_one = rcond + 1.0;
3552 
3553  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
3554  {
3555  err = -2;
3556 
3557  if (sing_handler)
3558  {
3559  sing_handler (rcond);
3560  mattype.mark_as_rectangular ();
3561  }
3562  else
3564  }
3565  }
3566 
3567  return retval;
3568 }
3569 
3570 Matrix
3572  octave_idx_type& err, double& rcond,
3573  solve_singularity_handler sing_handler,
3574  bool calc_cond) const
3575 {
3576  Matrix retval;
3577 
3578  octave_idx_type nr = rows ();
3579  octave_idx_type nc = cols ();
3580  err = 0;
3581 
3582  if (nr != nc || nr != b.rows ())
3584  ("matrix dimension mismatch solution of linear equations");
3585 
3586  if (nr == 0 || b.cols () == 0)
3587  retval = Matrix (nc, b.cols (), 0.0);
3588  else if (calc_cond)
3589  (*current_liboctave_error_handler)
3590  ("calculation of condition number not implemented");
3591  else
3592  {
3593  // Print spparms("spumoni") info if requested
3594  volatile int typ = mattype.type ();
3595  mattype.info ();
3596 
3598  {
3599  OCTAVE_LOCAL_BUFFER (double, D, nr);
3600  OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
3601 
3602  if (mattype.is_dense ())
3603  {
3604  octave_idx_type ii = 0;
3605 
3606  for (octave_idx_type j = 0; j < nc-1; j++)
3607  {
3608  D[j] = data (ii++);
3609  DL[j] = data (ii);
3610  ii += 2;
3611  }
3612  D[nc-1] = data (ii);
3613  }
3614  else
3615  {
3616  D[0] = 0.;
3617  for (octave_idx_type i = 0; i < nr - 1; i++)
3618  {
3619  D[i+1] = 0.;
3620  DL[i] = 0.;
3621  }
3622 
3623  for (octave_idx_type j = 0; j < nc; j++)
3624  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3625  {
3626  if (ridx (i) == j)
3627  D[j] = data (i);
3628  else if (ridx (i) == j + 1)
3629  DL[j] = data (i);
3630  }
3631  }
3632 
3633  F77_INT tmp_nr = octave::to_f77_int (nr);
3634 
3635  F77_INT b_nr = octave::to_f77_int (b.rows ());
3636  F77_INT b_nc = octave::to_f77_int (b.cols ());
3637 
3638  retval = b;
3639  double *result = retval.fortran_vec ();
3640 
3641  F77_INT tmp_err = 0;
3642 
3643  F77_XFCN (dptsv, DPTSV, (tmp_nr, b_nc, D, DL, result,
3644  b_nr, tmp_err));
3645 
3646  err = tmp_err;
3647 
3648  if (err != 0)
3649  {
3650  err = 0;
3651  mattype.mark_as_unsymmetric ();
3653  }
3654  else
3655  rcond = 1.;
3656  }
3657 
3658  if (typ == MatrixType::Tridiagonal)
3659  {
3660  OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
3661  OCTAVE_LOCAL_BUFFER (double, D, nr);
3662  OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
3663 
3664  if (mattype.is_dense ())
3665  {
3666  octave_idx_type ii = 0;
3667 
3668  for (octave_idx_type j = 0; j < nc-1; j++)
3669  {
3670  D[j] = data (ii++);
3671  DL[j] = data (ii++);
3672  DU[j] = data (ii++);
3673  }
3674  D[nc-1] = data (ii);
3675  }
3676  else
3677  {
3678  D[0] = 0.;
3679  for (octave_idx_type i = 0; i < nr - 1; i++)
3680  {
3681  D[i+1] = 0.;
3682  DL[i] = 0.;
3683  DU[i] = 0.;
3684  }
3685 
3686  for (octave_idx_type j = 0; j < nc; j++)
3687  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3688  {
3689  if (ridx (i) == j)
3690  D[j] = data (i);
3691  else if (ridx (i) == j + 1)
3692  DL[j] = data (i);
3693  else if (ridx (i) == j - 1)
3694  DU[j-1] = data (i);
3695  }
3696  }
3697 
3698  F77_INT tmp_nr = octave::to_f77_int (nr);
3699 
3700  F77_INT b_nr = octave::to_f77_int (b.rows ());
3701  F77_INT b_nc = octave::to_f77_int (b.cols ());
3702 
3703  retval = b;
3704  double *result = retval.fortran_vec ();
3705 
3706  F77_INT tmp_err = 0;
3707 
3708  F77_XFCN (dgtsv, DGTSV, (tmp_nr, b_nc, DL, D, DU, result,
3709  b_nr, tmp_err));
3710 
3711  err = tmp_err;
3712 
3713  if (err != 0)
3714  {
3715  rcond = 0.;
3716  err = -2;
3717 
3718  if (sing_handler)
3719  {
3720  sing_handler (rcond);
3721  mattype.mark_as_rectangular ();
3722  }
3723  else
3725  }
3726  else
3727  rcond = 1.;
3728  }
3729  else if (typ != MatrixType::Tridiagonal_Hermitian)
3730  (*current_liboctave_error_handler) ("incorrect matrix type");
3731  }
3732 
3733  return retval;
3734 }
3735 
3738  octave_idx_type& err, double& rcond,
3739  solve_singularity_handler sing_handler,
3740  bool calc_cond) const
3741 {
3743 
3744  octave_idx_type nr = rows ();
3745  octave_idx_type nc = cols ();
3746  err = 0;
3747 
3748  if (nr != nc || nr != b.rows ())
3750  ("matrix dimension mismatch solution of linear equations");
3751 
3752  if (nr == 0 || b.cols () == 0)
3753  retval = SparseMatrix (nc, b.cols ());
3754  else if (calc_cond)
3755  (*current_liboctave_error_handler)
3756  ("calculation of condition number not implemented");
3757  else
3758  {
3759  // Print spparms("spumoni") info if requested
3760  int typ = mattype.type ();
3761  mattype.info ();
3762 
3763  // Note can't treat symmetric case as there is no dpttrf function
3764  if (typ == MatrixType::Tridiagonal
3766  {
3767  OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
3768  OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
3769  OCTAVE_LOCAL_BUFFER (double, D, nr);
3770  OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
3771  Array<F77_INT> ipvt (dim_vector (nr, 1));
3772  F77_INT *pipvt = ipvt.fortran_vec ();
3773 
3774  if (mattype.is_dense ())
3775  {
3776  octave_idx_type ii = 0;
3777 
3778  for (octave_idx_type j = 0; j < nc-1; j++)
3779  {
3780  D[j] = data (ii++);
3781  DL[j] = data (ii++);
3782  DU[j] = data (ii++);
3783  }
3784  D[nc-1] = data (ii);
3785  }
3786  else
3787  {
3788  D[0] = 0.;
3789  for (octave_idx_type i = 0; i < nr - 1; i++)
3790  {
3791  D[i+1] = 0.;
3792  DL[i] = 0.;
3793  DU[i] = 0.;
3794  }
3795 
3796  for (octave_idx_type j = 0; j < nc; j++)
3797  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3798  {
3799  if (ridx (i) == j)
3800  D[j] = data (i);
3801  else if (ridx (i) == j + 1)
3802  DL[j] = data (i);
3803  else if (ridx (i) == j - 1)
3804  DU[j-1] = data (i);
3805  }
3806  }
3807 
3808  F77_INT tmp_nr = octave::to_f77_int (nr);
3809 
3810  F77_INT tmp_err = 0;
3811 
3812  F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
3813 
3814  if (err != 0)
3815  {
3816  rcond = 0.0;
3817  err = -2;
3818 
3819  if (sing_handler)
3820  {
3821  sing_handler (rcond);
3822  mattype.mark_as_rectangular ();
3823  }
3824  else
3826  }
3827  else
3828  {
3829  rcond = 1.0;
3830  char job = 'N';
3831  volatile octave_idx_type x_nz = b.nnz ();
3832  octave_idx_type b_nc = b.cols ();
3833  retval = SparseMatrix (nr, b_nc, x_nz);
3834  retval.xcidx (0) = 0;
3835  volatile octave_idx_type ii = 0;
3836 
3837  OCTAVE_LOCAL_BUFFER (double, work, nr);
3838 
3839  for (volatile octave_idx_type j = 0; j < b_nc; j++)
3840  {
3841  for (octave_idx_type i = 0; i < nr; i++)
3842  work[i] = 0.;
3843  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3844  work[b.ridx (i)] = b.data (i);
3845 
3846  F77_INT b_nr = octave::to_f77_int (b.rows ());
3847 
3848  F77_XFCN (dgttrs, DGTTRS,
3849  (F77_CONST_CHAR_ARG2 (&job, 1),
3850  tmp_nr, 1, DL, D, DU, DU2, pipvt,
3851  work, b_nr, tmp_err
3852  F77_CHAR_ARG_LEN (1)));
3853 
3854  err = tmp_err;
3855 
3856  // Count nonzeros in work vector and adjust
3857  // space in retval if needed
3858  octave_idx_type new_nnz = 0;
3859  for (octave_idx_type i = 0; i < nr; i++)
3860  if (work[i] != 0.)
3861  new_nnz++;
3862 
3863  if (ii + new_nnz > x_nz)
3864  {
3865  // Resize the sparse matrix
3866  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3867  retval.change_capacity (sz);
3868  x_nz = sz;
3869  }
3870 
3871  for (octave_idx_type i = 0; i < nr; i++)
3872  if (work[i] != 0.)
3873  {
3874  retval.xridx (ii) = i;
3875  retval.xdata (ii++) = work[i];
3876  }
3877  retval.xcidx (j+1) = ii;
3878  }
3879 
3880  retval.maybe_compress ();
3881  }
3882  }
3883  else if (typ != MatrixType::Tridiagonal_Hermitian)
3884  (*current_liboctave_error_handler) ("incorrect matrix type");
3885  }
3886 
3887  return retval;
3888 }
3889 
3892  octave_idx_type& err, double& rcond,
3893  solve_singularity_handler sing_handler,
3894  bool calc_cond) const
3895 {
3897 
3898  octave_idx_type nr = rows ();
3899  octave_idx_type nc = cols ();
3900  err = 0;
3901 
3902  if (nr != nc || nr != b.rows ())
3904  ("matrix dimension mismatch solution of linear equations");
3905 
3906  if (nr == 0 || b.cols () == 0)
3907  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3908  else if (calc_cond)
3909  (*current_liboctave_error_handler)
3910  ("calculation of condition number not implemented");
3911  else
3912  {
3913  // Print spparms("spumoni") info if requested
3914  volatile int typ = mattype.type ();
3915  mattype.info ();
3916 
3918  {
3919  OCTAVE_LOCAL_BUFFER (double, D, nr);
3920  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3921 
3922  if (mattype.is_dense ())
3923  {
3924  octave_idx_type ii = 0;
3925 
3926  for (octave_idx_type j = 0; j < nc-1; j++)
3927  {
3928  D[j] = data (ii++);
3929  DL[j] = data (ii);
3930  ii += 2;
3931  }
3932  D[nc-1] = data (ii);
3933  }
3934  else
3935  {
3936  D[0] = 0.;
3937  for (octave_idx_type i = 0; i < nr - 1; i++)
3938  {
3939  D[i+1] = 0.;
3940  DL[i] = 0.;
3941  }
3942 
3943  for (octave_idx_type j = 0; j < nc; j++)
3944  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3945  {
3946  if (ridx (i) == j)
3947  D[j] = data (i);
3948  else if (ridx (i) == j + 1)
3949  DL[j] = data (i);
3950  }
3951  }
3952 
3953  F77_INT tmp_nr = octave::to_f77_int (nr);
3954 
3955  F77_INT b_nr = octave::to_f77_int (b.rows ());
3956  F77_INT b_nc = octave::to_f77_int (b.cols ());
3957 
3958  rcond = 1.;
3959 
3960  retval = b;
3961  Complex *result = retval.fortran_vec ();
3962 
3963  F77_INT tmp_err = 0;
3964 
3965  F77_XFCN (zptsv, ZPTSV, (tmp_nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL),
3967  b_nr, tmp_err));
3968 
3969  err = tmp_err;
3970 
3971  if (err != 0)
3972  {
3973  err = 0;
3974  mattype.mark_as_unsymmetric ();
3976  }
3977  }
3978 
3979  if (typ == MatrixType::Tridiagonal)
3980  {
3981  OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3982  OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3983  OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3984 
3985  if (mattype.is_dense ())
3986  {
3987  octave_idx_type ii = 0;
3988 
3989  for (octave_idx_type j = 0; j < nc-1; j++)
3990  {
3991  D[j] = data (ii++);
3992  DL[j] = data (ii++);
3993  DU[j] = data (ii++);
3994  }
3995  D[nc-1] = data (ii);
3996  }
3997  else
3998  {
3999  D[0] = 0.;
4000  for (octave_idx_type i = 0; i < nr - 1; i++)
4001  {
4002  D[i+1] = 0.;
4003  DL[i] = 0.;
4004  DU[i] = 0.;
4005  }
4006 
4007  for (octave_idx_type j = 0; j < nc; j++)
4008  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4009  {
4010  if (ridx (i) == j)
4011  D[j] = data (i);
4012  else if (ridx (i) == j + 1)
4013  DL[j] = data (i);
4014  else if (ridx (i) == j - 1)
4015  DU[j-1] = data (i);
4016  }
4017  }
4018 
4019  F77_INT tmp_nr = octave::to_f77_int (nr);
4020 
4021  F77_INT b_nr = octave::to_f77_int (b.rows ());
4022  F77_INT b_nc = octave::to_f77_int (b.cols ());
4023 
4024  rcond = 1.;
4025 
4026  retval = b;
4027  Complex *result = retval.fortran_vec ();
4028 
4029  F77_INT tmp_err = 0;
4030 
4031  F77_XFCN (zgtsv, ZGTSV, (tmp_nr, b_nc, F77_DBLE_CMPLX_ARG (DL),
4032  F77_DBLE_CMPLX_ARG (D),
4033  F77_DBLE_CMPLX_ARG (DU),
4035  b_nr, tmp_err));
4036 
4037  err = tmp_err;
4038 
4039  if (err != 0)
4040  {
4041  rcond = 0.;
4042  err = -2;
4043 
4044  if (sing_handler)
4045  {
4046  sing_handler (rcond);
4047  mattype.mark_as_rectangular ();
4048  }
4049  else
4051  }
4052  }
4053  else if (typ != MatrixType::Tridiagonal_Hermitian)
4054  (*current_liboctave_error_handler) ("incorrect matrix type");
4055  }
4056 
4057  return retval;
4058 }
4059 
4062  octave_idx_type& err, double& rcond,
4063  solve_singularity_handler sing_handler,
4064  bool calc_cond) const
4065 {
4067 
4068  octave_idx_type nr = rows ();
4069  octave_idx_type nc = cols ();
4070  err = 0;
4071 
4072  if (nr != nc || nr != b.rows ())
4074  ("matrix dimension mismatch solution of linear equations");
4075 
4076  if (nr == 0 || b.cols () == 0)
4077  retval = SparseComplexMatrix (nc, b.cols ());
4078  else if (calc_cond)
4079  (*current_liboctave_error_handler)
4080  ("calculation of condition number not implemented");
4081  else
4082  {
4083  // Print spparms("spumoni") info if requested
4084  int typ = mattype.type ();
4085  mattype.info ();
4086 
4087  // Note can't treat symmetric case as there is no dpttrf function
4088  if (typ == MatrixType::Tridiagonal
4090  {
4091  OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
4092  OCTAVE_LOCAL_BUFFER (double, DU, nr - 1);
4093  OCTAVE_LOCAL_BUFFER (double, D, nr);
4094  OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
4095  Array<F77_INT> ipvt (dim_vector (nr, 1));
4096  F77_INT *pipvt = ipvt.fortran_vec ();
4097 
4098  if (mattype.is_dense ())
4099  {
4100  octave_idx_type ii = 0;
4101 
4102  for (octave_idx_type j = 0; j < nc-1; j++)
4103  {
4104  D[j] = data (ii++);
4105  DL[j] = data (ii++);
4106  DU[j] = data (ii++);
4107  }
4108  D[nc-1] = data (ii);
4109  }
4110  else
4111  {
4112  D[0] = 0.;
4113  for (octave_idx_type i = 0; i < nr - 1; i++)
4114  {
4115  D[i+1] = 0.;
4116  DL[i] = 0.;
4117  DU[i] = 0.;
4118  }
4119 
4120  for (octave_idx_type j = 0; j < nc; j++)
4121  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4122  {
4123  if (ridx (i) == j)
4124  D[j] = data (i);
4125  else if (ridx (i) == j + 1)
4126  DL[j] = data (i);
4127  else if (ridx (i) == j - 1)
4128  DU[j-1] = data (i);
4129  }
4130  }
4131 
4132  F77_INT tmp_nr = octave::to_f77_int (nr);
4133 
4134  F77_INT tmp_err = 0;
4135 
4136  F77_XFCN (dgttrf, DGTTRF, (tmp_nr, DL, D, DU, DU2, pipvt, tmp_err));
4137 
4138  err = tmp_err;
4139 
4140  if (err != 0)
4141  {
4142  rcond = 0.0;
4143  err = -2;
4144 
4145  if (sing_handler)
4146  {
4147  sing_handler (rcond);
4148  mattype.mark_as_rectangular ();
4149  }
4150  else
4152  }
4153  else
4154  {
4155  rcond = 1.;
4156  char job = 'N';
4157  F77_INT b_nr = octave::to_f77_int (b.rows ());
4158  octave_idx_type b_nc = b.cols ();
4159  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4160  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
4161 
4162  // Take a first guess that the number of nonzero terms
4163  // will be as many as in b
4164  volatile octave_idx_type x_nz = b.nnz ();
4165  volatile octave_idx_type ii = 0;
4166  retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4167 
4168  retval.xcidx (0) = 0;
4169  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4170  {
4171 
4172  for (F77_INT i = 0; i < b_nr; i++)
4173  {
4174  Complex c = b(i,j);
4175  Bx[i] = c.real ();
4176  Bz[i] = c.imag ();
4177  }
4178 
4179  F77_XFCN (dgttrs, DGTTRS,
4180  (F77_CONST_CHAR_ARG2 (&job, 1),
4181  tmp_nr, 1, DL, D, DU, DU2, pipvt,
4182  Bx, b_nr, tmp_err
4183  F77_CHAR_ARG_LEN (1)));
4184 
4185  err = tmp_err;
4186 
4187  if (err != 0)
4188  {
4189  // FIXME: Should this be a warning?
4190  (*current_liboctave_error_handler)
4191  ("SparseMatrix::solve solve failed");
4192 
4193  err = -1;
4194  break;
4195  }
4196 
4197  F77_XFCN (dgttrs, DGTTRS,
4198  (F77_CONST_CHAR_ARG2 (&job, 1),
4199  tmp_nr, 1, DL, D, DU, DU2, pipvt,
4200  Bz, b_nr, tmp_err
4201  F77_CHAR_ARG_LEN (1)));
4202 
4203  err = tmp_err;
4204 
4205  if (err != 0)
4206  {
4207  // FIXME: Should this be a warning?
4208  (*current_liboctave_error_handler)
4209  ("SparseMatrix::solve solve failed");
4210 
4211  err = -1;
4212  break;
4213  }
4214 
4215  // Count nonzeros in work vector and adjust
4216  // space in retval if needed
4217  octave_idx_type new_nnz = 0;
4218  for (octave_idx_type i = 0; i < nr; i++)
4219  if (Bx[i] != 0. || Bz[i] != 0.)
4220  new_nnz++;
4221 
4222  if (ii + new_nnz > x_nz)
4223  {
4224  // Resize the sparse matrix
4225  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4226  retval.change_capacity (sz);
4227  x_nz = sz;
4228  }
4229 
4230  for (octave_idx_type i = 0; i < nr; i++)
4231  if (Bx[i] != 0. || Bz[i] != 0.)
4232  {
4233  retval.xridx (ii) = i;
4234  retval.xdata (ii++) =
4235  Complex (Bx[i], Bz[i]);
4236  }
4237 
4238  retval.xcidx (j+1) = ii;
4239  }
4240 
4241  retval.maybe_compress ();
4242  }
4243  }
4244  else if (typ != MatrixType::Tridiagonal_Hermitian)
4245  (*current_liboctave_error_handler) ("incorrect matrix type");
4246  }
4247 
4248  return retval;
4249 }
4250 
4251 Matrix
4253  octave_idx_type& err, double& rcond,
4254  solve_singularity_handler sing_handler,
4255  bool calc_cond) const
4256 {
4257  Matrix retval;
4258 
4259  octave_idx_type nr = rows ();
4260  octave_idx_type nc = cols ();
4261  err = 0;
4262 
4263  if (nr != nc || nr != b.rows ())
4265  ("matrix dimension mismatch solution of linear equations");
4266 
4267  if (nr == 0 || b.cols () == 0)
4268  retval = Matrix (nc, b.cols (), 0.0);
4269  else
4270  {
4271  // Print spparms("spumoni") info if requested
4272  volatile int typ = mattype.type ();
4273  mattype.info ();
4274 
4275  if (typ == MatrixType::Banded_Hermitian)
4276  {
4277  F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4278  F77_INT ldm = n_lower + 1;
4279  Matrix m_band (ldm, nc);
4280  double *tmp_data = m_band.fortran_vec ();
4281 
4282  if (! mattype.is_dense ())
4283  {
4284  octave_idx_type ii = 0;
4285 
4286  for (octave_idx_type j = 0; j < ldm; j++)
4287  for (octave_idx_type i = 0; i < nc; i++)
4288  tmp_data[ii++] = 0.;
4289  }
4290 
4291  for (octave_idx_type j = 0; j < nc; j++)
4292  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4293  {
4294  octave_idx_type ri = ridx (i);
4295  if (ri >= j)
4296  m_band(ri - j, j) = data (i);
4297  }
4298 
4299  // Calculate the norm of the matrix, for later use.
4300  double anorm;
4301  if (calc_cond)
4302  anorm = m_band.abs ().sum ().row (0).max ();
4303 
4304  F77_INT tmp_nr = octave::to_f77_int (nr);
4305 
4306  F77_INT tmp_err = 0;
4307 
4308  char job = 'L';
4309  F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4310  tmp_nr, n_lower, tmp_data, ldm, tmp_err
4311  F77_CHAR_ARG_LEN (1)));
4312 
4313  err = tmp_err;
4314 
4315  if (err != 0)
4316  {
4317  // Matrix is not positive definite!! Fall through to
4318  // unsymmetric banded solver.
4319  mattype.mark_as_unsymmetric ();
4320  typ = MatrixType::Banded;
4321  rcond = 0.0;
4322  err = 0;
4323  }
4324  else
4325  {
4326  if (calc_cond)
4327  {
4328  Array<double> z (dim_vector (3 * nr, 1));
4329  double *pz = z.fortran_vec ();
4330  Array<F77_INT> iz (dim_vector (nr, 1));
4331  F77_INT *piz = iz.fortran_vec ();
4332 
4333  F77_XFCN (dpbcon, DPBCON,
4334  (F77_CONST_CHAR_ARG2 (&job, 1),
4335  tmp_nr, n_lower, tmp_data, ldm,
4336  anorm, rcond, pz, piz, tmp_err
4337  F77_CHAR_ARG_LEN (1)));
4338 
4339  err = tmp_err;
4340 
4341  if (err != 0)
4342  err = -2;
4343 
4344  volatile double rcond_plus_one = rcond + 1.0;
4345 
4346  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4347  {
4348  err = -2;
4349 
4350  if (sing_handler)
4351  {
4352  sing_handler (rcond);
4353  mattype.mark_as_rectangular ();
4354  }
4355  else
4357  }
4358  }
4359  else
4360  rcond = 1.;
4361 
4362  if (err == 0)
4363  {
4364  retval = b;
4365  double *result = retval.fortran_vec ();
4366 
4367  F77_INT b_nr = octave::to_f77_int (b.rows ());
4368  F77_INT b_nc = octave::to_f77_int (b.cols ());
4369 
4370  F77_XFCN (dpbtrs, DPBTRS,
4371  (F77_CONST_CHAR_ARG2 (&job, 1),
4372  tmp_nr, n_lower, b_nc, tmp_data,
4373  ldm, result, b_nr, tmp_err
4374  F77_CHAR_ARG_LEN (1)));
4375 
4376  err = tmp_err;
4377 
4378  if (err != 0)
4379  {
4380  // FIXME: Should this be a warning?
4381  (*current_liboctave_error_handler)
4382  ("SparseMatrix::solve solve failed");
4383  err = -1;
4384  }
4385  }
4386  }
4387  }
4388 
4389  if (typ == MatrixType::Banded)
4390  {
4391  // Create the storage for the banded form of the sparse matrix
4392  F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
4393  F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4394  F77_INT ldm = n_upper + 2 * n_lower + 1;
4395 
4396  Matrix m_band (ldm, nc);
4397  double *tmp_data = m_band.fortran_vec ();
4398 
4399  if (! mattype.is_dense ())
4400  {
4401  octave_idx_type ii = 0;
4402 
4403  for (F77_INT j = 0; j < ldm; j++)
4404  for (octave_idx_type i = 0; i < nc; i++)
4405  tmp_data[ii++] = 0.;
4406  }
4407 
4408  for (octave_idx_type j = 0; j < nc; j++)
4409  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4410  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4411 
4412  // Calculate the norm of the matrix, for later use.
4413  double anorm = 0.0;
4414  if (calc_cond)
4415  {
4416  for (octave_idx_type j = 0; j < nr; j++)
4417  {
4418  double atmp = 0.;
4419  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4420  atmp += fabs (data (i));
4421  if (atmp > anorm)
4422  anorm = atmp;
4423  }
4424  }
4425 
4426  F77_INT tmp_nr = octave::to_f77_int (nr);
4427 
4428  Array<F77_INT> ipvt (dim_vector (nr, 1));
4429  F77_INT *pipvt = ipvt.fortran_vec ();
4430 
4431  F77_INT tmp_err = 0;
4432 
4433  F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4434  tmp_data, ldm, pipvt, tmp_err));
4435 
4436  err = tmp_err;
4437 
4438  // Throw away extra info LAPACK gives so as to not
4439  // change output.
4440  if (err != 0)
4441  {
4442  err = -2;
4443  rcond = 0.0;
4444 
4445  if (sing_handler)
4446  {
4447  sing_handler (rcond);
4448  mattype.mark_as_rectangular ();
4449  }
4450  else
4452  }
4453  else
4454  {
4455  if (calc_cond)
4456  {
4457  char job = '1';
4458  Array<double> z (dim_vector (3 * nr, 1));
4459  double *pz = z.fortran_vec ();
4460  Array<F77_INT> iz (dim_vector (nr, 1));
4461  F77_INT *piz = iz.fortran_vec ();
4462 
4463  F77_INT tmp_nc = octave::to_f77_int (nc);
4464 
4465  F77_XFCN (dgbcon, DGBCON,
4466  (F77_CONST_CHAR_ARG2 (&job, 1),
4467  tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4468  anorm, rcond, pz, piz, tmp_err
4469  F77_CHAR_ARG_LEN (1)));
4470 
4471  err = tmp_err;
4472 
4473  if (err != 0)
4474  err = -2;
4475 
4476  volatile double rcond_plus_one = rcond + 1.0;
4477 
4478  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4479  {
4480  err = -2;
4481 
4482  if (sing_handler)
4483  {
4484  sing_handler (rcond);
4485  mattype.mark_as_rectangular ();
4486  }
4487  else
4489  }
4490  }
4491  else
4492  rcond = 1.;
4493 
4494  if (err == 0)
4495  {
4496  retval = b;
4497  double *result = retval.fortran_vec ();
4498 
4499  F77_INT b_nr = octave::to_f77_int (b.rows ());
4500  F77_INT b_nc = octave::to_f77_int (b.cols ());
4501 
4502  char job = 'N';
4503  F77_XFCN (dgbtrs, DGBTRS,
4504  (F77_CONST_CHAR_ARG2 (&job, 1),
4505  tmp_nr, n_lower, n_upper, b_nc, tmp_data,
4506  ldm, pipvt, result, b_nr, tmp_err
4507  F77_CHAR_ARG_LEN (1)));
4508 
4509  err = tmp_err;
4510  }
4511  }
4512  }
4513  else if (typ != MatrixType::Banded_Hermitian)
4514  (*current_liboctave_error_handler) ("incorrect matrix type");
4515  }
4516 
4517  return retval;
4518 }
4519 
4522  octave_idx_type& err, double& rcond,
4523  solve_singularity_handler sing_handler,
4524  bool calc_cond) const
4525 {
4527 
4528  octave_idx_type nr = rows ();
4529  octave_idx_type nc = cols ();
4530  err = 0;
4531 
4532  if (nr != nc || nr != b.rows ())
4534  ("matrix dimension mismatch solution of linear equations");
4535 
4536  if (nr == 0 || b.cols () == 0)
4537  retval = SparseMatrix (nc, b.cols ());
4538  else
4539  {
4540  // Print spparms("spumoni") info if requested
4541  volatile int typ = mattype.type ();
4542  mattype.info ();
4543 
4544  if (typ == MatrixType::Banded_Hermitian)
4545  {
4546  F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4547  F77_INT ldm = octave::to_f77_int (n_lower + 1);
4548 
4549  Matrix m_band (ldm, nc);
4550  double *tmp_data = m_band.fortran_vec ();
4551 
4552  if (! mattype.is_dense ())
4553  {
4554  octave_idx_type ii = 0;
4555 
4556  for (F77_INT j = 0; j < ldm; j++)
4557  for (octave_idx_type i = 0; i < nc; i++)
4558  tmp_data[ii++] = 0.;
4559  }
4560 
4561  for (octave_idx_type j = 0; j < nc; j++)
4562  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4563  {
4564  octave_idx_type ri = ridx (i);
4565  if (ri >= j)
4566  m_band(ri - j, j) = data (i);
4567  }
4568 
4569  // Calculate the norm of the matrix, for later use.
4570  double anorm;
4571  if (calc_cond)
4572  anorm = m_band.abs ().sum ().row (0).max ();
4573 
4574  F77_INT tmp_nr = octave::to_f77_int (nr);
4575 
4576  F77_INT tmp_err = 0;
4577 
4578  char job = 'L';
4579  F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4580  tmp_nr, n_lower, tmp_data, ldm, tmp_err
4581  F77_CHAR_ARG_LEN (1)));
4582 
4583  err = tmp_err;
4584 
4585  if (err != 0)
4586  {
4587  mattype.mark_as_unsymmetric ();
4588  typ = MatrixType::Banded;
4589  rcond = 0.0;
4590  err = 0;
4591  }
4592  else
4593  {
4594  if (calc_cond)
4595  {
4596  Array<double> z (dim_vector (3 * nr, 1));
4597  double *pz = z.fortran_vec ();
4598  Array<F77_INT> iz (dim_vector (nr, 1));
4599  F77_INT *piz = iz.fortran_vec ();
4600 
4601  F77_XFCN (dpbcon, DPBCON,
4602  (F77_CONST_CHAR_ARG2 (&job, 1),
4603  tmp_nr, n_lower, tmp_data, ldm,
4604  anorm, rcond, pz, piz, tmp_err
4605  F77_CHAR_ARG_LEN (1)));
4606 
4607  err = tmp_err;
4608 
4609  if (err != 0)
4610  err = -2;
4611 
4612  volatile double rcond_plus_one = rcond + 1.0;
4613 
4614  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4615  {
4616  err = -2;
4617 
4618  if (sing_handler)
4619  {
4620  sing_handler (rcond);
4621  mattype.mark_as_rectangular ();
4622  }
4623  else
4625  }
4626  }
4627  else
4628  rcond = 1.;
4629 
4630  if (err == 0)
4631  {
4632  F77_INT b_nr = octave::to_f77_int (b.rows ());
4633  octave_idx_type b_nc = b.cols ();
4634  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4635 
4636  // Take a first guess that the number of nonzero terms
4637  // will be as many as in b
4638  volatile octave_idx_type x_nz = b.nnz ();
4639  volatile octave_idx_type ii = 0;
4640  retval = SparseMatrix (b_nr, b_nc, x_nz);
4641 
4642  retval.xcidx (0) = 0;
4643  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4644  {
4645  for (F77_INT i = 0; i < b_nr; i++)
4646  Bx[i] = b.elem (i, j);
4647 
4648  F77_XFCN (dpbtrs, DPBTRS,
4649  (F77_CONST_CHAR_ARG2 (&job, 1),
4650  tmp_nr, n_lower, 1, tmp_data,
4651  ldm, Bx, b_nr, tmp_err
4652  F77_CHAR_ARG_LEN (1)));
4653 
4654  err = tmp_err;
4655 
4656  if (err != 0)
4657  {
4658  // FIXME: Should this be a warning?
4659  (*current_liboctave_error_handler)
4660  ("SparseMatrix::solve solve failed");
4661  err = -1;
4662  break;
4663  }
4664 
4665  for (F77_INT i = 0; i < b_nr; i++)
4666  {
4667  double tmp = Bx[i];
4668  if (tmp != 0.0)
4669  {
4670  if (ii == x_nz)
4671  {
4672  // Resize the sparse matrix
4674  sz = (static_cast<double> (b_nc) - j) / b_nc
4675  * x_nz;
4676  sz = x_nz + (sz > 100 ? sz : 100);
4677  retval.change_capacity (sz);
4678  x_nz = sz;
4679  }
4680  retval.xdata (ii) = tmp;
4681  retval.xridx (ii++) = i;
4682  }
4683  }
4684  retval.xcidx (j+1) = ii;
4685  }
4686 
4687  retval.maybe_compress ();
4688  }
4689  }
4690  }
4691 
4692  if (typ == MatrixType::Banded)
4693  {
4694  // Create the storage for the banded form of the sparse matrix
4695  F77_INT n_upper = octave::to_f77_int (mattype.nupper ());
4696  F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4697  F77_INT ldm = octave::to_f77_int (n_upper + 2 * n_lower + 1);
4698 
4699  Matrix m_band (ldm, nc);
4700  double *tmp_data = m_band.fortran_vec ();
4701 
4702  if (! mattype.is_dense ())
4703  {
4704  octave_idx_type ii = 0;
4705 
4706  for (octave_idx_type j = 0; j < ldm; j++)
4707  for (octave_idx_type i = 0; i < nc; i++)
4708  tmp_data[ii++] = 0.;
4709  }
4710 
4711  for (octave_idx_type j = 0; j < nc; j++)
4712  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4713  m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4714 
4715  // Calculate the norm of the matrix, for later use.
4716  double anorm;
4717  if (calc_cond)
4718  {
4719  for (octave_idx_type j = 0; j < nr; j++)
4720  {
4721  double atmp = 0.;
4722  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4723  atmp += fabs (data (i));
4724  if (atmp > anorm)
4725  anorm = atmp;
4726  }
4727  }
4728 
4729  F77_INT tmp_nr = octave::to_f77_int (nr);
4730 
4731  Array<F77_INT> ipvt (dim_vector (nr, 1));
4732  F77_INT *pipvt = ipvt.fortran_vec ();
4733 
4734  F77_INT tmp_err = 0;
4735 
4736  F77_XFCN (dgbtrf, DGBTRF, (tmp_nr, tmp_nr, n_lower, n_upper,
4737  tmp_data, ldm, pipvt, tmp_err));
4738 
4739  err = tmp_err;
4740 
4741  if (err != 0)
4742  {
4743  err = -2;
4744  rcond = 0.0;
4745 
4746  if (sing_handler)
4747  {
4748  sing_handler (rcond);
4749  mattype.mark_as_rectangular ();
4750  }
4751  else
4753  }
4754  else
4755  {
4756  if (calc_cond)
4757  {
4758  char job = '1';
4759  Array<double> z (dim_vector (3 * nr, 1));
4760  double *pz = z.fortran_vec ();
4761  Array<F77_INT> iz (dim_vector (nr, 1));
4762  F77_INT *piz = iz.fortran_vec ();
4763 
4764  F77_INT tmp_nc = octave::to_f77_int (nc);
4765 
4766  F77_XFCN (dgbcon, DGBCON,
4767  (F77_CONST_CHAR_ARG2 (&job, 1),
4768  tmp_nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4769  anorm, rcond, pz, piz, tmp_err
4770  F77_CHAR_ARG_LEN (1)));
4771 
4772  err = tmp_err;
4773 
4774  if (err != 0)
4775  err = -2;
4776 
4777  volatile double rcond_plus_one = rcond + 1.0;
4778 
4779  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4780  {
4781  err = -2;
4782 
4783  if (sing_handler)
4784  {
4785  sing_handler (rcond);
4786  mattype.mark_as_rectangular ();
4787  }
4788  else
4790  }
4791  }
4792  else
4793  rcond = 1.;
4794 
4795  if (err == 0)
4796  {
4797  char job = 'N';
4798  volatile octave_idx_type x_nz = b.nnz ();
4799  octave_idx_type b_nc = b.cols ();
4800  retval = SparseMatrix (nr, b_nc, x_nz);
4801  retval.xcidx (0) = 0;
4802  volatile octave_idx_type ii = 0;
4803 
4804  OCTAVE_LOCAL_BUFFER (double, work, nr);
4805 
4806  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4807  {
4808  for (octave_idx_type i = 0; i < nr; i++)
4809  work[i] = 0.;
4810  for (octave_idx_type i = b.cidx (j);
4811  i < b.cidx (j+1); i++)
4812  work[b.ridx (i)] = b.data (i);
4813 
4814  F77_INT b_nr = octave::to_f77_int (b.rows ());
4815 
4816  F77_XFCN (dgbtrs, DGBTRS,
4817  (F77_CONST_CHAR_ARG2 (&job, 1),
4818  tmp_nr, n_lower, n_upper, 1, tmp_data,
4819  ldm, pipvt, work, b_nr, tmp_err
4820  F77_CHAR_ARG_LEN (1)));
4821 
4822  err = tmp_err;
4823 
4824  // Count nonzeros in work vector and adjust
4825  // space in retval if needed
4826  octave_idx_type new_nnz = 0;
4827  for (octave_idx_type i = 0; i < nr; i++)
4828  if (work[i] != 0.)
4829  new_nnz++;
4830 
4831  if (ii + new_nnz > x_nz)
4832  {
4833  // Resize the sparse matrix
4834  octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4835  retval.change_capacity (sz);
4836  x_nz = sz;
4837  }
4838 
4839  for (octave_idx_type i = 0; i < nr; i++)
4840  if (work[i] != 0.)
4841  {
4842  retval.xridx (ii) = i;
4843  retval.xdata (ii++) = work[i];
4844  }
4845  retval.xcidx (j+1) = ii;
4846  }
4847 
4848  retval.maybe_compress ();
4849  }
4850  }
4851  }
4852  else if (typ != MatrixType::Banded_Hermitian)
4853  (*current_liboctave_error_handler) ("incorrect matrix type");
4854  }
4855 
4856  return retval;
4857 }
4858 
4861  octave_idx_type& err, double& rcond,
4862  solve_singularity_handler sing_handler,
4863  bool calc_cond) const
4864 {
4866 
4867  octave_idx_type nr = rows ();
4868  octave_idx_type nc = cols ();
4869  err = 0;
4870 
4871  if (nr != nc || nr != b.rows ())
4873  ("matrix dimension mismatch solution of linear equations");
4874 
4875  if (nr == 0 || b.cols () == 0)
4876  retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4877  else
4878  {
4879  // Print spparms("spumoni") info if requested
4880  volatile int typ = mattype.type ();
4881  mattype.info ();
4882 
4883  if (typ == MatrixType::Banded_Hermitian)
4884  {
4885  F77_INT n_lower = octave::to_f77_int (mattype.nlower ());
4886  F77_INT ldm = n_lower + 1;
4887 
4888  Matrix m_band (ldm, nc);
4889  double *tmp_data = m_band.fortran_vec ();
4890 
4891  if (! mattype.is_dense ())
4892  {
4893  octave_idx_type ii = 0;
4894 
4895  for (F77_INT j = 0; j < ldm; j++)
4896  for (octave_idx_type i = 0; i < nc; i++)
4897  tmp_data[ii++] = 0.;
4898  }
4899 
4900  for (octave_idx_type j = 0; j < nc; j++)
4901  for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4902  {
4903  octave_idx_type ri = ridx (i);
4904  if (ri >= j)
4905  m_band(ri - j, j) = data (i);
4906  }
4907 
4908  // Calculate the norm of the matrix, for later use.
4909  double anorm;
4910  if (calc_cond)
4911  anorm = m_band.abs ().sum ().row (0).max ();
4912 
4913  F77_INT tmp_nr = octave::to_f77_int (nr);
4914 
4915  F77_INT tmp_err = 0;
4916 
4917  char job = 'L';
4918  F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4919  tmp_nr, n_lower, tmp_data, ldm, tmp_err
4920  F77_CHAR_ARG_LEN (1)));
4921 
4922  err = tmp_err;
4923 
4924  if (err != 0)
4925  {
4926  // Matrix is not positive definite!! Fall through to
4927  // unsymmetric banded solver.
4928  mattype.mark_as_unsymmetric ();
4929  typ = MatrixType::Banded;
4930  rcond = 0.0;
4931  err = 0;
4932  }
4933  else
4934  {
4935  if (calc_cond)
4936  {
4937  Array<double> z (dim_vector (3 * nr, 1));
4938  double *pz = z.fortran_vec ();
4939  Array<F77_INT> iz (dim_vector (nr, 1));
4940  F77_INT *piz = iz.fortran_vec ();
4941 
4942  F77_XFCN (dpbcon, DPBCON,
4943  (F77_CONST_CHAR_ARG2 (&job, 1),
4944  tmp_nr, n_lower, tmp_data, ldm,
4945  anorm, rcond, pz, piz, tmp_err
4946  F77_CHAR_ARG_LEN (1)));
4947 
4948  err = tmp_err;
4949 
4950  if (err != 0)
4951  err = -2;
4952 
4953  volatile double rcond_plus_one = rcond + 1.0;
4954 
4955  if (rcond_plus_one == 1.0 || octave::math::isnan (rcond))
4956  {
4957  err = -2;
4958 
4959  if (sing_handler)
4960  {
4961  sing_handler (rcond);
4962  mattype.mark_as_rectangular ();
4963  }
4964  else
4966  }
4967  }
4968  else
4969  rcond = 1.;
4970 
4971  if (err == 0)
4972  {
4973  F77_INT b_nr = octave::to_f77_int (b.rows ());
4974  octave_idx_type b_nc = b.cols ();
4975 
4976  OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
4977  OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
4978 
4979  retval.resize (b_nr, b_nc);
4980 
4981  for (volatile octave_idx_type j = 0; j < b_nc; j++)
4982  {
4983  for (F77_INT i = 0; i < b_nr; i++)
4984  {
4985  Complex c = b(i,j);
4986  Bx[i] = c.real ();
4987  Bz[i] = c.imag ();
4988  }
4989 
4990  F77_XFCN (dpbtrs, DPBTRS,
4991  (F77_CONST_CHAR_ARG2 (&job, 1),
4992  tmp_nr, n_lower, 1, tmp_data,
4993  ldm, Bx, b_nr, tmp_err
4994  F77_CHAR_ARG_LEN (1)));
4995 
4996  err = tmp_err;
4997 
4998