3
Copyright (C) 2004-2013 David Bateman
4
Copyright (C) 1998-2004 Andy Adler
5
Copyright (C) 2010 VZLU Prague
7
This file is part of Octave.
9
Octave is free software; you can redistribute it and/or modify it
10
under the terms of the GNU General Public License as published by the
11
Free Software Foundation; either version 3 of the License, or (at your
12
option) any later version.
14
Octave is distributed in the hope that it will be useful, but WITHOUT
15
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19
You should have received a copy of the GNU General Public License
20
along with Octave; see the file COPYING. If not, see
21
<http://www.gnu.org/licenses/>.
36
#include "lo-mappers.h"
38
#include "dRowVector.h"
39
#include "oct-locbuf.h"
41
#include "dDiagMatrix.h"
42
#include "CDiagMatrix.h"
44
#include "boolSparse.h"
47
#include "oct-spparms.h"
48
#include "SparseCmplxLU.h"
49
#include "oct-sparse.h"
50
#include "sparse-util.h"
51
#include "SparseCmplxCHOL.h"
52
#include "SparseCmplxQR.h"
54
#include "Sparse-diag-op-defs.h"
56
#include "Sparse-perm-op-defs.h"
57
#include "mx-inlines.cc"
59
// Define whether to use a basic QR solver or one that uses a Dulmange
60
// Mendelsohn factorization to seperate the problem into under-determined,
61
// well-determined and over-determined parts and solves them seperately
63
#include "sparse-dmsolve.cc"
66
// Fortran functions we call.
70
F77_FUNC (zgbtrf, ZGBTRF) (const octave_idx_type&, const octave_idx_type&,
71
const octave_idx_type&, const octave_idx_type&,
72
Complex*, const octave_idx_type&,
73
octave_idx_type*, octave_idx_type&);
76
F77_FUNC (zgbtrs, ZGBTRS) (F77_CONST_CHAR_ARG_DECL,
77
const octave_idx_type&, const octave_idx_type&,
78
const octave_idx_type&, const octave_idx_type&,
79
const Complex*, const octave_idx_type&,
80
const octave_idx_type*, Complex*,
81
const octave_idx_type&, octave_idx_type&
82
F77_CHAR_ARG_LEN_DECL);
85
F77_FUNC (zgbcon, ZGBCON) (F77_CONST_CHAR_ARG_DECL,
86
const octave_idx_type&, const octave_idx_type&,
87
const octave_idx_type&, Complex*,
88
const octave_idx_type&, const octave_idx_type*,
89
const double&, double&, Complex*, double*,
91
F77_CHAR_ARG_LEN_DECL);
94
F77_FUNC (zpbtrf, ZPBTRF) (F77_CONST_CHAR_ARG_DECL,
95
const octave_idx_type&, const octave_idx_type&,
96
Complex*, const octave_idx_type&, octave_idx_type&
97
F77_CHAR_ARG_LEN_DECL);
100
F77_FUNC (zpbtrs, ZPBTRS) (F77_CONST_CHAR_ARG_DECL,
101
const octave_idx_type&, const octave_idx_type&,
102
const octave_idx_type&, Complex*,
103
const octave_idx_type&, Complex*,
104
const octave_idx_type&, octave_idx_type&
105
F77_CHAR_ARG_LEN_DECL);
108
F77_FUNC (zpbcon, ZPBCON) (F77_CONST_CHAR_ARG_DECL,
109
const octave_idx_type&, const octave_idx_type&,
110
Complex*, const octave_idx_type&, const double&,
111
double&, Complex*, double*, octave_idx_type&
112
F77_CHAR_ARG_LEN_DECL);
115
F77_FUNC (zgttrf, ZGTTRF) (const octave_idx_type&, Complex*, Complex*,
116
Complex*, Complex*, octave_idx_type*,
120
F77_FUNC (zgttrs, ZGTTRS) (F77_CONST_CHAR_ARG_DECL,
121
const octave_idx_type&, const octave_idx_type&,
122
const Complex*, const Complex*, const Complex*,
123
const Complex*, const octave_idx_type*,
124
Complex *, const octave_idx_type&,
126
F77_CHAR_ARG_LEN_DECL);
129
F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&,
130
double*, Complex*, Complex*,
131
const octave_idx_type&, octave_idx_type&);
134
F77_FUNC (zgtsv, ZGTSV) (const octave_idx_type&, const octave_idx_type&,
135
Complex*, Complex*, Complex*, Complex*,
136
const octave_idx_type&, octave_idx_type&);
139
SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a)
140
: MSparse<Complex> (a)
144
SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a)
145
: MSparse<Complex> (a.rows (), a.cols (), a.nnz ())
147
octave_idx_type nc = cols ();
148
octave_idx_type nz = a.nnz ();
150
for (octave_idx_type i = 0; i < nc + 1; i++)
151
cidx (i) = a.cidx (i);
153
for (octave_idx_type i = 0; i < nz; i++)
155
data (i) = Complex (a.data (i));
156
ridx (i) = a.ridx (i);
160
SparseComplexMatrix::SparseComplexMatrix (const ComplexDiagMatrix& a)
161
: MSparse<Complex> (a.rows (), a.cols (), a.length ())
163
octave_idx_type j = 0, l = a.length ();
164
for (octave_idx_type i = 0; i < l; i++)
174
for (octave_idx_type i = l; i <= a.cols (); i++)
178
SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const
180
octave_idx_type nr = rows ();
181
octave_idx_type nc = cols ();
182
octave_idx_type nz = nnz ();
183
octave_idx_type nr_a = a.rows ();
184
octave_idx_type nc_a = a.cols ();
185
octave_idx_type nz_a = a.nnz ();
187
if (nr != nr_a || nc != nc_a || nz != nz_a)
190
for (octave_idx_type i = 0; i < nc + 1; i++)
191
if (cidx (i) != a.cidx (i))
194
for (octave_idx_type i = 0; i < nz; i++)
195
if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
202
SparseComplexMatrix::operator != (const SparseComplexMatrix& a) const
204
return !(*this == a);
208
SparseComplexMatrix::is_hermitian (void) const
210
octave_idx_type nr = rows ();
211
octave_idx_type nc = cols ();
213
if (nr == nc && nr > 0)
215
for (octave_idx_type j = 0; j < nc; j++)
217
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
219
octave_idx_type ri = ridx (i);
225
for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++)
229
if (data (i) == conj (data (k)))
247
static const Complex Complex_NaN_result (octave_NaN, octave_NaN);
250
SparseComplexMatrix::max (int dim) const
252
Array<octave_idx_type> dummy_idx;
253
return max (dummy_idx, dim);
257
SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const
259
SparseComplexMatrix result;
260
dim_vector dv = dims ();
261
octave_idx_type nr = dv(0);
262
octave_idx_type nc = dv(1);
265
if (dim >= dv.length ())
267
idx_arg.resize (dim_vector (nr, nc), 0);
272
dim = dv.first_non_singleton ();
276
idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
278
if (nr == 0 || nc == 0 || dim >= dv.length ())
279
return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
281
octave_idx_type nel = 0;
282
for (octave_idx_type j = 0; j < nc; j++)
285
double abs_max = octave_NaN;
286
octave_idx_type idx_j = 0;
287
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
289
if (ridx (i) != idx_j)
301
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
303
Complex tmp = data (i);
308
double abs_tmp = std::abs (tmp);
310
if (xisnan (abs_max) || abs_tmp > abs_max)
318
idx_arg.elem (j) = xisnan (tmp_max) ? 0 : idx_j;
323
result = SparseComplexMatrix (1, nc, nel);
325
octave_idx_type ii = 0;
326
result.xcidx (0) = 0;
327
for (octave_idx_type j = 0; j < nc; j++)
329
Complex tmp = elem (idx_arg(j), j);
332
result.xdata (ii) = tmp;
333
result.xridx (ii++) = 0;
335
result.xcidx (j+1) = ii;
340
idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
342
if (nr == 0 || nc == 0 || dim >= dv.length ())
343
return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
345
OCTAVE_LOCAL_BUFFER (octave_idx_type, found, nr);
347
for (octave_idx_type i = 0; i < nr; i++)
350
for (octave_idx_type j = 0; j < nc; j++)
351
for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
352
if (found[ridx (i)] == -j)
353
found[ridx (i)] = -j - 1;
355
for (octave_idx_type i = 0; i < nr; i++)
356
if (found[i] > -nc && found[i] < 0)
357
idx_arg.elem (i) = -found[i];
359
for (octave_idx_type j = 0; j < nc; j++)
361
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
363
octave_idx_type ir = ridx (i);
364
octave_idx_type ix = idx_arg.elem (ir);
365
Complex tmp = data (i);
369
else if (ix == -1 || std::abs (tmp) > std::abs (elem (ir, ix)))
370
idx_arg.elem (ir) = j;
374
octave_idx_type nel = 0;
375
for (octave_idx_type j = 0; j < nr; j++)
376
if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
379
result = SparseComplexMatrix (nr, 1, nel);
381
octave_idx_type ii = 0;
382
result.xcidx (0) = 0;
383
result.xcidx (1) = nel;
384
for (octave_idx_type j = 0; j < nr; j++)
386
if (idx_arg(j) == -1)
389
result.xdata (ii) = Complex_NaN_result;
390
result.xridx (ii++) = j;
394
Complex tmp = elem (j, idx_arg(j));
397
result.xdata (ii) = tmp;
398
result.xridx (ii++) = j;
408
SparseComplexMatrix::min (int dim) const
410
Array<octave_idx_type> dummy_idx;
411
return min (dummy_idx, dim);
415
SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const
417
SparseComplexMatrix result;
418
dim_vector dv = dims ();
419
octave_idx_type nr = dv(0);
420
octave_idx_type nc = dv(1);
422
if (dim >= dv.length ())
424
idx_arg.resize (dim_vector (nr, nc), 0);
429
dim = dv.first_non_singleton ();
433
idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0);
435
if (nr == 0 || nc == 0 || dim >= dv.length ())
436
return SparseComplexMatrix (nr == 0 ? 0 : 1, nc);
438
octave_idx_type nel = 0;
439
for (octave_idx_type j = 0; j < nc; j++)
442
double abs_min = octave_NaN;
443
octave_idx_type idx_j = 0;
444
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
446
if (ridx (i) != idx_j)
458
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
460
Complex tmp = data (i);
465
double abs_tmp = std::abs (tmp);
467
if (xisnan (abs_min) || abs_tmp < abs_min)
475
idx_arg.elem (j) = xisnan (tmp_min) ? 0 : idx_j;
480
result = SparseComplexMatrix (1, nc, nel);
482
octave_idx_type ii = 0;
483
result.xcidx (0) = 0;
484
for (octave_idx_type j = 0; j < nc; j++)
486
Complex tmp = elem (idx_arg(j), j);
489
result.xdata (ii) = tmp;
490
result.xridx (ii++) = 0;
492
result.xcidx (j+1) = ii;
497
idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0);
499
if (nr == 0 || nc == 0 || dim >= dv.length ())
500
return SparseComplexMatrix (nr, nc == 0 ? 0 : 1);
502
OCTAVE_LOCAL_BUFFER (octave_idx_type, found, nr);
504
for (octave_idx_type i = 0; i < nr; i++)
507
for (octave_idx_type j = 0; j < nc; j++)
508
for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
509
if (found[ridx (i)] == -j)
510
found[ridx (i)] = -j - 1;
512
for (octave_idx_type i = 0; i < nr; i++)
513
if (found[i] > -nc && found[i] < 0)
514
idx_arg.elem (i) = -found[i];
516
for (octave_idx_type j = 0; j < nc; j++)
518
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
520
octave_idx_type ir = ridx (i);
521
octave_idx_type ix = idx_arg.elem (ir);
522
Complex tmp = data (i);
526
else if (ix == -1 || std::abs (tmp) < std::abs (elem (ir, ix)))
527
idx_arg.elem (ir) = j;
531
octave_idx_type nel = 0;
532
for (octave_idx_type j = 0; j < nr; j++)
533
if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.)
536
result = SparseComplexMatrix (nr, 1, nel);
538
octave_idx_type ii = 0;
539
result.xcidx (0) = 0;
540
result.xcidx (1) = nel;
541
for (octave_idx_type j = 0; j < nr; j++)
543
if (idx_arg(j) == -1)
546
result.xdata (ii) = Complex_NaN_result;
547
result.xridx (ii++) = j;
551
Complex tmp = elem (j, idx_arg(j));
554
result.xdata (ii) = tmp;
555
result.xridx (ii++) = j;
566
%!assert (max (max (speye (65536) * 1i)), sparse (1i))
567
%!assert (min (min (speye (65536) * 1i)), sparse (0))
568
%!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
569
%!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
570
%!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
571
%!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
572
%!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
573
%!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
574
%!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
575
%!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
580
SparseComplexMatrix::row (octave_idx_type i) const
582
octave_idx_type nc = columns ();
583
ComplexRowVector retval (nc, 0);
585
for (octave_idx_type j = 0; j < nc; j++)
586
for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
590
retval(j) = data (k);
599
SparseComplexMatrix::column (octave_idx_type i) const
601
octave_idx_type nr = rows ();
602
ComplexColumnVector retval (nr, 0);
604
for (octave_idx_type k = cidx (i); k < cidx (i+1); k++)
605
retval(ridx (k)) = data (k);
610
// destructive insert/delete/reorder operations
613
SparseComplexMatrix::insert (const SparseMatrix& a,
614
octave_idx_type r, octave_idx_type c)
616
SparseComplexMatrix tmp (a);
617
return insert (tmp /*a*/, r, c);
621
SparseComplexMatrix::insert (const SparseComplexMatrix& a,
622
octave_idx_type r, octave_idx_type c)
624
MSparse<Complex>::insert (a, r, c);
629
SparseComplexMatrix::insert (const SparseMatrix& a,
630
const Array<octave_idx_type>& indx)
632
SparseComplexMatrix tmp (a);
633
return insert (tmp /*a*/, indx);
637
SparseComplexMatrix::insert (const SparseComplexMatrix& a,
638
const Array<octave_idx_type>& indx)
640
MSparse<Complex>::insert (a, indx);
645
SparseComplexMatrix::concat (const SparseComplexMatrix& rb,
646
const Array<octave_idx_type>& ra_idx)
648
// Don't use numel to avoid all possiblity of an overflow
649
if (rb.rows () > 0 && rb.cols () > 0)
650
insert (rb, ra_idx(0), ra_idx(1));
655
SparseComplexMatrix::concat (const SparseMatrix& rb,
656
const Array<octave_idx_type>& ra_idx)
658
SparseComplexMatrix tmp (rb);
659
if (rb.rows () > 0 && rb.cols () > 0)
660
insert (tmp, ra_idx(0), ra_idx(1));
665
SparseComplexMatrix::matrix_value (void) const
667
return Sparse<Complex>::array_value ();
671
SparseComplexMatrix::hermitian (void) const
673
octave_idx_type nr = rows ();
674
octave_idx_type nc = cols ();
675
octave_idx_type nz = nnz ();
676
SparseComplexMatrix retval (nc, nr, nz);
678
for (octave_idx_type i = 0; i < nz; i++)
679
retval.xcidx (ridx (i) + 1)++;
680
// retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1)
682
for (octave_idx_type i = 1; i <= nr; i++)
684
const octave_idx_type tmp = retval.xcidx (i);
685
retval.xcidx (i) = nz;
688
// retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1)
690
for (octave_idx_type j = 0; j < nc; j++)
691
for (octave_idx_type k = cidx (j); k < cidx (j+1); k++)
693
octave_idx_type q = retval.xcidx (ridx (k) + 1)++;
694
retval.xridx (q) = j;
695
retval.xdata (q) = conj (data (k));
697
assert (nnz () == retval.xcidx (nr));
698
// retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1)
699
// and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets
705
conj (const SparseComplexMatrix& a)
707
octave_idx_type nr = a.rows ();
708
octave_idx_type nc = a.cols ();
709
octave_idx_type nz = a.nnz ();
710
SparseComplexMatrix retval (nc, nr, nz);
712
for (octave_idx_type i = 0; i < nc + 1; i++)
713
retval.cidx (i) = a.cidx (i);
715
for (octave_idx_type i = 0; i < nz; i++)
717
retval.data (i) = conj (a.data (i));
718
retval.ridx (i) = a.ridx (i);
725
SparseComplexMatrix::inverse (void) const
727
octave_idx_type info;
729
MatrixType mattype (*this);
730
return inverse (mattype, info, rcond, 0, 0);
734
SparseComplexMatrix::inverse (MatrixType& mattype) const
736
octave_idx_type info;
738
return inverse (mattype, info, rcond, 0, 0);
742
SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const
745
return inverse (mattype, info, rcond, 0, 0);
749
SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info,
750
double& rcond, const bool,
751
const bool calccond) const
753
SparseComplexMatrix retval;
755
octave_idx_type nr = rows ();
756
octave_idx_type nc = cols ();
759
if (nr == 0 || nc == 0 || nr != nc)
760
(*current_liboctave_error_handler) ("inverse requires square matrix");
763
// Print spparms("spumoni") info if requested
764
int typ = mattyp.type ();
767
if (typ == MatrixType::Diagonal ||
768
typ == MatrixType::Permuted_Diagonal)
770
if (typ == MatrixType::Permuted_Diagonal)
771
retval = transpose ();
775
// Force make_unique to be called
776
Complex *v = retval.data ();
780
double dmax = 0., dmin = octave_Inf;
781
for (octave_idx_type i = 0; i < nr; i++)
783
double tmp = std::abs (v[i]);
792
for (octave_idx_type i = 0; i < nr; i++)
796
(*current_liboctave_error_handler) ("incorrect matrix type");
803
SparseComplexMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info,
804
double& rcond, const bool,
805
const bool calccond) const
807
SparseComplexMatrix retval;
809
octave_idx_type nr = rows ();
810
octave_idx_type nc = cols ();
813
if (nr == 0 || nc == 0 || nr != nc)
814
(*current_liboctave_error_handler) ("inverse requires square matrix");
817
// Print spparms("spumoni") info if requested
818
int typ = mattyp.type ();
821
if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper ||
822
typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
825
double ainvnorm = 0.;
829
// Calculate the 1-norm of matrix for rcond calculation
830
for (octave_idx_type j = 0; j < nr; j++)
833
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
834
atmp += std::abs (data (i));
840
if (typ == MatrixType::Upper || typ == MatrixType::Lower)
842
octave_idx_type nz = nnz ();
843
octave_idx_type cx = 0;
844
octave_idx_type nz2 = nz;
845
retval = SparseComplexMatrix (nr, nc, nz2);
847
for (octave_idx_type i = 0; i < nr; i++)
850
// place the 1 in the identity position
851
octave_idx_type cx_colstart = cx;
856
retval.change_capacity (nz2);
859
retval.xcidx (i) = cx;
860
retval.xridx (cx) = i;
861
retval.xdata (cx) = 1.0;
864
// iterate accross columns of input matrix
865
for (octave_idx_type j = i+1; j < nr; j++)
868
// iterate to calculate sum
869
octave_idx_type colXp = retval.xcidx (i);
870
octave_idx_type colUp = cidx (j);
871
octave_idx_type rpX, rpU;
873
if (cidx (j) == cidx (j+1))
875
(*current_liboctave_error_handler)
876
("division by zero");
877
goto inverse_singular;
883
rpX = retval.xridx (colXp);
892
v -= retval.xdata (colXp) * data (colUp);
896
} while ((rpX<j) && (rpU<j) &&
897
(colXp<cx) && (colUp<nz));
901
if (typ == MatrixType::Upper)
902
colUp = cidx (j+1) - 1;
905
Complex pivot = data (colUp);
906
if (pivot == 0. || ridx (colUp) != j)
908
(*current_liboctave_error_handler)
909
("division by zero");
910
goto inverse_singular;
918
retval.change_capacity (nz2);
921
retval.xridx (cx) = j;
922
retval.xdata (cx) = v / pivot;
928
octave_idx_type colUp;
929
if (typ == MatrixType::Upper)
930
colUp = cidx (i+1) - 1;
933
Complex pivot = data (colUp);
934
if (pivot == 0. || ridx (colUp) != i)
936
(*current_liboctave_error_handler) ("division by zero");
937
goto inverse_singular;
941
for (octave_idx_type j = cx_colstart; j < cx; j++)
942
retval.xdata (j) /= pivot;
944
retval.xcidx (nr) = cx;
945
retval.maybe_compress ();
949
octave_idx_type nz = nnz ();
950
octave_idx_type cx = 0;
951
octave_idx_type nz2 = nz;
952
retval = SparseComplexMatrix (nr, nc, nz2);
954
OCTAVE_LOCAL_BUFFER (Complex, work, nr);
955
OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
957
octave_idx_type *perm = mattyp.triangular_perm ();
958
if (typ == MatrixType::Permuted_Upper)
960
for (octave_idx_type i = 0; i < nr; i++)
965
for (octave_idx_type i = 0; i < nr; i++)
967
for (octave_idx_type i = 0; i < nr; i++)
971
for (octave_idx_type i = 0; i < nr; i++)
974
octave_idx_type iidx = rperm[i];
976
for (octave_idx_type j = 0; j < nr; j++)
979
// place the 1 in the identity position
982
// iterate accross columns of input matrix
983
for (octave_idx_type j = iidx+1; j < nr; j++)
986
octave_idx_type jidx = perm[j];
987
// iterate to calculate sum
988
for (octave_idx_type k = cidx (jidx);
989
k < cidx (jidx+1); k++)
992
v -= work[ridx (k)] * data (k);
997
if (typ == MatrixType::Permuted_Upper)
998
pivot = data (cidx (jidx+1) - 1);
1000
pivot = data (cidx (jidx));
1003
(*current_liboctave_error_handler)
1004
("division by zero");
1005
goto inverse_singular;
1008
work[j] = v / pivot;
1012
octave_idx_type colUp;
1013
if (typ == MatrixType::Permuted_Upper)
1014
colUp = cidx (perm[iidx]+1) - 1;
1016
colUp = cidx (perm[iidx]);
1018
Complex pivot = data (colUp);
1021
(*current_liboctave_error_handler)
1022
("division by zero");
1023
goto inverse_singular;
1026
octave_idx_type new_cx = cx;
1027
for (octave_idx_type j = iidx; j < nr; j++)
1037
nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1038
retval.change_capacity (nz2);
1041
retval.xcidx (i) = cx;
1042
for (octave_idx_type j = iidx; j < nr; j++)
1045
retval.xridx (cx) = j;
1046
retval.xdata (cx++) = work[j];
1050
retval.xcidx (nr) = cx;
1051
retval.maybe_compress ();
1056
// Calculate the 1-norm of inverse matrix for rcond calculation
1057
for (octave_idx_type j = 0; j < nr; j++)
1060
for (octave_idx_type i = retval.cidx (j);
1061
i < retval.cidx (j+1); i++)
1062
atmp += std::abs (retval.data (i));
1063
if (atmp > ainvnorm)
1067
rcond = 1. / ainvnorm / anorm;
1071
(*current_liboctave_error_handler) ("incorrect matrix type");
1077
return SparseComplexMatrix ();
1081
SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info,
1082
double& rcond, int, int calc_cond) const
1084
int typ = mattype.type (false);
1085
SparseComplexMatrix ret;
1087
if (typ == MatrixType::Unknown)
1088
typ = mattype.type (*this);
1090
if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
1091
ret = dinverse (mattype, info, rcond, true, calc_cond);
1092
else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1093
ret = tinverse (mattype, info, rcond, true, calc_cond).transpose ();
1094
else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1096
MatrixType newtype = mattype.transpose ();
1097
ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond);
1101
if (mattype.is_hermitian ())
1103
MatrixType tmp_typ (MatrixType::Upper);
1104
SparseComplexCHOL fact (*this, info, false);
1105
rcond = fact.rcond ();
1109
SparseMatrix Q = fact.Q ();
1110
SparseComplexMatrix InvL = fact.L ().transpose ().
1111
tinverse (tmp_typ, info, rcond2,
1113
ret = Q * InvL.hermitian () * InvL * Q.transpose ();
1117
// Matrix is either singular or not positive definite
1118
mattype.mark_as_unsymmetric ();
1119
typ = MatrixType::Full;
1123
if (!mattype.is_hermitian ())
1125
octave_idx_type n = rows ();
1126
ColumnVector Qinit(n);
1127
for (octave_idx_type i = 0; i < n; i++)
1130
MatrixType tmp_typ (MatrixType::Upper);
1131
SparseComplexLU fact (*this, Qinit, Matrix (), false, false);
1132
rcond = fact.rcond ();
1134
SparseComplexMatrix InvL = fact.L ().transpose ().
1135
tinverse (tmp_typ, info, rcond2,
1137
SparseComplexMatrix InvU = fact.U ().
1138
tinverse (tmp_typ, info, rcond2,
1139
true, false).transpose ();
1140
ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1148
SparseComplexMatrix::determinant (void) const
1150
octave_idx_type info;
1152
return determinant (info, rcond, 0);
1156
SparseComplexMatrix::determinant (octave_idx_type& info) const
1159
return determinant (info, rcond, 0);
1163
SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond,
1169
octave_idx_type nr = rows ();
1170
octave_idx_type nc = cols ();
1172
if (nr == 0 || nc == 0 || nr != nc)
1174
retval = ComplexDET (1.0);
1180
// Setup the control parameters
1181
Matrix Control (UMFPACK_CONTROL, 1);
1182
double *control = Control.fortran_vec ();
1183
UMFPACK_ZNAME (defaults) (control);
1185
double tmp = octave_sparse_params::get_key ("spumoni");
1187
Control (UMFPACK_PRL) = tmp;
1189
tmp = octave_sparse_params::get_key ("piv_tol");
1192
Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
1193
Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
1196
// Set whether we are allowed to modify Q or not
1197
tmp = octave_sparse_params::get_key ("autoamd");
1199
Control (UMFPACK_FIXQ) = tmp;
1201
// Turn-off UMFPACK scaling for LU
1202
Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
1204
UMFPACK_ZNAME (report_control) (control);
1206
const octave_idx_type *Ap = cidx ();
1207
const octave_idx_type *Ai = ridx ();
1208
const Complex *Ax = data ();
1210
UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
1211
reinterpret_cast<const double *> (Ax),
1215
Matrix Info (1, UMFPACK_INFO);
1216
double *info = Info.fortran_vec ();
1217
int status = UMFPACK_ZNAME (qsymbolic)
1218
(nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0,
1219
0, &Symbolic, control, info);
1223
(*current_liboctave_error_handler)
1224
("SparseComplexMatrix::determinant symbolic factorization failed");
1226
UMFPACK_ZNAME (report_status) (control, status);
1227
UMFPACK_ZNAME (report_info) (control, info);
1229
UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1233
UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
1237
= UMFPACK_ZNAME (numeric) (Ap, Ai,
1238
reinterpret_cast<const double *> (Ax),
1239
0, Symbolic, &Numeric, control, info);
1240
UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1242
rcond = Info (UMFPACK_RCOND);
1246
(*current_liboctave_error_handler)
1247
("SparseComplexMatrix::determinant numeric factorization failed");
1249
UMFPACK_ZNAME (report_status) (control, status);
1250
UMFPACK_ZNAME (report_info) (control, info);
1252
UMFPACK_ZNAME (free_numeric) (&Numeric);
1256
UMFPACK_ZNAME (report_numeric) (Numeric, control);
1260
status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10,
1265
(*current_liboctave_error_handler)
1266
("SparseComplexMatrix::determinant error calculating determinant");
1268
UMFPACK_ZNAME (report_status) (control, status);
1269
UMFPACK_ZNAME (report_info) (control, info);
1272
retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10);
1274
UMFPACK_ZNAME (free_numeric) (&Numeric);
1279
(*current_liboctave_error_handler) ("UMFPACK not installed");
1286
SparseComplexMatrix::dsolve (MatrixType &mattype, const Matrix& b,
1287
octave_idx_type& err, double& rcond,
1288
solve_singularity_handler, bool calc_cond) const
1290
ComplexMatrix retval;
1292
octave_idx_type nr = rows ();
1293
octave_idx_type nc = cols ();
1294
octave_idx_type nm = (nc < nr ? nc : nr);
1297
if (nr != b.rows ())
1298
(*current_liboctave_error_handler)
1299
("matrix dimension mismatch solution of linear equations");
1300
else if (nr == 0 || nc == 0 || b.cols () == 0)
1301
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1304
// Print spparms("spumoni") info if requested
1305
int typ = mattype.type ();
1308
if (typ == MatrixType::Diagonal ||
1309
typ == MatrixType::Permuted_Diagonal)
1311
retval.resize (nc, b.cols (), Complex (0.,0.));
1312
if (typ == MatrixType::Diagonal)
1313
for (octave_idx_type j = 0; j < b.cols (); j++)
1314
for (octave_idx_type i = 0; i < nm; i++)
1315
retval(i,j) = b(i,j) / data (i);
1317
for (octave_idx_type j = 0; j < b.cols (); j++)
1318
for (octave_idx_type k = 0; k < nc; k++)
1319
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1320
retval(k,j) = b(ridx (i),j) / data (i);
1324
double dmax = 0., dmin = octave_Inf;
1325
for (octave_idx_type i = 0; i < nm; i++)
1327
double tmp = std::abs (data (i));
1333
rcond = dmin / dmax;
1339
(*current_liboctave_error_handler) ("incorrect matrix type");
1346
SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b,
1347
octave_idx_type& err, double& rcond,
1348
solve_singularity_handler,
1349
bool calc_cond) const
1351
SparseComplexMatrix retval;
1353
octave_idx_type nr = rows ();
1354
octave_idx_type nc = cols ();
1355
octave_idx_type nm = (nc < nr ? nc : nr);
1358
if (nr != b.rows ())
1359
(*current_liboctave_error_handler)
1360
("matrix dimension mismatch solution of linear equations");
1361
else if (nr == 0 || nc == 0 || b.cols () == 0)
1362
retval = SparseComplexMatrix (nc, b.cols ());
1365
// Print spparms("spumoni") info if requested
1366
int typ = mattype.type ();
1369
if (typ == MatrixType::Diagonal ||
1370
typ == MatrixType::Permuted_Diagonal)
1372
octave_idx_type b_nc = b.cols ();
1373
octave_idx_type b_nz = b.nnz ();
1374
retval = SparseComplexMatrix (nc, b_nc, b_nz);
1376
retval.xcidx (0) = 0;
1377
octave_idx_type ii = 0;
1378
if (typ == MatrixType::Diagonal)
1379
for (octave_idx_type j = 0; j < b.cols (); j++)
1381
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1383
if (b.ridx (i) >= nm)
1385
retval.xridx (ii) = b.ridx (i);
1386
retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1388
retval.xcidx (j+1) = ii;
1391
for (octave_idx_type j = 0; j < b.cols (); j++)
1393
for (octave_idx_type l = 0; l < nc; l++)
1394
for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1398
for (k = b.cidx (j); k < b.cidx (j+1); k++)
1399
if (ridx (i) == b.ridx (k))
1406
retval.xridx (ii) = l;
1407
retval.xdata (ii++) = b.data (k) / data (i);
1410
retval.xcidx (j+1) = ii;
1415
double dmax = 0., dmin = octave_Inf;
1416
for (octave_idx_type i = 0; i < nm; i++)
1418
double tmp = std::abs (data (i));
1424
rcond = dmin / dmax;
1430
(*current_liboctave_error_handler) ("incorrect matrix type");
1437
SparseComplexMatrix::dsolve (MatrixType &mattype, const ComplexMatrix& b,
1438
octave_idx_type& err, double& rcond,
1439
solve_singularity_handler,
1440
bool calc_cond) const
1442
ComplexMatrix retval;
1444
octave_idx_type nr = rows ();
1445
octave_idx_type nc = cols ();
1446
octave_idx_type nm = (nc < nr ? nc : nr);
1449
if (nr != b.rows ())
1450
(*current_liboctave_error_handler)
1451
("matrix dimension mismatch solution of linear equations");
1452
else if (nr == 0 || nc == 0 || b.cols () == 0)
1453
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1456
// Print spparms("spumoni") info if requested
1457
int typ = mattype.type ();
1460
if (typ == MatrixType::Diagonal ||
1461
typ == MatrixType::Permuted_Diagonal)
1463
retval.resize (nc, b.cols (), Complex (0.,0.));
1464
if (typ == MatrixType::Diagonal)
1465
for (octave_idx_type j = 0; j < b.cols (); j++)
1466
for (octave_idx_type i = 0; i < nm; i++)
1467
retval(i,j) = b(i,j) / data (i);
1469
for (octave_idx_type j = 0; j < b.cols (); j++)
1470
for (octave_idx_type k = 0; k < nc; k++)
1471
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1472
retval(k,j) = b(ridx (i),j) / data (i);
1476
double dmax = 0., dmin = octave_Inf;
1477
for (octave_idx_type i = 0; i < nr; i++)
1479
double tmp = std::abs (data (i));
1485
rcond = dmin / dmax;
1491
(*current_liboctave_error_handler) ("incorrect matrix type");
1498
SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b,
1499
octave_idx_type& err, double& rcond,
1500
solve_singularity_handler,
1501
bool calc_cond) const
1503
SparseComplexMatrix retval;
1505
octave_idx_type nr = rows ();
1506
octave_idx_type nc = cols ();
1507
octave_idx_type nm = (nc < nr ? nc : nr);
1510
if (nr != b.rows ())
1511
(*current_liboctave_error_handler)
1512
("matrix dimension mismatch solution of linear equations");
1513
else if (nr == 0 || nc == 0 || b.cols () == 0)
1514
retval = SparseComplexMatrix (nc, b.cols ());
1517
// Print spparms("spumoni") info if requested
1518
int typ = mattype.type ();
1521
if (typ == MatrixType::Diagonal ||
1522
typ == MatrixType::Permuted_Diagonal)
1524
octave_idx_type b_nc = b.cols ();
1525
octave_idx_type b_nz = b.nnz ();
1526
retval = SparseComplexMatrix (nc, b_nc, b_nz);
1528
retval.xcidx (0) = 0;
1529
octave_idx_type ii = 0;
1530
if (typ == MatrixType::Diagonal)
1531
for (octave_idx_type j = 0; j < b.cols (); j++)
1533
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1535
if (b.ridx (i) >= nm)
1537
retval.xridx (ii) = b.ridx (i);
1538
retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1540
retval.xcidx (j+1) = ii;
1543
for (octave_idx_type j = 0; j < b.cols (); j++)
1545
for (octave_idx_type l = 0; l < nc; l++)
1546
for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1550
for (k = b.cidx (j); k < b.cidx (j+1); k++)
1551
if (ridx (i) == b.ridx (k))
1558
retval.xridx (ii) = l;
1559
retval.xdata (ii++) = b.data (k) / data (i);
1562
retval.xcidx (j+1) = ii;
1567
double dmax = 0., dmin = octave_Inf;
1568
for (octave_idx_type i = 0; i < nm; i++)
1570
double tmp = std::abs (data (i));
1576
rcond = dmin / dmax;
1582
(*current_liboctave_error_handler) ("incorrect matrix type");
1589
SparseComplexMatrix::utsolve (MatrixType &mattype, const Matrix& b,
1590
octave_idx_type& err, double& rcond,
1591
solve_singularity_handler sing_handler,
1592
bool calc_cond) const
1594
ComplexMatrix retval;
1596
octave_idx_type nr = rows ();
1597
octave_idx_type nc = cols ();
1598
octave_idx_type nm = (nc > nr ? nc : nr);
1601
if (nr != b.rows ())
1602
(*current_liboctave_error_handler)
1603
("matrix dimension mismatch solution of linear equations");
1604
else if (nr == 0 || nc == 0 || b.cols () == 0)
1605
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1608
// Print spparms("spumoni") info if requested
1609
int typ = mattype.type ();
1612
if (typ == MatrixType::Permuted_Upper ||
1613
typ == MatrixType::Upper)
1616
double ainvnorm = 0.;
1617
octave_idx_type b_nc = b.cols ();
1622
// Calculate the 1-norm of matrix for rcond calculation
1623
for (octave_idx_type j = 0; j < nc; j++)
1626
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1627
atmp += std::abs (data (i));
1633
if (typ == MatrixType::Permuted_Upper)
1635
retval.resize (nc, b_nc);
1636
octave_idx_type *perm = mattype.triangular_perm ();
1637
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1639
for (octave_idx_type j = 0; j < b_nc; j++)
1641
for (octave_idx_type i = 0; i < nr; i++)
1643
for (octave_idx_type i = nr; i < nc; i++)
1646
for (octave_idx_type k = nc-1; k >= 0; k--)
1648
octave_idx_type kidx = perm[k];
1652
if (ridx (cidx (kidx+1)-1) != k ||
1653
data (cidx (kidx+1)-1) == 0.)
1656
goto triangular_error;
1659
Complex tmp = work[k] / data (cidx (kidx+1)-1);
1661
for (octave_idx_type i = cidx (kidx);
1662
i < cidx (kidx+1)-1; i++)
1664
octave_idx_type iidx = ridx (i);
1665
work[iidx] = work[iidx] - tmp * data (i);
1670
for (octave_idx_type i = 0; i < nc; i++)
1671
retval(perm[i], j) = work[i];
1676
// Calculation of 1-norm of inv(*this)
1677
for (octave_idx_type i = 0; i < nm; i++)
1680
for (octave_idx_type j = 0; j < nr; j++)
1684
for (octave_idx_type k = j; k >= 0; k--)
1686
octave_idx_type iidx = perm[k];
1690
Complex tmp = work[k] / data (cidx (iidx+1)-1);
1692
for (octave_idx_type i = cidx (iidx);
1693
i < cidx (iidx+1)-1; i++)
1695
octave_idx_type idx2 = ridx (i);
1696
work[idx2] = work[idx2] - tmp * data (i);
1701
for (octave_idx_type i = 0; i < j+1; i++)
1703
atmp += std::abs (work[i]);
1706
if (atmp > ainvnorm)
1709
rcond = 1. / ainvnorm / anorm;
1714
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1715
retval.resize (nc, b_nc);
1717
for (octave_idx_type j = 0; j < b_nc; j++)
1719
for (octave_idx_type i = 0; i < nr; i++)
1721
for (octave_idx_type i = nr; i < nc; i++)
1724
for (octave_idx_type k = nc-1; k >= 0; k--)
1728
if (ridx (cidx (k+1)-1) != k ||
1729
data (cidx (k+1)-1) == 0.)
1732
goto triangular_error;
1735
Complex tmp = work[k] / data (cidx (k+1)-1);
1737
for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1739
octave_idx_type iidx = ridx (i);
1740
work[iidx] = work[iidx] - tmp * data (i);
1745
for (octave_idx_type i = 0; i < nc; i++)
1746
retval.xelem (i, j) = work[i];
1751
// Calculation of 1-norm of inv(*this)
1752
for (octave_idx_type i = 0; i < nm; i++)
1755
for (octave_idx_type j = 0; j < nr; j++)
1759
for (octave_idx_type k = j; k >= 0; k--)
1763
Complex tmp = work[k] / data (cidx (k+1)-1);
1765
for (octave_idx_type i = cidx (k);
1766
i < cidx (k+1)-1; i++)
1768
octave_idx_type iidx = ridx (i);
1769
work[iidx] = work[iidx] - tmp * data (i);
1774
for (octave_idx_type i = 0; i < j+1; i++)
1776
atmp += std::abs (work[i]);
1779
if (atmp > ainvnorm)
1782
rcond = 1. / ainvnorm / anorm;
1791
sing_handler (rcond);
1792
mattype.mark_as_rectangular ();
1795
(*current_liboctave_error_handler)
1796
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
1800
volatile double rcond_plus_one = rcond + 1.0;
1802
if (rcond_plus_one == 1.0 || xisnan (rcond))
1808
sing_handler (rcond);
1809
mattype.mark_as_rectangular ();
1812
(*current_liboctave_error_handler)
1813
("matrix singular to machine precision, rcond = %g",
1818
(*current_liboctave_error_handler) ("incorrect matrix type");
1825
SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b,
1826
octave_idx_type& err, double& rcond,
1827
solve_singularity_handler sing_handler,
1828
bool calc_cond) const
1830
SparseComplexMatrix retval;
1832
octave_idx_type nr = rows ();
1833
octave_idx_type nc = cols ();
1834
octave_idx_type nm = (nc > nr ? nc : nr);
1837
if (nr != b.rows ())
1838
(*current_liboctave_error_handler)
1839
("matrix dimension mismatch solution of linear equations");
1840
else if (nr == 0 || nc == 0 || b.cols () == 0)
1841
retval = SparseComplexMatrix (nc, b.cols ());
1844
// Print spparms("spumoni") info if requested
1845
int typ = mattype.type ();
1848
if (typ == MatrixType::Permuted_Upper ||
1849
typ == MatrixType::Upper)
1852
double ainvnorm = 0.;
1857
// Calculate the 1-norm of matrix for rcond calculation
1858
for (octave_idx_type j = 0; j < nc; j++)
1861
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1862
atmp += std::abs (data (i));
1868
octave_idx_type b_nc = b.cols ();
1869
octave_idx_type b_nz = b.nnz ();
1870
retval = SparseComplexMatrix (nc, b_nc, b_nz);
1871
retval.xcidx (0) = 0;
1872
octave_idx_type ii = 0;
1873
octave_idx_type x_nz = b_nz;
1875
if (typ == MatrixType::Permuted_Upper)
1877
octave_idx_type *perm = mattype.triangular_perm ();
1878
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1880
OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1881
for (octave_idx_type i = 0; i < nc; i++)
1884
for (octave_idx_type j = 0; j < b_nc; j++)
1886
for (octave_idx_type i = 0; i < nm; i++)
1888
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1889
work[b.ridx (i)] = b.data (i);
1891
for (octave_idx_type k = nc-1; k >= 0; k--)
1893
octave_idx_type kidx = perm[k];
1897
if (ridx (cidx (kidx+1)-1) != k ||
1898
data (cidx (kidx+1)-1) == 0.)
1901
goto triangular_error;
1904
Complex tmp = work[k] / data (cidx (kidx+1)-1);
1906
for (octave_idx_type i = cidx (kidx);
1907
i < cidx (kidx+1)-1; i++)
1909
octave_idx_type iidx = ridx (i);
1910
work[iidx] = work[iidx] - tmp * data (i);
1915
// Count non-zeros in work vector and adjust space in
1917
octave_idx_type new_nnz = 0;
1918
for (octave_idx_type i = 0; i < nc; i++)
1922
if (ii + new_nnz > x_nz)
1924
// Resize the sparse matrix
1925
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1926
retval.change_capacity (sz);
1930
for (octave_idx_type i = 0; i < nc; i++)
1931
if (work[rperm[i]] != 0.)
1933
retval.xridx (ii) = i;
1934
retval.xdata (ii++) = work[rperm[i]];
1936
retval.xcidx (j+1) = ii;
1939
retval.maybe_compress ();
1943
// Calculation of 1-norm of inv(*this)
1944
for (octave_idx_type i = 0; i < nm; i++)
1947
for (octave_idx_type j = 0; j < nr; j++)
1951
for (octave_idx_type k = j; k >= 0; k--)
1953
octave_idx_type iidx = perm[k];
1957
Complex tmp = work[k] / data (cidx (iidx+1)-1);
1959
for (octave_idx_type i = cidx (iidx);
1960
i < cidx (iidx+1)-1; i++)
1962
octave_idx_type idx2 = ridx (i);
1963
work[idx2] = work[idx2] - tmp * data (i);
1968
for (octave_idx_type i = 0; i < j+1; i++)
1970
atmp += std::abs (work[i]);
1973
if (atmp > ainvnorm)
1976
rcond = 1. / ainvnorm / anorm;
1981
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
1983
for (octave_idx_type j = 0; j < b_nc; j++)
1985
for (octave_idx_type i = 0; i < nm; i++)
1987
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1988
work[b.ridx (i)] = b.data (i);
1990
for (octave_idx_type k = nc-1; k >= 0; k--)
1994
if (ridx (cidx (k+1)-1) != k ||
1995
data (cidx (k+1)-1) == 0.)
1998
goto triangular_error;
2001
Complex tmp = work[k] / data (cidx (k+1)-1);
2003
for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2005
octave_idx_type iidx = ridx (i);
2006
work[iidx] = work[iidx] - tmp * data (i);
2011
// Count non-zeros in work vector and adjust space in
2013
octave_idx_type new_nnz = 0;
2014
for (octave_idx_type i = 0; i < nc; i++)
2018
if (ii + new_nnz > x_nz)
2020
// Resize the sparse matrix
2021
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2022
retval.change_capacity (sz);
2026
for (octave_idx_type i = 0; i < nc; i++)
2029
retval.xridx (ii) = i;
2030
retval.xdata (ii++) = work[i];
2032
retval.xcidx (j+1) = ii;
2035
retval.maybe_compress ();
2039
// Calculation of 1-norm of inv(*this)
2040
for (octave_idx_type i = 0; i < nm; i++)
2043
for (octave_idx_type j = 0; j < nr; j++)
2047
for (octave_idx_type k = j; k >= 0; k--)
2051
Complex tmp = work[k] / data (cidx (k+1)-1);
2053
for (octave_idx_type i = cidx (k);
2054
i < cidx (k+1)-1; i++)
2056
octave_idx_type iidx = ridx (i);
2057
work[iidx] = work[iidx] - tmp * data (i);
2062
for (octave_idx_type i = 0; i < j+1; i++)
2064
atmp += std::abs (work[i]);
2067
if (atmp > ainvnorm)
2070
rcond = 1. / ainvnorm / anorm;
2079
sing_handler (rcond);
2080
mattype.mark_as_rectangular ();
2083
(*current_liboctave_error_handler)
2084
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
2088
volatile double rcond_plus_one = rcond + 1.0;
2090
if (rcond_plus_one == 1.0 || xisnan (rcond))
2096
sing_handler (rcond);
2097
mattype.mark_as_rectangular ();
2100
(*current_liboctave_error_handler)
2101
("matrix singular to machine precision, rcond = %g",
2106
(*current_liboctave_error_handler) ("incorrect matrix type");
2112
SparseComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
2113
octave_idx_type& err, double& rcond,
2114
solve_singularity_handler sing_handler,
2115
bool calc_cond) const
2117
ComplexMatrix retval;
2119
octave_idx_type nr = rows ();
2120
octave_idx_type nc = cols ();
2121
octave_idx_type nm = (nc > nr ? nc : nr);
2124
if (nr != b.rows ())
2125
(*current_liboctave_error_handler)
2126
("matrix dimension mismatch solution of linear equations");
2127
else if (nr == 0 || nc == 0 || b.cols () == 0)
2128
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2131
// Print spparms("spumoni") info if requested
2132
int typ = mattype.type ();
2135
if (typ == MatrixType::Permuted_Upper ||
2136
typ == MatrixType::Upper)
2139
double ainvnorm = 0.;
2140
octave_idx_type b_nc = b.cols ();
2145
// Calculate the 1-norm of matrix for rcond calculation
2146
for (octave_idx_type j = 0; j < nc; j++)
2149
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2150
atmp += std::abs (data (i));
2156
if (typ == MatrixType::Permuted_Upper)
2158
retval.resize (nc, b_nc);
2159
octave_idx_type *perm = mattype.triangular_perm ();
2160
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2162
for (octave_idx_type j = 0; j < b_nc; j++)
2164
for (octave_idx_type i = 0; i < nr; i++)
2166
for (octave_idx_type i = nr; i < nc; i++)
2169
for (octave_idx_type k = nc-1; k >= 0; k--)
2171
octave_idx_type kidx = perm[k];
2175
if (ridx (cidx (kidx+1)-1) != k ||
2176
data (cidx (kidx+1)-1) == 0.)
2179
goto triangular_error;
2182
Complex tmp = work[k] / data (cidx (kidx+1)-1);
2184
for (octave_idx_type i = cidx (kidx);
2185
i < cidx (kidx+1)-1; i++)
2187
octave_idx_type iidx = ridx (i);
2188
work[iidx] = work[iidx] - tmp * data (i);
2193
for (octave_idx_type i = 0; i < nc; i++)
2194
retval(perm[i], j) = work[i];
2199
// Calculation of 1-norm of inv(*this)
2200
for (octave_idx_type i = 0; i < nm; i++)
2203
for (octave_idx_type j = 0; j < nr; j++)
2207
for (octave_idx_type k = j; k >= 0; k--)
2209
octave_idx_type iidx = perm[k];
2213
Complex tmp = work[k] / data (cidx (iidx+1)-1);
2215
for (octave_idx_type i = cidx (iidx);
2216
i < cidx (iidx+1)-1; i++)
2218
octave_idx_type idx2 = ridx (i);
2219
work[idx2] = work[idx2] - tmp * data (i);
2224
for (octave_idx_type i = 0; i < j+1; i++)
2226
atmp += std::abs (work[i]);
2229
if (atmp > ainvnorm)
2232
rcond = 1. / ainvnorm / anorm;
2237
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2238
retval.resize (nc, b_nc);
2240
for (octave_idx_type j = 0; j < b_nc; j++)
2242
for (octave_idx_type i = 0; i < nr; i++)
2244
for (octave_idx_type i = nr; i < nc; i++)
2247
for (octave_idx_type k = nc-1; k >= 0; k--)
2251
if (ridx (cidx (k+1)-1) != k ||
2252
data (cidx (k+1)-1) == 0.)
2255
goto triangular_error;
2258
Complex tmp = work[k] / data (cidx (k+1)-1);
2260
for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2262
octave_idx_type iidx = ridx (i);
2263
work[iidx] = work[iidx] - tmp * data (i);
2268
for (octave_idx_type i = 0; i < nc; i++)
2269
retval.xelem (i, j) = work[i];
2274
// Calculation of 1-norm of inv(*this)
2275
for (octave_idx_type i = 0; i < nm; i++)
2278
for (octave_idx_type j = 0; j < nr; j++)
2282
for (octave_idx_type k = j; k >= 0; k--)
2286
Complex tmp = work[k] / data (cidx (k+1)-1);
2288
for (octave_idx_type i = cidx (k);
2289
i < cidx (k+1)-1; i++)
2291
octave_idx_type iidx = ridx (i);
2292
work[iidx] = work[iidx] - tmp * data (i);
2297
for (octave_idx_type i = 0; i < j+1; i++)
2299
atmp += std::abs (work[i]);
2302
if (atmp > ainvnorm)
2305
rcond = 1. / ainvnorm / anorm;
2314
sing_handler (rcond);
2315
mattype.mark_as_rectangular ();
2318
(*current_liboctave_error_handler)
2319
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
2323
volatile double rcond_plus_one = rcond + 1.0;
2325
if (rcond_plus_one == 1.0 || xisnan (rcond))
2331
sing_handler (rcond);
2332
mattype.mark_as_rectangular ();
2335
(*current_liboctave_error_handler)
2336
("matrix singular to machine precision, rcond = %g",
2341
(*current_liboctave_error_handler) ("incorrect matrix type");
2348
SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b,
2349
octave_idx_type& err, double& rcond,
2350
solve_singularity_handler sing_handler,
2351
bool calc_cond) const
2353
SparseComplexMatrix retval;
2355
octave_idx_type nr = rows ();
2356
octave_idx_type nc = cols ();
2357
octave_idx_type nm = (nc > nr ? nc : nr);
2360
if (nr != b.rows ())
2361
(*current_liboctave_error_handler)
2362
("matrix dimension mismatch solution of linear equations");
2363
else if (nr == 0 || nc == 0 || b.cols () == 0)
2364
retval = SparseComplexMatrix (nc, b.cols ());
2367
// Print spparms("spumoni") info if requested
2368
int typ = mattype.type ();
2371
if (typ == MatrixType::Permuted_Upper ||
2372
typ == MatrixType::Upper)
2375
double ainvnorm = 0.;
2380
// Calculate the 1-norm of matrix for rcond calculation
2381
for (octave_idx_type j = 0; j < nc; j++)
2384
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2385
atmp += std::abs (data (i));
2391
octave_idx_type b_nc = b.cols ();
2392
octave_idx_type b_nz = b.nnz ();
2393
retval = SparseComplexMatrix (nc, b_nc, b_nz);
2394
retval.xcidx (0) = 0;
2395
octave_idx_type ii = 0;
2396
octave_idx_type x_nz = b_nz;
2398
if (typ == MatrixType::Permuted_Upper)
2400
octave_idx_type *perm = mattype.triangular_perm ();
2401
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2403
OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
2404
for (octave_idx_type i = 0; i < nc; i++)
2407
for (octave_idx_type j = 0; j < b_nc; j++)
2409
for (octave_idx_type i = 0; i < nm; i++)
2411
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2412
work[b.ridx (i)] = b.data (i);
2414
for (octave_idx_type k = nc-1; k >= 0; k--)
2416
octave_idx_type kidx = perm[k];
2420
if (ridx (cidx (kidx+1)-1) != k ||
2421
data (cidx (kidx+1)-1) == 0.)
2424
goto triangular_error;
2427
Complex tmp = work[k] / data (cidx (kidx+1)-1);
2429
for (octave_idx_type i = cidx (kidx);
2430
i < cidx (kidx+1)-1; i++)
2432
octave_idx_type iidx = ridx (i);
2433
work[iidx] = work[iidx] - tmp * data (i);
2438
// Count non-zeros in work vector and adjust space in
2440
octave_idx_type new_nnz = 0;
2441
for (octave_idx_type i = 0; i < nc; i++)
2445
if (ii + new_nnz > x_nz)
2447
// Resize the sparse matrix
2448
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2449
retval.change_capacity (sz);
2453
for (octave_idx_type i = 0; i < nc; i++)
2454
if (work[rperm[i]] != 0.)
2456
retval.xridx (ii) = i;
2457
retval.xdata (ii++) = work[rperm[i]];
2459
retval.xcidx (j+1) = ii;
2462
retval.maybe_compress ();
2466
// Calculation of 1-norm of inv(*this)
2467
for (octave_idx_type i = 0; i < nm; i++)
2470
for (octave_idx_type j = 0; j < nr; j++)
2474
for (octave_idx_type k = j; k >= 0; k--)
2476
octave_idx_type iidx = perm[k];
2480
Complex tmp = work[k] / data (cidx (iidx+1)-1);
2482
for (octave_idx_type i = cidx (iidx);
2483
i < cidx (iidx+1)-1; i++)
2485
octave_idx_type idx2 = ridx (i);
2486
work[idx2] = work[idx2] - tmp * data (i);
2491
for (octave_idx_type i = 0; i < j+1; i++)
2493
atmp += std::abs (work[i]);
2496
if (atmp > ainvnorm)
2499
rcond = 1. / ainvnorm / anorm;
2504
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2506
for (octave_idx_type j = 0; j < b_nc; j++)
2508
for (octave_idx_type i = 0; i < nm; i++)
2510
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2511
work[b.ridx (i)] = b.data (i);
2513
for (octave_idx_type k = nr-1; k >= 0; k--)
2517
if (ridx (cidx (k+1)-1) != k ||
2518
data (cidx (k+1)-1) == 0.)
2521
goto triangular_error;
2524
Complex tmp = work[k] / data (cidx (k+1)-1);
2526
for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2528
octave_idx_type iidx = ridx (i);
2529
work[iidx] = work[iidx] - tmp * data (i);
2534
// Count non-zeros in work vector and adjust space in
2536
octave_idx_type new_nnz = 0;
2537
for (octave_idx_type i = 0; i < nc; i++)
2541
if (ii + new_nnz > x_nz)
2543
// Resize the sparse matrix
2544
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2545
retval.change_capacity (sz);
2549
for (octave_idx_type i = 0; i < nc; i++)
2552
retval.xridx (ii) = i;
2553
retval.xdata (ii++) = work[i];
2555
retval.xcidx (j+1) = ii;
2558
retval.maybe_compress ();
2562
// Calculation of 1-norm of inv(*this)
2563
for (octave_idx_type i = 0; i < nm; i++)
2566
for (octave_idx_type j = 0; j < nr; j++)
2570
for (octave_idx_type k = j; k >= 0; k--)
2574
Complex tmp = work[k] / data (cidx (k+1)-1);
2576
for (octave_idx_type i = cidx (k);
2577
i < cidx (k+1)-1; i++)
2579
octave_idx_type iidx = ridx (i);
2580
work[iidx] = work[iidx] - tmp * data (i);
2585
for (octave_idx_type i = 0; i < j+1; i++)
2587
atmp += std::abs (work[i]);
2590
if (atmp > ainvnorm)
2593
rcond = 1. / ainvnorm / anorm;
2602
sing_handler (rcond);
2603
mattype.mark_as_rectangular ();
2606
(*current_liboctave_error_handler)
2607
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
2611
volatile double rcond_plus_one = rcond + 1.0;
2613
if (rcond_plus_one == 1.0 || xisnan (rcond))
2619
sing_handler (rcond);
2620
mattype.mark_as_rectangular ();
2623
(*current_liboctave_error_handler)
2624
("matrix singular to machine precision, rcond = %g",
2629
(*current_liboctave_error_handler) ("incorrect matrix type");
2636
SparseComplexMatrix::ltsolve (MatrixType &mattype, const Matrix& b,
2637
octave_idx_type& err, double& rcond,
2638
solve_singularity_handler sing_handler,
2639
bool calc_cond) const
2641
ComplexMatrix retval;
2643
octave_idx_type nr = rows ();
2644
octave_idx_type nc = cols ();
2645
octave_idx_type nm = (nc > nr ? nc : nr);
2648
if (nr != b.rows ())
2649
(*current_liboctave_error_handler)
2650
("matrix dimension mismatch solution of linear equations");
2651
else if (nr == 0 || nc == 0 || b.cols () == 0)
2652
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2655
// Print spparms("spumoni") info if requested
2656
int typ = mattype.type ();
2659
if (typ == MatrixType::Permuted_Lower ||
2660
typ == MatrixType::Lower)
2663
double ainvnorm = 0.;
2664
octave_idx_type b_nc = b.cols ();
2669
// Calculate the 1-norm of matrix for rcond calculation
2670
for (octave_idx_type j = 0; j < nc; j++)
2673
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2674
atmp += std::abs (data (i));
2680
if (typ == MatrixType::Permuted_Lower)
2682
retval.resize (nc, b_nc);
2683
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2684
octave_idx_type *perm = mattype.triangular_perm ();
2686
for (octave_idx_type j = 0; j < b_nc; j++)
2688
for (octave_idx_type i = 0; i < nm; i++)
2690
for (octave_idx_type i = 0; i < nr; i++)
2691
work[perm[i]] = b(i,j);
2693
for (octave_idx_type k = 0; k < nc; k++)
2697
octave_idx_type minr = nr;
2698
octave_idx_type mini = 0;
2700
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2701
if (perm[ridx (i)] < minr)
2703
minr = perm[ridx (i)];
2707
if (minr != k || data (mini) == 0.)
2710
goto triangular_error;
2713
Complex tmp = work[k] / data (mini);
2715
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2720
octave_idx_type iidx = perm[ridx (i)];
2721
work[iidx] = work[iidx] - tmp * data (i);
2726
for (octave_idx_type i = 0; i < nc; i++)
2727
retval(i, j) = work[i];
2732
// Calculation of 1-norm of inv(*this)
2733
for (octave_idx_type i = 0; i < nm; i++)
2736
for (octave_idx_type j = 0; j < nr; j++)
2740
for (octave_idx_type k = 0; k < nc; k++)
2744
octave_idx_type minr = nr;
2745
octave_idx_type mini = 0;
2747
for (octave_idx_type i = cidx (k);
2748
i < cidx (k+1); i++)
2749
if (perm[ridx (i)] < minr)
2751
minr = perm[ridx (i)];
2755
Complex tmp = work[k] / data (mini);
2757
for (octave_idx_type i = cidx (k);
2758
i < cidx (k+1); i++)
2763
octave_idx_type iidx = perm[ridx (i)];
2764
work[iidx] = work[iidx] - tmp * data (i);
2770
for (octave_idx_type i = j; i < nc; i++)
2772
atmp += std::abs (work[i]);
2775
if (atmp > ainvnorm)
2778
rcond = 1. / ainvnorm / anorm;
2783
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2784
retval.resize (nc, b_nc, 0.);
2786
for (octave_idx_type j = 0; j < b_nc; j++)
2788
for (octave_idx_type i = 0; i < nr; i++)
2790
for (octave_idx_type i = nr; i < nc; i++)
2792
for (octave_idx_type k = 0; k < nc; k++)
2796
if (ridx (cidx (k)) != k ||
2797
data (cidx (k)) == 0.)
2800
goto triangular_error;
2803
Complex tmp = work[k] / data (cidx (k));
2805
for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
2807
octave_idx_type iidx = ridx (i);
2808
work[iidx] = work[iidx] - tmp * data (i);
2812
for (octave_idx_type i = 0; i < nc; i++)
2813
retval.xelem (i, j) = work[i];
2818
// Calculation of 1-norm of inv(*this)
2819
for (octave_idx_type i = 0; i < nm; i++)
2822
for (octave_idx_type j = 0; j < nr; j++)
2826
for (octave_idx_type k = j; k < nc; k++)
2831
Complex tmp = work[k] / data (cidx (k));
2833
for (octave_idx_type i = cidx (k)+1;
2834
i < cidx (k+1); i++)
2836
octave_idx_type iidx = ridx (i);
2837
work[iidx] = work[iidx] - tmp * data (i);
2842
for (octave_idx_type i = j; i < nc; i++)
2844
atmp += std::abs (work[i]);
2847
if (atmp > ainvnorm)
2850
rcond = 1. / ainvnorm / anorm;
2858
sing_handler (rcond);
2859
mattype.mark_as_rectangular ();
2862
(*current_liboctave_error_handler)
2863
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
2867
volatile double rcond_plus_one = rcond + 1.0;
2869
if (rcond_plus_one == 1.0 || xisnan (rcond))
2875
sing_handler (rcond);
2876
mattype.mark_as_rectangular ();
2879
(*current_liboctave_error_handler)
2880
("matrix singular to machine precision, rcond = %g",
2885
(*current_liboctave_error_handler) ("incorrect matrix type");
2892
SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b,
2893
octave_idx_type& err, double& rcond,
2894
solve_singularity_handler sing_handler,
2895
bool calc_cond) const
2897
SparseComplexMatrix retval;
2899
octave_idx_type nr = rows ();
2900
octave_idx_type nc = cols ();
2901
octave_idx_type nm = (nc > nr ? nc : nr);
2905
if (nr != b.rows ())
2906
(*current_liboctave_error_handler)
2907
("matrix dimension mismatch solution of linear equations");
2908
else if (nr == 0 || nc == 0 || b.cols () == 0)
2909
retval = SparseComplexMatrix (nc, b.cols ());
2912
// Print spparms("spumoni") info if requested
2913
int typ = mattype.type ();
2916
if (typ == MatrixType::Permuted_Lower ||
2917
typ == MatrixType::Lower)
2920
double ainvnorm = 0.;
2925
// Calculate the 1-norm of matrix for rcond calculation
2926
for (octave_idx_type j = 0; j < nc; j++)
2929
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2930
atmp += std::abs (data (i));
2936
octave_idx_type b_nc = b.cols ();
2937
octave_idx_type b_nz = b.nnz ();
2938
retval = SparseComplexMatrix (nc, b_nc, b_nz);
2939
retval.xcidx (0) = 0;
2940
octave_idx_type ii = 0;
2941
octave_idx_type x_nz = b_nz;
2943
if (typ == MatrixType::Permuted_Lower)
2945
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
2946
octave_idx_type *perm = mattype.triangular_perm ();
2948
for (octave_idx_type j = 0; j < b_nc; j++)
2950
for (octave_idx_type i = 0; i < nm; i++)
2952
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2953
work[perm[b.ridx (i)]] = b.data (i);
2955
for (octave_idx_type k = 0; k < nc; k++)
2959
octave_idx_type minr = nr;
2960
octave_idx_type mini = 0;
2962
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2963
if (perm[ridx (i)] < minr)
2965
minr = perm[ridx (i)];
2969
if (minr != k || data (mini) == 0.)
2972
goto triangular_error;
2975
Complex tmp = work[k] / data (mini);
2977
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2982
octave_idx_type iidx = perm[ridx (i)];
2983
work[iidx] = work[iidx] - tmp * data (i);
2988
// Count non-zeros in work vector and adjust space in
2990
octave_idx_type new_nnz = 0;
2991
for (octave_idx_type i = 0; i < nc; i++)
2995
if (ii + new_nnz > x_nz)
2997
// Resize the sparse matrix
2998
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2999
retval.change_capacity (sz);
3003
for (octave_idx_type i = 0; i < nc; i++)
3006
retval.xridx (ii) = i;
3007
retval.xdata (ii++) = work[i];
3009
retval.xcidx (j+1) = ii;
3012
retval.maybe_compress ();
3016
// Calculation of 1-norm of inv(*this)
3017
for (octave_idx_type i = 0; i < nm; i++)
3020
for (octave_idx_type j = 0; j < nr; j++)
3024
for (octave_idx_type k = 0; k < nc; k++)
3028
octave_idx_type minr = nr;
3029
octave_idx_type mini = 0;
3031
for (octave_idx_type i = cidx (k);
3032
i < cidx (k+1); i++)
3033
if (perm[ridx (i)] < minr)
3035
minr = perm[ridx (i)];
3039
Complex tmp = work[k] / data (mini);
3041
for (octave_idx_type i = cidx (k);
3042
i < cidx (k+1); i++)
3047
octave_idx_type iidx = perm[ridx (i)];
3048
work[iidx] = work[iidx] - tmp * data (i);
3054
for (octave_idx_type i = j; i < nc; i++)
3056
atmp += std::abs (work[i]);
3059
if (atmp > ainvnorm)
3062
rcond = 1. / ainvnorm / anorm;
3067
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3069
for (octave_idx_type j = 0; j < b_nc; j++)
3071
for (octave_idx_type i = 0; i < nm; i++)
3073
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3074
work[b.ridx (i)] = b.data (i);
3076
for (octave_idx_type k = 0; k < nc; k++)
3080
if (ridx (cidx (k)) != k ||
3081
data (cidx (k)) == 0.)
3084
goto triangular_error;
3087
Complex tmp = work[k] / data (cidx (k));
3089
for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3091
octave_idx_type iidx = ridx (i);
3092
work[iidx] = work[iidx] - tmp * data (i);
3097
// Count non-zeros in work vector and adjust space in
3099
octave_idx_type new_nnz = 0;
3100
for (octave_idx_type i = 0; i < nc; i++)
3104
if (ii + new_nnz > x_nz)
3106
// Resize the sparse matrix
3107
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3108
retval.change_capacity (sz);
3112
for (octave_idx_type i = 0; i < nc; i++)
3115
retval.xridx (ii) = i;
3116
retval.xdata (ii++) = work[i];
3118
retval.xcidx (j+1) = ii;
3121
retval.maybe_compress ();
3125
// Calculation of 1-norm of inv(*this)
3126
for (octave_idx_type i = 0; i < nm; i++)
3129
for (octave_idx_type j = 0; j < nr; j++)
3133
for (octave_idx_type k = j; k < nc; k++)
3138
Complex tmp = work[k] / data (cidx (k));
3140
for (octave_idx_type i = cidx (k)+1;
3141
i < cidx (k+1); i++)
3143
octave_idx_type iidx = ridx (i);
3144
work[iidx] = work[iidx] - tmp * data (i);
3149
for (octave_idx_type i = j; i < nc; i++)
3151
atmp += std::abs (work[i]);
3154
if (atmp > ainvnorm)
3157
rcond = 1. / ainvnorm / anorm;
3166
sing_handler (rcond);
3167
mattype.mark_as_rectangular ();
3170
(*current_liboctave_error_handler)
3171
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
3175
volatile double rcond_plus_one = rcond + 1.0;
3177
if (rcond_plus_one == 1.0 || xisnan (rcond))
3183
sing_handler (rcond);
3184
mattype.mark_as_rectangular ();
3187
(*current_liboctave_error_handler)
3188
("matrix singular to machine precision, rcond = %g",
3193
(*current_liboctave_error_handler) ("incorrect matrix type");
3200
SparseComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
3201
octave_idx_type& err, double& rcond,
3202
solve_singularity_handler sing_handler,
3203
bool calc_cond) const
3205
ComplexMatrix retval;
3207
octave_idx_type nr = rows ();
3208
octave_idx_type nc = cols ();
3209
octave_idx_type nm = (nc > nr ? nc : nr);
3212
if (nr != b.rows ())
3213
(*current_liboctave_error_handler)
3214
("matrix dimension mismatch solution of linear equations");
3215
else if (nr == 0 || nc == 0 || b.cols () == 0)
3216
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3219
// Print spparms("spumoni") info if requested
3220
int typ = mattype.type ();
3223
if (typ == MatrixType::Permuted_Lower ||
3224
typ == MatrixType::Lower)
3227
double ainvnorm = 0.;
3228
octave_idx_type b_nc = b.cols ();
3233
// Calculate the 1-norm of matrix for rcond calculation
3234
for (octave_idx_type j = 0; j < nc; j++)
3237
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3238
atmp += std::abs (data (i));
3244
if (typ == MatrixType::Permuted_Lower)
3246
retval.resize (nc, b_nc);
3247
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3248
octave_idx_type *perm = mattype.triangular_perm ();
3250
for (octave_idx_type j = 0; j < b_nc; j++)
3252
for (octave_idx_type i = 0; i < nm; i++)
3254
for (octave_idx_type i = 0; i < nr; i++)
3255
work[perm[i]] = b(i,j);
3257
for (octave_idx_type k = 0; k < nc; k++)
3261
octave_idx_type minr = nr;
3262
octave_idx_type mini = 0;
3264
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3265
if (perm[ridx (i)] < minr)
3267
minr = perm[ridx (i)];
3271
if (minr != k || data (mini) == 0.)
3274
goto triangular_error;
3277
Complex tmp = work[k] / data (mini);
3279
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3284
octave_idx_type iidx = perm[ridx (i)];
3285
work[iidx] = work[iidx] - tmp * data (i);
3290
for (octave_idx_type i = 0; i < nc; i++)
3291
retval(i, j) = work[i];
3296
// Calculation of 1-norm of inv(*this)
3297
for (octave_idx_type i = 0; i < nm; i++)
3300
for (octave_idx_type j = 0; j < nr; j++)
3304
for (octave_idx_type k = 0; k < nc; k++)
3308
octave_idx_type minr = nr;
3309
octave_idx_type mini = 0;
3311
for (octave_idx_type i = cidx (k);
3312
i < cidx (k+1); i++)
3313
if (perm[ridx (i)] < minr)
3315
minr = perm[ridx (i)];
3319
Complex tmp = work[k] / data (mini);
3321
for (octave_idx_type i = cidx (k);
3322
i < cidx (k+1); i++)
3327
octave_idx_type iidx = perm[ridx (i)];
3328
work[iidx] = work[iidx] - tmp * data (i);
3334
for (octave_idx_type i = j; i < nc; i++)
3336
atmp += std::abs (work[i]);
3339
if (atmp > ainvnorm)
3342
rcond = 1. / ainvnorm / anorm;
3347
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3348
retval.resize (nc, b_nc, 0.);
3351
for (octave_idx_type j = 0; j < b_nc; j++)
3353
for (octave_idx_type i = 0; i < nr; i++)
3355
for (octave_idx_type i = nr; i < nc; i++)
3358
for (octave_idx_type k = 0; k < nc; k++)
3362
if (ridx (cidx (k)) != k ||
3363
data (cidx (k)) == 0.)
3366
goto triangular_error;
3369
Complex tmp = work[k] / data (cidx (k));
3371
for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3373
octave_idx_type iidx = ridx (i);
3374
work[iidx] = work[iidx] - tmp * data (i);
3379
for (octave_idx_type i = 0; i < nc; i++)
3380
retval.xelem (i, j) = work[i];
3385
// Calculation of 1-norm of inv(*this)
3386
for (octave_idx_type i = 0; i < nm; i++)
3389
for (octave_idx_type j = 0; j < nr; j++)
3393
for (octave_idx_type k = j; k < nc; k++)
3398
Complex tmp = work[k] / data (cidx (k));
3400
for (octave_idx_type i = cidx (k)+1;
3401
i < cidx (k+1); i++)
3403
octave_idx_type iidx = ridx (i);
3404
work[iidx] = work[iidx] - tmp * data (i);
3409
for (octave_idx_type i = j; i < nc; i++)
3411
atmp += std::abs (work[i]);
3414
if (atmp > ainvnorm)
3417
rcond = 1. / ainvnorm / anorm;
3426
sing_handler (rcond);
3427
mattype.mark_as_rectangular ();
3430
(*current_liboctave_error_handler)
3431
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
3435
volatile double rcond_plus_one = rcond + 1.0;
3437
if (rcond_plus_one == 1.0 || xisnan (rcond))
3443
sing_handler (rcond);
3444
mattype.mark_as_rectangular ();
3447
(*current_liboctave_error_handler)
3448
("matrix singular to machine precision, rcond = %g",
3453
(*current_liboctave_error_handler) ("incorrect matrix type");
3460
SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseComplexMatrix& b,
3461
octave_idx_type& err, double& rcond,
3462
solve_singularity_handler sing_handler,
3463
bool calc_cond) const
3465
SparseComplexMatrix retval;
3467
octave_idx_type nr = rows ();
3468
octave_idx_type nc = cols ();
3469
octave_idx_type nm = (nc > nr ? nc : nr);
3472
if (nr != b.rows ())
3473
(*current_liboctave_error_handler)
3474
("matrix dimension mismatch solution of linear equations");
3475
else if (nr == 0 || nc == 0 || b.cols () == 0)
3476
retval = SparseComplexMatrix (nc, b.cols ());
3479
// Print spparms("spumoni") info if requested
3480
int typ = mattype.type ();
3483
if (typ == MatrixType::Permuted_Lower ||
3484
typ == MatrixType::Lower)
3487
double ainvnorm = 0.;
3492
// Calculate the 1-norm of matrix for rcond calculation
3493
for (octave_idx_type j = 0; j < nc; j++)
3496
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3497
atmp += std::abs (data (i));
3503
octave_idx_type b_nc = b.cols ();
3504
octave_idx_type b_nz = b.nnz ();
3505
retval = SparseComplexMatrix (nc, b_nc, b_nz);
3506
retval.xcidx (0) = 0;
3507
octave_idx_type ii = 0;
3508
octave_idx_type x_nz = b_nz;
3510
if (typ == MatrixType::Permuted_Lower)
3512
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3513
octave_idx_type *perm = mattype.triangular_perm ();
3515
for (octave_idx_type j = 0; j < b_nc; j++)
3517
for (octave_idx_type i = 0; i < nm; i++)
3519
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3520
work[perm[b.ridx (i)]] = b.data (i);
3522
for (octave_idx_type k = 0; k < nc; k++)
3526
octave_idx_type minr = nr;
3527
octave_idx_type mini = 0;
3529
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3530
if (perm[ridx (i)] < minr)
3532
minr = perm[ridx (i)];
3536
if (minr != k || data (mini) == 0.)
3539
goto triangular_error;
3542
Complex tmp = work[k] / data (mini);
3544
for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3549
octave_idx_type iidx = perm[ridx (i)];
3550
work[iidx] = work[iidx] - tmp * data (i);
3555
// Count non-zeros in work vector and adjust space in
3557
octave_idx_type new_nnz = 0;
3558
for (octave_idx_type i = 0; i < nc; i++)
3562
if (ii + new_nnz > x_nz)
3564
// Resize the sparse matrix
3565
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3566
retval.change_capacity (sz);
3570
for (octave_idx_type i = 0; i < nc; i++)
3573
retval.xridx (ii) = i;
3574
retval.xdata (ii++) = work[i];
3576
retval.xcidx (j+1) = ii;
3579
retval.maybe_compress ();
3583
// Calculation of 1-norm of inv(*this)
3584
for (octave_idx_type i = 0; i < nm; i++)
3587
for (octave_idx_type j = 0; j < nr; j++)
3591
for (octave_idx_type k = 0; k < nc; k++)
3595
octave_idx_type minr = nr;
3596
octave_idx_type mini = 0;
3598
for (octave_idx_type i = cidx (k);
3599
i < cidx (k+1); i++)
3600
if (perm[ridx (i)] < minr)
3602
minr = perm[ridx (i)];
3606
Complex tmp = work[k] / data (mini);
3608
for (octave_idx_type i = cidx (k);
3609
i < cidx (k+1); i++)
3614
octave_idx_type iidx = perm[ridx (i)];
3615
work[iidx] = work[iidx] - tmp * data (i);
3621
for (octave_idx_type i = j; i < nc; i++)
3623
atmp += std::abs (work[i]);
3626
if (atmp > ainvnorm)
3629
rcond = 1. / ainvnorm / anorm;
3634
OCTAVE_LOCAL_BUFFER (Complex, work, nm);
3636
for (octave_idx_type j = 0; j < b_nc; j++)
3638
for (octave_idx_type i = 0; i < nm; i++)
3640
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3641
work[b.ridx (i)] = b.data (i);
3643
for (octave_idx_type k = 0; k < nc; k++)
3647
if (ridx (cidx (k)) != k ||
3648
data (cidx (k)) == 0.)
3651
goto triangular_error;
3654
Complex tmp = work[k] / data (cidx (k));
3656
for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3658
octave_idx_type iidx = ridx (i);
3659
work[iidx] = work[iidx] - tmp * data (i);
3664
// Count non-zeros in work vector and adjust space in
3666
octave_idx_type new_nnz = 0;
3667
for (octave_idx_type i = 0; i < nc; i++)
3671
if (ii + new_nnz > x_nz)
3673
// Resize the sparse matrix
3674
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3675
retval.change_capacity (sz);
3679
for (octave_idx_type i = 0; i < nc; i++)
3682
retval.xridx (ii) = i;
3683
retval.xdata (ii++) = work[i];
3685
retval.xcidx (j+1) = ii;
3688
retval.maybe_compress ();
3692
// Calculation of 1-norm of inv(*this)
3693
for (octave_idx_type i = 0; i < nm; i++)
3696
for (octave_idx_type j = 0; j < nr; j++)
3700
for (octave_idx_type k = j; k < nc; k++)
3705
Complex tmp = work[k] / data (cidx (k));
3707
for (octave_idx_type i = cidx (k)+1;
3708
i < cidx (k+1); i++)
3710
octave_idx_type iidx = ridx (i);
3711
work[iidx] = work[iidx] - tmp * data (i);
3716
for (octave_idx_type i = j; i < nc; i++)
3718
atmp += std::abs (work[i]);
3721
if (atmp > ainvnorm)
3724
rcond = 1. / ainvnorm / anorm;
3733
sing_handler (rcond);
3734
mattype.mark_as_rectangular ();
3737
(*current_liboctave_error_handler)
3738
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
3742
volatile double rcond_plus_one = rcond + 1.0;
3744
if (rcond_plus_one == 1.0 || xisnan (rcond))
3750
sing_handler (rcond);
3751
mattype.mark_as_rectangular ();
3754
(*current_liboctave_error_handler)
3755
("matrix singular to machine precision, rcond = %g",
3760
(*current_liboctave_error_handler) ("incorrect matrix type");
3767
SparseComplexMatrix::trisolve (MatrixType &mattype, const Matrix& b,
3768
octave_idx_type& err, double& rcond,
3769
solve_singularity_handler sing_handler,
3770
bool calc_cond) const
3772
ComplexMatrix retval;
3774
octave_idx_type nr = rows ();
3775
octave_idx_type nc = cols ();
3778
if (nr != nc || nr != b.rows ())
3779
(*current_liboctave_error_handler)
3780
("matrix dimension mismatch solution of linear equations");
3781
else if (nr == 0 || b.cols () == 0)
3782
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3784
(*current_liboctave_error_handler)
3785
("calculation of condition number not implemented");
3788
// Print spparms("spumoni") info if requested
3789
volatile int typ = mattype.type ();
3792
if (typ == MatrixType::Tridiagonal_Hermitian)
3794
OCTAVE_LOCAL_BUFFER (double, D, nr);
3795
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3797
if (mattype.is_dense ())
3799
octave_idx_type ii = 0;
3801
for (octave_idx_type j = 0; j < nc-1; j++)
3803
D[j] = std::real (data (ii++));
3807
D[nc-1] = std::real (data (ii));
3812
for (octave_idx_type i = 0; i < nr - 1; i++)
3818
for (octave_idx_type j = 0; j < nc; j++)
3819
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3822
D[j] = std::real (data (i));
3823
else if (ridx (i) == j + 1)
3828
octave_idx_type b_nc = b.cols ();
3829
retval = ComplexMatrix (b);
3830
Complex *result = retval.fortran_vec ();
3832
F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
3838
mattype.mark_as_unsymmetric ();
3839
typ = MatrixType::Tridiagonal;
3845
if (typ == MatrixType::Tridiagonal)
3847
OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3848
OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3849
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3851
if (mattype.is_dense ())
3853
octave_idx_type ii = 0;
3855
for (octave_idx_type j = 0; j < nc-1; j++)
3858
DL[j] = data (ii++);
3859
DU[j] = data (ii++);
3861
D[nc-1] = data (ii);
3866
for (octave_idx_type i = 0; i < nr - 1; i++)
3873
for (octave_idx_type j = 0; j < nc; j++)
3874
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3878
else if (ridx (i) == j + 1)
3880
else if (ridx (i) == j - 1)
3885
octave_idx_type b_nc = b.cols ();
3886
retval = ComplexMatrix (b);
3887
Complex *result = retval.fortran_vec ();
3889
F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
3899
sing_handler (rcond);
3900
mattype.mark_as_rectangular ();
3903
(*current_liboctave_error_handler)
3904
("matrix singular to machine precision");
3910
else if (typ != MatrixType::Tridiagonal_Hermitian)
3911
(*current_liboctave_error_handler) ("incorrect matrix type");
3918
SparseComplexMatrix::trisolve (MatrixType &mattype, const SparseMatrix& b,
3919
octave_idx_type& err, double& rcond,
3920
solve_singularity_handler sing_handler,
3921
bool calc_cond) const
3923
SparseComplexMatrix retval;
3925
octave_idx_type nr = rows ();
3926
octave_idx_type nc = cols ();
3929
if (nr != nc || nr != b.rows ())
3930
(*current_liboctave_error_handler)
3931
("matrix dimension mismatch solution of linear equations");
3932
else if (nr == 0 || b.cols () == 0)
3933
retval = SparseComplexMatrix (nc, b.cols ());
3935
(*current_liboctave_error_handler)
3936
("calculation of condition number not implemented");
3939
// Print spparms("spumoni") info if requested
3940
int typ = mattype.type ();
3943
// Note can't treat symmetric case as there is no dpttrf function
3944
if (typ == MatrixType::Tridiagonal ||
3945
typ == MatrixType::Tridiagonal_Hermitian)
3947
OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
3948
OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
3949
OCTAVE_LOCAL_BUFFER (Complex, D, nr);
3950
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
3951
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
3952
octave_idx_type *pipvt = ipvt.fortran_vec ();
3954
if (mattype.is_dense ())
3956
octave_idx_type ii = 0;
3958
for (octave_idx_type j = 0; j < nc-1; j++)
3961
DL[j] = data (ii++);
3962
DU[j] = data (ii++);
3964
D[nc-1] = data (ii);
3969
for (octave_idx_type i = 0; i < nr - 1; i++)
3976
for (octave_idx_type j = 0; j < nc; j++)
3977
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3981
else if (ridx (i) == j + 1)
3983
else if (ridx (i) == j - 1)
3988
F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
3997
sing_handler (rcond);
3998
mattype.mark_as_rectangular ();
4001
(*current_liboctave_error_handler)
4002
("matrix singular to machine precision");
4008
volatile octave_idx_type x_nz = b.nnz ();
4009
octave_idx_type b_nc = b.cols ();
4010
retval = SparseComplexMatrix (nr, b_nc, x_nz);
4011
retval.xcidx (0) = 0;
4012
volatile octave_idx_type ii = 0;
4015
OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4017
for (volatile octave_idx_type j = 0; j < b_nc; j++)
4019
for (octave_idx_type i = 0; i < nr; i++)
4021
for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
4022
work[b.ridx (i)] = b.data (i);
4024
F77_XFCN (zgttrs, ZGTTRS,
4025
(F77_CONST_CHAR_ARG2 (&job, 1),
4026
nr, 1, DL, D, DU, DU2, pipvt,
4027
work, b.rows (), err
4028
F77_CHAR_ARG_LEN (1)));
4030
// Count non-zeros in work vector and adjust
4031
// space in retval if needed
4032
octave_idx_type new_nnz = 0;
4033
for (octave_idx_type i = 0; i < nr; i++)
4037
if (ii + new_nnz > x_nz)
4039
// Resize the sparse matrix
4040
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4041
retval.change_capacity (sz);
4045
for (octave_idx_type i = 0; i < nr; i++)
4048
retval.xridx (ii) = i;
4049
retval.xdata (ii++) = work[i];
4051
retval.xcidx (j+1) = ii;
4054
retval.maybe_compress ();
4057
else if (typ != MatrixType::Tridiagonal_Hermitian)
4058
(*current_liboctave_error_handler) ("incorrect matrix type");
4065
SparseComplexMatrix::trisolve (MatrixType &mattype, const ComplexMatrix& b,
4066
octave_idx_type& err, double& rcond,
4067
solve_singularity_handler sing_handler,
4068
bool calc_cond) const
4070
ComplexMatrix retval;
4072
octave_idx_type nr = rows ();
4073
octave_idx_type nc = cols ();
4076
if (nr != nc || nr != b.rows ())
4077
(*current_liboctave_error_handler)
4078
("matrix dimension mismatch solution of linear equations");
4079
else if (nr == 0 || b.cols () == 0)
4080
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4082
(*current_liboctave_error_handler)
4083
("calculation of condition number not implemented");
4086
// Print spparms("spumoni") info if requested
4087
volatile int typ = mattype.type ();
4090
if (typ == MatrixType::Tridiagonal_Hermitian)
4092
OCTAVE_LOCAL_BUFFER (double, D, nr);
4093
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4095
if (mattype.is_dense ())
4097
octave_idx_type ii = 0;
4099
for (octave_idx_type j = 0; j < nc-1; j++)
4101
D[j] = std::real (data (ii++));
4105
D[nc-1] = std::real (data (ii));
4110
for (octave_idx_type i = 0; i < nr - 1; i++)
4116
for (octave_idx_type j = 0; j < nc; j++)
4117
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4120
D[j] = std::real (data (i));
4121
else if (ridx (i) == j + 1)
4126
octave_idx_type b_nr = b.rows ();
4127
octave_idx_type b_nc = b.cols ();
4130
retval = ComplexMatrix (b);
4131
Complex *result = retval.fortran_vec ();
4133
F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
4139
mattype.mark_as_unsymmetric ();
4140
typ = MatrixType::Tridiagonal;
4144
if (typ == MatrixType::Tridiagonal)
4146
OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4147
OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4148
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4150
if (mattype.is_dense ())
4152
octave_idx_type ii = 0;
4154
for (octave_idx_type j = 0; j < nc-1; j++)
4157
DL[j] = data (ii++);
4158
DU[j] = data (ii++);
4160
D[nc-1] = data (ii);
4165
for (octave_idx_type i = 0; i < nr - 1; i++)
4172
for (octave_idx_type j = 0; j < nc; j++)
4173
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4177
else if (ridx (i) == j + 1)
4179
else if (ridx (i) == j - 1)
4184
octave_idx_type b_nr = b.rows ();
4185
octave_idx_type b_nc = b.cols ();
4188
retval = ComplexMatrix (b);
4189
Complex *result = retval.fortran_vec ();
4191
F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
4201
sing_handler (rcond);
4202
mattype.mark_as_rectangular ();
4205
(*current_liboctave_error_handler)
4206
("matrix singular to machine precision");
4209
else if (typ != MatrixType::Tridiagonal_Hermitian)
4210
(*current_liboctave_error_handler) ("incorrect matrix type");
4217
SparseComplexMatrix::trisolve (MatrixType &mattype,
4218
const SparseComplexMatrix& b,
4219
octave_idx_type& err, double& rcond,
4220
solve_singularity_handler sing_handler,
4221
bool calc_cond) const
4223
SparseComplexMatrix retval;
4225
octave_idx_type nr = rows ();
4226
octave_idx_type nc = cols ();
4229
if (nr != nc || nr != b.rows ())
4230
(*current_liboctave_error_handler)
4231
("matrix dimension mismatch solution of linear equations");
4232
else if (nr == 0 || b.cols () == 0)
4233
retval = SparseComplexMatrix (nc, b.cols ());
4235
(*current_liboctave_error_handler)
4236
("calculation of condition number not implemented");
4239
// Print spparms("spumoni") info if requested
4240
int typ = mattype.type ();
4243
// Note can't treat symmetric case as there is no dpttrf function
4244
if (typ == MatrixType::Tridiagonal ||
4245
typ == MatrixType::Tridiagonal_Hermitian)
4247
OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2);
4248
OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1);
4249
OCTAVE_LOCAL_BUFFER (Complex, D, nr);
4250
OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
4251
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4252
octave_idx_type *pipvt = ipvt.fortran_vec ();
4254
if (mattype.is_dense ())
4256
octave_idx_type ii = 0;
4258
for (octave_idx_type j = 0; j < nc-1; j++)
4261
DL[j] = data (ii++);
4262
DU[j] = data (ii++);
4264
D[nc-1] = data (ii);
4269
for (octave_idx_type i = 0; i < nr - 1; i++)
4276
for (octave_idx_type j = 0; j < nc; j++)
4277
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4281
else if (ridx (i) == j + 1)
4283
else if (ridx (i) == j - 1)
4288
F77_XFCN (zgttrf, ZGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
4297
sing_handler (rcond);
4298
mattype.mark_as_rectangular ();
4301
(*current_liboctave_error_handler)
4302
("matrix singular to machine precision");
4308
octave_idx_type b_nr = b.rows ();
4309
octave_idx_type b_nc = b.cols ();
4310
OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4312
// Take a first guess that the number of non-zero terms
4313
// will be as many as in b
4314
volatile octave_idx_type x_nz = b.nnz ();
4315
volatile octave_idx_type ii = 0;
4316
retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4318
retval.xcidx (0) = 0;
4319
for (volatile octave_idx_type j = 0; j < b_nc; j++)
4322
for (octave_idx_type i = 0; i < b_nr; i++)
4325
F77_XFCN (zgttrs, ZGTTRS,
4326
(F77_CONST_CHAR_ARG2 (&job, 1),
4327
nr, 1, DL, D, DU, DU2, pipvt,
4329
F77_CHAR_ARG_LEN (1)));
4333
(*current_liboctave_error_handler)
4334
("SparseComplexMatrix::solve solve failed");
4340
// Count non-zeros in work vector and adjust
4341
// space in retval if needed
4342
octave_idx_type new_nnz = 0;
4343
for (octave_idx_type i = 0; i < nr; i++)
4347
if (ii + new_nnz > x_nz)
4349
// Resize the sparse matrix
4350
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4351
retval.change_capacity (sz);
4355
for (octave_idx_type i = 0; i < nr; i++)
4358
retval.xridx (ii) = i;
4359
retval.xdata (ii++) = Bx[i];
4362
retval.xcidx (j+1) = ii;
4365
retval.maybe_compress ();
4368
else if (typ != MatrixType::Tridiagonal_Hermitian)
4369
(*current_liboctave_error_handler) ("incorrect matrix type");
4376
SparseComplexMatrix::bsolve (MatrixType &mattype, const Matrix& b,
4377
octave_idx_type& err, double& rcond,
4378
solve_singularity_handler sing_handler,
4379
bool calc_cond) const
4381
ComplexMatrix retval;
4383
octave_idx_type nr = rows ();
4384
octave_idx_type nc = cols ();
4387
if (nr != nc || nr != b.rows ())
4388
(*current_liboctave_error_handler)
4389
("matrix dimension mismatch solution of linear equations");
4390
else if (nr == 0 || b.cols () == 0)
4391
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4394
// Print spparms("spumoni") info if requested
4395
volatile int typ = mattype.type ();
4398
if (typ == MatrixType::Banded_Hermitian)
4400
octave_idx_type n_lower = mattype.nlower ();
4401
octave_idx_type ldm = n_lower + 1;
4402
ComplexMatrix m_band (ldm, nc);
4403
Complex *tmp_data = m_band.fortran_vec ();
4405
if (! mattype.is_dense ())
4407
octave_idx_type ii = 0;
4409
for (octave_idx_type j = 0; j < ldm; j++)
4410
for (octave_idx_type i = 0; i < nc; i++)
4411
tmp_data[ii++] = 0.;
4414
for (octave_idx_type j = 0; j < nc; j++)
4415
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4417
octave_idx_type ri = ridx (i);
4419
m_band(ri - j, j) = data (i);
4422
// Calculate the norm of the matrix, for later use.
4425
anorm = m_band.abs ().sum ().row (0).max ();
4428
F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4429
nr, n_lower, tmp_data, ldm, err
4430
F77_CHAR_ARG_LEN (1)));
4435
// Matrix is not positive definite!! Fall through to
4436
// unsymmetric banded solver.
4437
mattype.mark_as_unsymmetric ();
4438
typ = MatrixType::Banded;
4445
Array<Complex> z (dim_vector (2 * nr, 1));
4446
Complex *pz = z.fortran_vec ();
4447
Array<double> iz (dim_vector (nr, 1));
4448
double *piz = iz.fortran_vec ();
4450
F77_XFCN (zpbcon, ZPBCON,
4451
(F77_CONST_CHAR_ARG2 (&job, 1),
4452
nr, n_lower, tmp_data, ldm,
4453
anorm, rcond, pz, piz, err
4454
F77_CHAR_ARG_LEN (1)));
4459
volatile double rcond_plus_one = rcond + 1.0;
4461
if (rcond_plus_one == 1.0 || xisnan (rcond))
4467
sing_handler (rcond);
4468
mattype.mark_as_rectangular ();
4471
(*current_liboctave_error_handler)
4472
("matrix singular to machine precision, rcond = %g",
4481
retval = ComplexMatrix (b);
4482
Complex *result = retval.fortran_vec ();
4484
octave_idx_type b_nc = b.cols ();
4486
F77_XFCN (zpbtrs, ZPBTRS,
4487
(F77_CONST_CHAR_ARG2 (&job, 1),
4488
nr, n_lower, b_nc, tmp_data,
4489
ldm, result, b.rows (), err
4490
F77_CHAR_ARG_LEN (1)));
4494
(*current_liboctave_error_handler)
4495
("SparseMatrix::solve solve failed");
4502
if (typ == MatrixType::Banded)
4504
// Create the storage for the banded form of the sparse matrix
4505
octave_idx_type n_upper = mattype.nupper ();
4506
octave_idx_type n_lower = mattype.nlower ();
4507
octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4509
ComplexMatrix m_band (ldm, nc);
4510
Complex *tmp_data = m_band.fortran_vec ();
4512
if (! mattype.is_dense ())
4514
octave_idx_type ii = 0;
4516
for (octave_idx_type j = 0; j < ldm; j++)
4517
for (octave_idx_type i = 0; i < nc; i++)
4518
tmp_data[ii++] = 0.;
4521
for (octave_idx_type j = 0; j < nc; j++)
4522
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4523
m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4525
// Calculate the norm of the matrix, for later use.
4529
for (octave_idx_type j = 0; j < nr; j++)
4532
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4533
atmp += std::abs (data (i));
4539
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4540
octave_idx_type *pipvt = ipvt.fortran_vec ();
4542
F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data,
4545
// Throw-away extra info LAPACK gives so as to not
4554
sing_handler (rcond);
4555
mattype.mark_as_rectangular ();
4558
(*current_liboctave_error_handler)
4559
("matrix singular to machine precision");
4566
Array<Complex> z (dim_vector (2 * nr, 1));
4567
Complex *pz = z.fortran_vec ();
4568
Array<double> iz (dim_vector (nr, 1));
4569
double *piz = iz.fortran_vec ();
4571
F77_XFCN (zgbcon, ZGBCON,
4572
(F77_CONST_CHAR_ARG2 (&job, 1),
4573
nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4574
anorm, rcond, pz, piz, err
4575
F77_CHAR_ARG_LEN (1)));
4580
volatile double rcond_plus_one = rcond + 1.0;
4582
if (rcond_plus_one == 1.0 || xisnan (rcond))
4588
sing_handler (rcond);
4589
mattype.mark_as_rectangular ();
4592
(*current_liboctave_error_handler)
4593
("matrix singular to machine precision, rcond = %g",
4602
retval = ComplexMatrix (b);
4603
Complex *result = retval.fortran_vec ();
4605
octave_idx_type b_nc = b.cols ();
4608
F77_XFCN (zgbtrs, ZGBTRS,
4609
(F77_CONST_CHAR_ARG2 (&job, 1),
4610
nr, n_lower, n_upper, b_nc, tmp_data,
4611
ldm, pipvt, result, b.rows (), err
4612
F77_CHAR_ARG_LEN (1)));
4616
else if (typ != MatrixType::Banded_Hermitian)
4617
(*current_liboctave_error_handler) ("incorrect matrix type");
4624
SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseMatrix& b,
4625
octave_idx_type& err, double& rcond,
4626
solve_singularity_handler sing_handler,
4627
bool calc_cond) const
4629
SparseComplexMatrix retval;
4631
octave_idx_type nr = rows ();
4632
octave_idx_type nc = cols ();
4635
if (nr != nc || nr != b.rows ())
4636
(*current_liboctave_error_handler)
4637
("matrix dimension mismatch solution of linear equations");
4638
else if (nr == 0 || b.cols () == 0)
4639
retval = SparseComplexMatrix (nc, b.cols ());
4642
// Print spparms("spumoni") info if requested
4643
volatile int typ = mattype.type ();
4646
if (typ == MatrixType::Banded_Hermitian)
4648
octave_idx_type n_lower = mattype.nlower ();
4649
octave_idx_type ldm = n_lower + 1;
4651
ComplexMatrix m_band (ldm, nc);
4652
Complex *tmp_data = m_band.fortran_vec ();
4654
if (! mattype.is_dense ())
4656
octave_idx_type ii = 0;
4658
for (octave_idx_type j = 0; j < ldm; j++)
4659
for (octave_idx_type i = 0; i < nc; i++)
4660
tmp_data[ii++] = 0.;
4663
for (octave_idx_type j = 0; j < nc; j++)
4664
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4666
octave_idx_type ri = ridx (i);
4668
m_band(ri - j, j) = data (i);
4671
// Calculate the norm of the matrix, for later use.
4674
anorm = m_band.abs ().sum ().row (0).max ();
4677
F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4678
nr, n_lower, tmp_data, ldm, err
4679
F77_CHAR_ARG_LEN (1)));
4684
mattype.mark_as_unsymmetric ();
4685
typ = MatrixType::Banded;
4692
Array<Complex> z (dim_vector (2 * nr, 1));
4693
Complex *pz = z.fortran_vec ();
4694
Array<double> iz (dim_vector (nr, 1));
4695
double *piz = iz.fortran_vec ();
4697
F77_XFCN (zpbcon, ZPBCON,
4698
(F77_CONST_CHAR_ARG2 (&job, 1),
4699
nr, n_lower, tmp_data, ldm,
4700
anorm, rcond, pz, piz, err
4701
F77_CHAR_ARG_LEN (1)));
4706
volatile double rcond_plus_one = rcond + 1.0;
4708
if (rcond_plus_one == 1.0 || xisnan (rcond))
4714
sing_handler (rcond);
4715
mattype.mark_as_rectangular ();
4718
(*current_liboctave_error_handler)
4719
("matrix singular to machine precision, rcond = %g",
4728
octave_idx_type b_nr = b.rows ();
4729
octave_idx_type b_nc = b.cols ();
4730
OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
4732
// Take a first guess that the number of non-zero terms
4733
// will be as many as in b
4734
volatile octave_idx_type x_nz = b.nnz ();
4735
volatile octave_idx_type ii = 0;
4736
retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
4738
retval.xcidx (0) = 0;
4739
for (volatile octave_idx_type j = 0; j < b_nc; j++)
4741
for (octave_idx_type i = 0; i < b_nr; i++)
4742
Bx[i] = b.elem (i, j);
4744
F77_XFCN (zpbtrs, ZPBTRS,
4745
(F77_CONST_CHAR_ARG2 (&job, 1),
4746
nr, n_lower, 1, tmp_data,
4748
F77_CHAR_ARG_LEN (1)));
4752
(*current_liboctave_error_handler)
4753
("SparseComplexMatrix::solve solve failed");
4758
for (octave_idx_type i = 0; i < b_nr; i++)
4760
Complex tmp = Bx[i];
4765
// Resize the sparse matrix
4766
octave_idx_type sz = x_nz *
4768
sz = (sz > 10 ? sz : 10) + x_nz;
4769
retval.change_capacity (sz);
4772
retval.xdata (ii) = tmp;
4773
retval.xridx (ii++) = i;
4776
retval.xcidx (j+1) = ii;
4779
retval.maybe_compress ();
4784
if (typ == MatrixType::Banded)
4786
// Create the storage for the banded form of the sparse matrix
4787
octave_idx_type n_upper = mattype.nupper ();
4788
octave_idx_type n_lower = mattype.nlower ();
4789
octave_idx_type ldm = n_upper + 2 * n_lower + 1;
4791
ComplexMatrix m_band (ldm, nc);
4792
Complex *tmp_data = m_band.fortran_vec ();
4794
if (! mattype.is_dense ())
4796
octave_idx_type ii = 0;
4798
for (octave_idx_type j = 0; j < ldm; j++)
4799
for (octave_idx_type i = 0; i < nc; i++)
4800
tmp_data[ii++] = 0.;
4803
for (octave_idx_type j = 0; j < nc; j++)
4804
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4805
m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
4807
// Calculate the norm of the matrix, for later use.
4811
for (octave_idx_type j = 0; j < nr; j++)
4814
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4815
atmp += std::abs (data (i));
4821
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
4822
octave_idx_type *pipvt = ipvt.fortran_vec ();
4824
F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
4834
sing_handler (rcond);
4835
mattype.mark_as_rectangular ();
4838
(*current_liboctave_error_handler)
4839
("matrix singular to machine precision");
4847
Array<Complex> z (dim_vector (2 * nr, 1));
4848
Complex *pz = z.fortran_vec ();
4849
Array<double> iz (dim_vector (nr, 1));
4850
double *piz = iz.fortran_vec ();
4852
F77_XFCN (zgbcon, ZGBCON,
4853
(F77_CONST_CHAR_ARG2 (&job, 1),
4854
nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4855
anorm, rcond, pz, piz, err
4856
F77_CHAR_ARG_LEN (1)));
4861
volatile double rcond_plus_one = rcond + 1.0;
4863
if (rcond_plus_one == 1.0 || xisnan (rcond))
4869
sing_handler (rcond);
4870
mattype.mark_as_rectangular ();
4873
(*current_liboctave_error_handler)
4874
("matrix singular to machine precision, rcond = %g",
4884
volatile octave_idx_type x_nz = b.nnz ();
4885
octave_idx_type b_nc = b.cols ();
4886
retval = SparseComplexMatrix (nr, b_nc, x_nz);
4887
retval.xcidx (0) = 0;
4888
volatile octave_idx_type ii = 0;
4890
OCTAVE_LOCAL_BUFFER (Complex, work, nr);
4892
for (volatile octave_idx_type j = 0; j < b_nc; j++)
4894
for (octave_idx_type i = 0; i < nr; i++)
4896
for (octave_idx_type i = b.cidx (j);
4897
i < b.cidx (j+1); i++)
4898
work[b.ridx (i)] = b.data (i);
4900
F77_XFCN (zgbtrs, ZGBTRS,
4901
(F77_CONST_CHAR_ARG2 (&job, 1),
4902
nr, n_lower, n_upper, 1, tmp_data,
4903
ldm, pipvt, work, b.rows (), err
4904
F77_CHAR_ARG_LEN (1)));
4906
// Count non-zeros in work vector and adjust
4907
// space in retval if needed
4908
octave_idx_type new_nnz = 0;
4909
for (octave_idx_type i = 0; i < nr; i++)
4913
if (ii + new_nnz > x_nz)
4915
// Resize the sparse matrix
4916
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
4917
retval.change_capacity (sz);
4921
for (octave_idx_type i = 0; i < nr; i++)
4924
retval.xridx (ii) = i;
4925
retval.xdata (ii++) = work[i];
4927
retval.xcidx (j+1) = ii;
4930
retval.maybe_compress ();
4934
else if (typ != MatrixType::Banded_Hermitian)
4935
(*current_liboctave_error_handler) ("incorrect matrix type");
4942
SparseComplexMatrix::bsolve (MatrixType &mattype, const ComplexMatrix& b,
4943
octave_idx_type& err, double& rcond,
4944
solve_singularity_handler sing_handler,
4945
bool calc_cond) const
4947
ComplexMatrix retval;
4949
octave_idx_type nr = rows ();
4950
octave_idx_type nc = cols ();
4953
if (nr != nc || nr != b.rows ())
4954
(*current_liboctave_error_handler)
4955
("matrix dimension mismatch solution of linear equations");
4956
else if (nr == 0 || b.cols () == 0)
4957
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4960
// Print spparms("spumoni") info if requested
4961
volatile int typ = mattype.type ();
4964
if (typ == MatrixType::Banded_Hermitian)
4966
octave_idx_type n_lower = mattype.nlower ();
4967
octave_idx_type ldm = n_lower + 1;
4969
ComplexMatrix m_band (ldm, nc);
4970
Complex *tmp_data = m_band.fortran_vec ();
4972
if (! mattype.is_dense ())
4974
octave_idx_type ii = 0;
4976
for (octave_idx_type j = 0; j < ldm; j++)
4977
for (octave_idx_type i = 0; i < nc; i++)
4978
tmp_data[ii++] = 0.;
4981
for (octave_idx_type j = 0; j < nc; j++)
4982
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
4984
octave_idx_type ri = ridx (i);
4986
m_band(ri - j, j) = data (i);
4989
// Calculate the norm of the matrix, for later use.
4992
anorm = m_band.abs ().sum ().row (0).max ();
4995
F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
4996
nr, n_lower, tmp_data, ldm, err
4997
F77_CHAR_ARG_LEN (1)));
5001
// Matrix is not positive definite!! Fall through to
5002
// unsymmetric banded solver.
5004
mattype.mark_as_unsymmetric ();
5005
typ = MatrixType::Banded;
5012
Array<Complex> z (dim_vector (2 * nr, 1));
5013
Complex *pz = z.fortran_vec ();
5014
Array<double> iz (dim_vector (nr, 1));
5015
double *piz = iz.fortran_vec ();
5017
F77_XFCN (zpbcon, ZPBCON,
5018
(F77_CONST_CHAR_ARG2 (&job, 1),
5019
nr, n_lower, tmp_data, ldm,
5020
anorm, rcond, pz, piz, err
5021
F77_CHAR_ARG_LEN (1)));
5026
volatile double rcond_plus_one = rcond + 1.0;
5028
if (rcond_plus_one == 1.0 || xisnan (rcond))
5034
sing_handler (rcond);
5035
mattype.mark_as_rectangular ();
5038
(*current_liboctave_error_handler)
5039
("matrix singular to machine precision, rcond = %g",
5048
octave_idx_type b_nr = b.rows ();
5049
octave_idx_type b_nc = b.cols ();
5050
retval = ComplexMatrix (b);
5051
Complex *result = retval.fortran_vec ();
5053
F77_XFCN (zpbtrs, ZPBTRS,
5054
(F77_CONST_CHAR_ARG2 (&job, 1),
5055
nr, n_lower, b_nc, tmp_data,
5056
ldm, result, b_nr, err
5057
F77_CHAR_ARG_LEN (1)));
5061
(*current_liboctave_error_handler)
5062
("SparseComplexMatrix::solve solve failed");
5069
if (typ == MatrixType::Banded)
5071
// Create the storage for the banded form of the sparse matrix
5072
octave_idx_type n_upper = mattype.nupper ();
5073
octave_idx_type n_lower = mattype.nlower ();
5074
octave_idx_type ldm = n_upper + 2 * n_lower + 1;
5076
ComplexMatrix m_band (ldm, nc);
5077
Complex *tmp_data = m_band.fortran_vec ();
5079
if (! mattype.is_dense ())
5081
octave_idx_type ii = 0;
5083
for (octave_idx_type j = 0; j < ldm; j++)
5084
for (octave_idx_type i = 0; i < nc; i++)
5085
tmp_data[ii++] = 0.;
5088
for (octave_idx_type j = 0; j < nc; j++)
5089
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5090
m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5092
// Calculate the norm of the matrix, for later use.
5096
for (octave_idx_type j = 0; j < nr; j++)
5099
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5100
atmp += std::abs (data (i));
5106
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
5107
octave_idx_type *pipvt = ipvt.fortran_vec ();
5109
F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5119
sing_handler (rcond);
5120
mattype.mark_as_rectangular ();
5123
(*current_liboctave_error_handler)
5124
("matrix singular to machine precision");
5131
Array<Complex> z (dim_vector (2 * nr, 1));
5132
Complex *pz = z.fortran_vec ();
5133
Array<double> iz (dim_vector (nr, 1));
5134
double *piz = iz.fortran_vec ();
5136
F77_XFCN (zgbcon, ZGBCON,
5137
(F77_CONST_CHAR_ARG2 (&job, 1),
5138
nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5139
anorm, rcond, pz, piz, err
5140
F77_CHAR_ARG_LEN (1)));
5145
volatile double rcond_plus_one = rcond + 1.0;
5147
if (rcond_plus_one == 1.0 || xisnan (rcond))
5153
sing_handler (rcond);
5154
mattype.mark_as_rectangular ();
5157
(*current_liboctave_error_handler)
5158
("matrix singular to machine precision, rcond = %g",
5168
octave_idx_type b_nc = b.cols ();
5169
retval = ComplexMatrix (b);
5170
Complex *result = retval.fortran_vec ();
5172
F77_XFCN (zgbtrs, ZGBTRS,
5173
(F77_CONST_CHAR_ARG2 (&job, 1),
5174
nr, n_lower, n_upper, b_nc, tmp_data,
5175
ldm, pipvt, result, b.rows (), err
5176
F77_CHAR_ARG_LEN (1)));
5180
else if (typ != MatrixType::Banded_Hermitian)
5181
(*current_liboctave_error_handler) ("incorrect matrix type");
5188
SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseComplexMatrix& b,
5189
octave_idx_type& err, double& rcond,
5190
solve_singularity_handler sing_handler,
5191
bool calc_cond) const
5193
SparseComplexMatrix retval;
5195
octave_idx_type nr = rows ();
5196
octave_idx_type nc = cols ();
5199
if (nr != nc || nr != b.rows ())
5200
(*current_liboctave_error_handler)
5201
("matrix dimension mismatch solution of linear equations");
5202
else if (nr == 0 || b.cols () == 0)
5203
retval = SparseComplexMatrix (nc, b.cols ());
5206
// Print spparms("spumoni") info if requested
5207
volatile int typ = mattype.type ();
5210
if (typ == MatrixType::Banded_Hermitian)
5212
octave_idx_type n_lower = mattype.nlower ();
5213
octave_idx_type ldm = n_lower + 1;
5215
ComplexMatrix m_band (ldm, nc);
5216
Complex *tmp_data = m_band.fortran_vec ();
5218
if (! mattype.is_dense ())
5220
octave_idx_type ii = 0;
5222
for (octave_idx_type j = 0; j < ldm; j++)
5223
for (octave_idx_type i = 0; i < nc; i++)
5224
tmp_data[ii++] = 0.;
5227
for (octave_idx_type j = 0; j < nc; j++)
5228
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5230
octave_idx_type ri = ridx (i);
5232
m_band(ri - j, j) = data (i);
5235
// Calculate the norm of the matrix, for later use.
5238
anorm = m_band.abs ().sum ().row (0).max ();
5241
F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
5242
nr, n_lower, tmp_data, ldm, err
5243
F77_CHAR_ARG_LEN (1)));
5247
// Matrix is not positive definite!! Fall through to
5248
// unsymmetric banded solver.
5249
mattype.mark_as_unsymmetric ();
5250
typ = MatrixType::Banded;
5259
Array<Complex> z (dim_vector (2 * nr, 1));
5260
Complex *pz = z.fortran_vec ();
5261
Array<double> iz (dim_vector (nr, 1));
5262
double *piz = iz.fortran_vec ();
5264
F77_XFCN (zpbcon, ZPBCON,
5265
(F77_CONST_CHAR_ARG2 (&job, 1),
5266
nr, n_lower, tmp_data, ldm,
5267
anorm, rcond, pz, piz, err
5268
F77_CHAR_ARG_LEN (1)));
5273
volatile double rcond_plus_one = rcond + 1.0;
5275
if (rcond_plus_one == 1.0 || xisnan (rcond))
5281
sing_handler (rcond);
5282
mattype.mark_as_rectangular ();
5285
(*current_liboctave_error_handler)
5286
("matrix singular to machine precision, rcond = %g",
5295
octave_idx_type b_nr = b.rows ();
5296
octave_idx_type b_nc = b.cols ();
5297
OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
5299
// Take a first guess that the number of non-zero terms
5300
// will be as many as in b
5301
volatile octave_idx_type x_nz = b.nnz ();
5302
volatile octave_idx_type ii = 0;
5303
retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
5305
retval.xcidx (0) = 0;
5306
for (volatile octave_idx_type j = 0; j < b_nc; j++)
5309
for (octave_idx_type i = 0; i < b_nr; i++)
5312
F77_XFCN (zpbtrs, ZPBTRS,
5313
(F77_CONST_CHAR_ARG2 (&job, 1),
5314
nr, n_lower, 1, tmp_data,
5316
F77_CHAR_ARG_LEN (1)));
5320
(*current_liboctave_error_handler)
5321
("SparseMatrix::solve solve failed");
5326
// Count non-zeros in work vector and adjust
5327
// space in retval if needed
5328
octave_idx_type new_nnz = 0;
5329
for (octave_idx_type i = 0; i < nr; i++)
5333
if (ii + new_nnz > x_nz)
5335
// Resize the sparse matrix
5336
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5337
retval.change_capacity (sz);
5341
for (octave_idx_type i = 0; i < nr; i++)
5344
retval.xridx (ii) = i;
5345
retval.xdata (ii++) = Bx[i];
5348
retval.xcidx (j+1) = ii;
5351
retval.maybe_compress ();
5356
if (typ == MatrixType::Banded)
5358
// Create the storage for the banded form of the sparse matrix
5359
octave_idx_type n_upper = mattype.nupper ();
5360
octave_idx_type n_lower = mattype.nlower ();
5361
octave_idx_type ldm = n_upper + 2 * n_lower + 1;
5363
ComplexMatrix m_band (ldm, nc);
5364
Complex *tmp_data = m_band.fortran_vec ();
5366
if (! mattype.is_dense ())
5368
octave_idx_type ii = 0;
5370
for (octave_idx_type j = 0; j < ldm; j++)
5371
for (octave_idx_type i = 0; i < nc; i++)
5372
tmp_data[ii++] = 0.;
5375
for (octave_idx_type j = 0; j < nc; j++)
5376
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5377
m_band(ridx (i) - j + n_lower + n_upper, j) = data (i);
5379
// Calculate the norm of the matrix, for later use.
5383
for (octave_idx_type j = 0; j < nr; j++)
5386
for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
5387
atmp += std::abs (data (i));
5393
Array<octave_idx_type> ipvt (dim_vector (nr, 1));
5394
octave_idx_type *pipvt = ipvt.fortran_vec ();
5396
F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
5406
sing_handler (rcond);
5407
mattype.mark_as_rectangular ();
5410
(*current_liboctave_error_handler)
5411
("matrix singular to machine precision");
5419
Array<Complex> z (dim_vector (2 * nr, 1));
5420
Complex *pz = z.fortran_vec ();
5421
Array<double> iz (dim_vector (nr, 1));
5422
double *piz = iz.fortran_vec ();
5424
F77_XFCN (zgbcon, ZGBCON,
5425
(F77_CONST_CHAR_ARG2 (&job, 1),
5426
nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5427
anorm, rcond, pz, piz, err
5428
F77_CHAR_ARG_LEN (1)));
5433
volatile double rcond_plus_one = rcond + 1.0;
5435
if (rcond_plus_one == 1.0 || xisnan (rcond))
5441
sing_handler (rcond);
5442
mattype.mark_as_rectangular ();
5445
(*current_liboctave_error_handler)
5446
("matrix singular to machine precision, rcond = %g",
5456
volatile octave_idx_type x_nz = b.nnz ();
5457
octave_idx_type b_nc = b.cols ();
5458
retval = SparseComplexMatrix (nr, b_nc, x_nz);
5459
retval.xcidx (0) = 0;
5460
volatile octave_idx_type ii = 0;
5462
OCTAVE_LOCAL_BUFFER (Complex, Bx, nr);
5464
for (volatile octave_idx_type j = 0; j < b_nc; j++)
5466
for (octave_idx_type i = 0; i < nr; i++)
5469
for (octave_idx_type i = b.cidx (j);
5470
i < b.cidx (j+1); i++)
5471
Bx[b.ridx (i)] = b.data (i);
5473
F77_XFCN (zgbtrs, ZGBTRS,
5474
(F77_CONST_CHAR_ARG2 (&job, 1),
5475
nr, n_lower, n_upper, 1, tmp_data,
5476
ldm, pipvt, Bx, b.rows (), err
5477
F77_CHAR_ARG_LEN (1)));
5479
// Count non-zeros in work vector and adjust
5480
// space in retval if needed
5481
octave_idx_type new_nnz = 0;
5482
for (octave_idx_type i = 0; i < nr; i++)
5486
if (ii + new_nnz > x_nz)
5488
// Resize the sparse matrix
5489
octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
5490
retval.change_capacity (sz);
5494
for (octave_idx_type i = 0; i < nr; i++)
5497
retval.xridx (ii) = i;
5498
retval.xdata (ii++) = Bx[i];
5500
retval.xcidx (j+1) = ii;
5503
retval.maybe_compress ();
5507
else if (typ != MatrixType::Banded_Hermitian)
5508
(*current_liboctave_error_handler) ("incorrect matrix type");
5515
SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond,
5516
Matrix &Control, Matrix &Info,
5517
solve_singularity_handler sing_handler,
5518
bool calc_cond) const
5520
// The return values
5525
// Setup the control parameters
5526
Control = Matrix (UMFPACK_CONTROL, 1);
5527
double *control = Control.fortran_vec ();
5528
UMFPACK_ZNAME (defaults) (control);
5530
double tmp = octave_sparse_params::get_key ("spumoni");
5532
Control (UMFPACK_PRL) = tmp;
5533
tmp = octave_sparse_params::get_key ("piv_tol");
5536
Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp;
5537
Control (UMFPACK_PIVOT_TOLERANCE) = tmp;
5540
// Set whether we are allowed to modify Q or not
5541
tmp = octave_sparse_params::get_key ("autoamd");
5543
Control (UMFPACK_FIXQ) = tmp;
5545
UMFPACK_ZNAME (report_control) (control);
5547
const octave_idx_type *Ap = cidx ();
5548
const octave_idx_type *Ai = ridx ();
5549
const Complex *Ax = data ();
5550
octave_idx_type nr = rows ();
5551
octave_idx_type nc = cols ();
5553
UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai,
5554
reinterpret_cast<const double *> (Ax),
5558
Info = Matrix (1, UMFPACK_INFO);
5559
double *info = Info.fortran_vec ();
5560
int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
5561
reinterpret_cast<const double *> (Ax),
5562
0, 0, &Symbolic, control, info);
5566
(*current_liboctave_error_handler)
5567
("SparseComplexMatrix::solve symbolic factorization failed");
5570
UMFPACK_ZNAME (report_status) (control, status);
5571
UMFPACK_ZNAME (report_info) (control, info);
5573
UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5577
UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
5579
status = UMFPACK_ZNAME (numeric) (Ap, Ai,
5580
reinterpret_cast<const double *> (Ax),
5581
0, Symbolic, &Numeric, control, info);
5582
UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5585
rcond = Info (UMFPACK_RCOND);
5588
volatile double rcond_plus_one = rcond + 1.0;
5590
if (status == UMFPACK_WARNING_singular_matrix ||
5591
rcond_plus_one == 1.0 || xisnan (rcond))
5593
UMFPACK_ZNAME (report_numeric) (Numeric, control);
5598
sing_handler (rcond);
5600
(*current_liboctave_error_handler)
5601
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
5605
else if (status < 0)
5607
(*current_liboctave_error_handler)
5608
("SparseComplexMatrix::solve numeric factorization failed");
5610
UMFPACK_ZNAME (report_status) (control, status);
5611
UMFPACK_ZNAME (report_info) (control, info);
5617
UMFPACK_ZNAME (report_numeric) (Numeric, control);
5622
UMFPACK_ZNAME (free_numeric) (&Numeric);
5624
(*current_liboctave_error_handler) ("UMFPACK not installed");
5631
SparseComplexMatrix::fsolve (MatrixType &mattype, const Matrix& b,
5632
octave_idx_type& err, double& rcond,
5633
solve_singularity_handler sing_handler,
5634
bool calc_cond) const
5636
ComplexMatrix retval;
5638
octave_idx_type nr = rows ();
5639
octave_idx_type nc = cols ();
5642
if (nr != nc || nr != b.rows ())
5643
(*current_liboctave_error_handler)
5644
("matrix dimension mismatch solution of linear equations");
5645
else if (nr == 0 || b.cols () == 0)
5646
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
5649
// Print spparms("spumoni") info if requested
5650
volatile int typ = mattype.type ();
5653
if (typ == MatrixType::Hermitian)
5656
cholmod_common Common;
5657
cholmod_common *cm = &Common;
5659
// Setup initial parameters
5660
CHOLMOD_NAME(start) (cm);
5661
cm->prefer_zomplex = false;
5663
double spu = octave_sparse_params::get_key ("spumoni");
5667
cm->print_function = 0;
5671
cm->print = static_cast<int> (spu) + 2;
5672
cm->print_function =&SparseCholPrint;
5675
cm->error_handler = &SparseCholError;
5676
cm->complex_divide = CHOLMOD_NAME(divcomplex);
5677
cm->hypotenuse = CHOLMOD_NAME(hypot);
5679
cm->final_ll = true;
5681
cholmod_sparse Astore;
5682
cholmod_sparse *A = &Astore;
5693
#ifdef USE_64_BIT_IDX_T
5694
A->itype = CHOLMOD_LONG;
5696
A->itype = CHOLMOD_INT;
5698
A->dtype = CHOLMOD_DOUBLE;
5700
A->xtype = CHOLMOD_COMPLEX;
5707
cholmod_dense Bstore;
5708
cholmod_dense *B = &Bstore;
5709
B->nrow = b.rows ();
5710
B->ncol = b.cols ();
5712
B->nzmax = B->nrow * B->ncol;
5713
B->dtype = CHOLMOD_DOUBLE;
5714
B->xtype = CHOLMOD_REAL;
5715
if (nc < 1 || b.cols () < 1)
5718
// We won't alter it, honest :-)
5719
B->x = const_cast<double *>(b.fortran_vec ());
5722
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5723
L = CHOLMOD_NAME(analyze) (A, cm);
5724
CHOLMOD_NAME(factorize) (A, L, cm);
5726
rcond = CHOLMOD_NAME(rcond)(L, cm);
5729
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5733
// Either its indefinite or singular. Try UMFPACK
5734
mattype.mark_as_unsymmetric ();
5735
typ = MatrixType::Full;
5739
volatile double rcond_plus_one = rcond + 1.0;
5741
if (rcond_plus_one == 1.0 || xisnan (rcond))
5747
sing_handler (rcond);
5748
mattype.mark_as_rectangular ();
5751
(*current_liboctave_error_handler)
5752
("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5759
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5760
X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5761
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5763
retval.resize (b.rows (), b.cols ());
5764
for (octave_idx_type j = 0; j < b.cols (); j++)
5766
octave_idx_type jr = j * b.rows ();
5767
for (octave_idx_type i = 0; i < b.rows (); i++)
5768
retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i];
5771
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5772
CHOLMOD_NAME(free_dense) (&X, cm);
5773
CHOLMOD_NAME(free_factor) (&L, cm);
5774
CHOLMOD_NAME(finish) (cm);
5775
static char tmp[] = " ";
5776
CHOLMOD_NAME(print_common) (tmp, cm);
5777
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5780
(*current_liboctave_warning_handler)
5781
("CHOLMOD not installed");
5783
mattype.mark_as_unsymmetric ();
5784
typ = MatrixType::Full;
5788
if (typ == MatrixType::Full)
5791
Matrix Control, Info;
5792
void *Numeric = factorize (err, rcond, Control, Info,
5793
sing_handler, calc_cond);
5797
octave_idx_type b_nr = b.rows ();
5798
octave_idx_type b_nc = b.cols ();
5800
double *control = Control.fortran_vec ();
5801
double *info = Info.fortran_vec ();
5802
const octave_idx_type *Ap = cidx ();
5803
const octave_idx_type *Ai = ridx ();
5804
const Complex *Ax = data ();
5805
#ifdef UMFPACK_SEPARATE_SPLIT
5806
const double *Bx = b.fortran_vec ();
5807
OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
5808
for (octave_idx_type i = 0; i < b_nr; i++)
5811
OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
5813
retval.resize (b_nr, b_nc);
5814
Complex *Xx = retval.fortran_vec ();
5816
for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
5818
#ifdef UMFPACK_SEPARATE_SPLIT
5819
status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5821
reinterpret_cast<const double *> (Ax),
5823
reinterpret_cast<double *> (&Xx[iidx]),
5825
&Bx[iidx], Bz, Numeric,
5828
for (octave_idx_type i = 0; i < b_nr; i++)
5829
Bz[i] = b.elem (i, j);
5831
status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5833
reinterpret_cast<const double *> (Ax),
5835
reinterpret_cast<double *> (&Xx[iidx]),
5837
reinterpret_cast<const double *> (Bz),
5844
(*current_liboctave_error_handler)
5845
("SparseComplexMatrix::solve solve failed");
5847
UMFPACK_ZNAME (report_status) (control, status);
5855
UMFPACK_ZNAME (report_info) (control, info);
5857
UMFPACK_ZNAME (free_numeric) (&Numeric);
5860
mattype.mark_as_rectangular ();
5863
(*current_liboctave_error_handler) ("UMFPACK not installed");
5866
else if (typ != MatrixType::Hermitian)
5867
(*current_liboctave_error_handler) ("incorrect matrix type");
5874
SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseMatrix& b,
5875
octave_idx_type& err, double& rcond,
5876
solve_singularity_handler sing_handler,
5877
bool calc_cond) const
5879
SparseComplexMatrix retval;
5881
octave_idx_type nr = rows ();
5882
octave_idx_type nc = cols ();
5885
if (nr != nc || nr != b.rows ())
5886
(*current_liboctave_error_handler)
5887
("matrix dimension mismatch solution of linear equations");
5888
else if (nr == 0 || b.cols () == 0)
5889
retval = SparseComplexMatrix (nc, b.cols ());
5892
// Print spparms("spumoni") info if requested
5893
volatile int typ = mattype.type ();
5896
if (typ == MatrixType::Hermitian)
5899
cholmod_common Common;
5900
cholmod_common *cm = &Common;
5902
// Setup initial parameters
5903
CHOLMOD_NAME(start) (cm);
5904
cm->prefer_zomplex = false;
5906
double spu = octave_sparse_params::get_key ("spumoni");
5910
cm->print_function = 0;
5914
cm->print = static_cast<int> (spu) + 2;
5915
cm->print_function =&SparseCholPrint;
5918
cm->error_handler = &SparseCholError;
5919
cm->complex_divide = CHOLMOD_NAME(divcomplex);
5920
cm->hypotenuse = CHOLMOD_NAME(hypot);
5922
cm->final_ll = true;
5924
cholmod_sparse Astore;
5925
cholmod_sparse *A = &Astore;
5936
#ifdef USE_64_BIT_IDX_T
5937
A->itype = CHOLMOD_LONG;
5939
A->itype = CHOLMOD_INT;
5941
A->dtype = CHOLMOD_DOUBLE;
5943
A->xtype = CHOLMOD_COMPLEX;
5950
cholmod_sparse Bstore;
5951
cholmod_sparse *B = &Bstore;
5952
B->nrow = b.rows ();
5953
B->ncol = b.cols ();
5956
B->nzmax = b.nnz ();
5960
#ifdef USE_64_BIT_IDX_T
5961
B->itype = CHOLMOD_LONG;
5963
B->itype = CHOLMOD_INT;
5965
B->dtype = CHOLMOD_DOUBLE;
5967
B->xtype = CHOLMOD_REAL;
5969
if (b.rows () < 1 || b.cols () < 1)
5975
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5976
L = CHOLMOD_NAME(analyze) (A, cm);
5977
CHOLMOD_NAME(factorize) (A, L, cm);
5979
rcond = CHOLMOD_NAME(rcond)(L, cm);
5982
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5986
// Either its indefinite or singular. Try UMFPACK
5987
mattype.mark_as_unsymmetric ();
5988
typ = MatrixType::Full;
5992
volatile double rcond_plus_one = rcond + 1.0;
5994
if (rcond_plus_one == 1.0 || xisnan (rcond))
6000
sing_handler (rcond);
6001
mattype.mark_as_rectangular ();
6004
(*current_liboctave_error_handler)
6005
("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
6012
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6013
X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6014
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6016
retval = SparseComplexMatrix
6017
(static_cast<octave_idx_type>(X->nrow),
6018
static_cast<octave_idx_type>(X->ncol),
6019
static_cast<octave_idx_type>(X->nzmax));
6020
for (octave_idx_type j = 0;
6021
j <= static_cast<octave_idx_type>(X->ncol); j++)
6022
retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6023
for (octave_idx_type j = 0;
6024
j < static_cast<octave_idx_type>(X->nzmax); j++)
6026
retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6027
retval.xdata (j) = static_cast<Complex *>(X->x)[j];
6030
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6031
CHOLMOD_NAME(free_sparse) (&X, cm);
6032
CHOLMOD_NAME(free_factor) (&L, cm);
6033
CHOLMOD_NAME(finish) (cm);
6034
static char tmp[] = " ";
6035
CHOLMOD_NAME(print_common) (tmp, cm);
6036
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6039
(*current_liboctave_warning_handler)
6040
("CHOLMOD not installed");
6042
mattype.mark_as_unsymmetric ();
6043
typ = MatrixType::Full;
6047
if (typ == MatrixType::Full)
6050
Matrix Control, Info;
6051
void *Numeric = factorize (err, rcond, Control, Info,
6052
sing_handler, calc_cond);
6056
octave_idx_type b_nr = b.rows ();
6057
octave_idx_type b_nc = b.cols ();
6059
double *control = Control.fortran_vec ();
6060
double *info = Info.fortran_vec ();
6061
const octave_idx_type *Ap = cidx ();
6062
const octave_idx_type *Ai = ridx ();
6063
const Complex *Ax = data ();
6065
#ifdef UMFPACK_SEPARATE_SPLIT
6066
OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
6067
OCTAVE_LOCAL_BUFFER (double, Bz, b_nr);
6068
for (octave_idx_type i = 0; i < b_nr; i++)
6071
OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr);
6074
// Take a first guess that the number of non-zero terms
6075
// will be as many as in b
6076
octave_idx_type x_nz = b.nnz ();
6077
octave_idx_type ii = 0;
6078
retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6080
OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
6082
retval.xcidx (0) = 0;
6083
for (octave_idx_type j = 0; j < b_nc; j++)
6086
#ifdef UMFPACK_SEPARATE_SPLIT
6087
for (octave_idx_type i = 0; i < b_nr; i++)
6088
Bx[i] = b.elem (i, j);
6090
status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
6092
reinterpret_cast<const double *> (Ax),
6094
reinterpret_cast<double *> (Xx),
6096
Bx, Bz, Numeric, control,
6099
for (octave_idx_type i = 0; i < b_nr; i++)
6100
Bz[i] = b.elem (i, j);
6102
status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
6103
reinterpret_cast<const double *> (Ax),
6105
reinterpret_cast<double *> (Xx),
6107
reinterpret_cast<double *> (Bz),
6114
(*current_liboctave_error_handler)
6115
("SparseComplexMatrix::solve solve failed");
6117
UMFPACK_ZNAME (report_status) (control, status);
6124
for (octave_idx_type i = 0; i < b_nr; i++)
6126
Complex tmp = Xx[i];
6131
// Resize the sparse matrix
6132
octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
6133
sz = (sz > 10 ? sz : 10) + x_nz;
6134
retval.change_capacity (sz);
6137
retval.xdata (ii) = tmp;
6138
retval.xridx (ii++) = i;
6141
retval.xcidx (j+1) = ii;
6144
retval.maybe_compress ();
6146
UMFPACK_ZNAME (report_info) (control, info);
6148
UMFPACK_ZNAME (free_numeric) (&Numeric);
6151
mattype.mark_as_rectangular ();
6154
(*current_liboctave_error_handler) ("UMFPACK not installed");
6157
else if (typ != MatrixType::Hermitian)
6158
(*current_liboctave_error_handler) ("incorrect matrix type");
6165
SparseComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
6166
octave_idx_type& err, double& rcond,
6167
solve_singularity_handler sing_handler,
6168
bool calc_cond) const
6170
ComplexMatrix retval;
6172
octave_idx_type nr = rows ();
6173
octave_idx_type nc = cols ();
6176
if (nr != nc || nr != b.rows ())
6177
(*current_liboctave_error_handler)
6178
("matrix dimension mismatch solution of linear equations");
6179
else if (nr == 0 || b.cols () == 0)
6180
retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
6183
// Print spparms("spumoni") info if requested
6184
volatile int typ = mattype.type ();
6187
if (typ == MatrixType::Hermitian)
6190
cholmod_common Common;
6191
cholmod_common *cm = &Common;
6193
// Setup initial parameters
6194
CHOLMOD_NAME(start) (cm);
6195
cm->prefer_zomplex = false;
6197
double spu = octave_sparse_params::get_key ("spumoni");
6201
cm->print_function = 0;
6205
cm->print = static_cast<int> (spu) + 2;
6206
cm->print_function =&SparseCholPrint;
6209
cm->error_handler = &SparseCholError;
6210
cm->complex_divide = CHOLMOD_NAME(divcomplex);
6211
cm->hypotenuse = CHOLMOD_NAME(hypot);
6213
cm->final_ll = true;
6215
cholmod_sparse Astore;
6216
cholmod_sparse *A = &Astore;
6227
#ifdef USE_64_BIT_IDX_T
6228
A->itype = CHOLMOD_LONG;
6230
A->itype = CHOLMOD_INT;
6232
A->dtype = CHOLMOD_DOUBLE;
6234
A->xtype = CHOLMOD_COMPLEX;
6241
cholmod_dense Bstore;
6242
cholmod_dense *B = &Bstore;
6243
B->nrow = b.rows ();
6244
B->ncol = b.cols ();
6246
B->nzmax = B->nrow * B->ncol;
6247
B->dtype = CHOLMOD_DOUBLE;
6248
B->xtype = CHOLMOD_COMPLEX;
6249
if (nc < 1 || b.cols () < 1)
6252
// We won't alter it, honest :-)
6253
B->x = const_cast<Complex *>(b.fortran_vec ());
6256
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6257
L = CHOLMOD_NAME(analyze) (A, cm);
6258
CHOLMOD_NAME(factorize) (A, L, cm);
6260
rcond = CHOLMOD_NAME(rcond)(L, cm);
6263
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6267
// Either its indefinite or singular. Try UMFPACK
6268
mattype.mark_as_unsymmetric ();
6269
typ = MatrixType::Full;
6273
volatile double rcond_plus_one = rcond + 1.0;
6275
if (rcond_plus_one == 1.0 || xisnan (rcond))
6281
sing_handler (rcond);
6282
mattype.mark_as_rectangular ();
6285
(*current_liboctave_error_handler)
6286
("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
6293
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6294
X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
6295
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6297
retval.resize (b.rows (), b.cols ());
6298
for (octave_idx_type j = 0; j < b.cols (); j++)
6300
octave_idx_type jr = j * b.rows ();
6301
for (octave_idx_type i = 0; i < b.rows (); i++)
6302
retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i];
6305
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6306
CHOLMOD_NAME(free_dense) (&X, cm);
6307
CHOLMOD_NAME(free_factor) (&L, cm);
6308
CHOLMOD_NAME(finish) (cm);
6309
static char tmp[] = " ";
6310
CHOLMOD_NAME(print_common) (tmp, cm);
6311
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6314
(*current_liboctave_warning_handler)
6315
("CHOLMOD not installed");
6317
mattype.mark_as_unsymmetric ();
6318
typ = MatrixType::Full;
6322
if (typ == MatrixType::Full)
6325
Matrix Control, Info;
6326
void *Numeric = factorize (err, rcond, Control, Info,
6327
sing_handler, calc_cond);
6331
octave_idx_type b_nr = b.rows ();
6332
octave_idx_type b_nc = b.cols ();
6334
double *control = Control.fortran_vec ();
6335
double *info = Info.fortran_vec ();
6336
const octave_idx_type *Ap = cidx ();
6337
const octave_idx_type *Ai = ridx ();
6338
const Complex *Ax = data ();
6339
const Complex *Bx = b.fortran_vec ();
6341
retval.resize (b_nr, b_nc);
6342
Complex *Xx = retval.fortran_vec ();
6344
for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
6347
UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
6348
reinterpret_cast<const double *> (Ax),
6350
reinterpret_cast<double *> (&Xx[iidx]),
6352
reinterpret_cast<const double *> (&Bx[iidx]),
6353
0, Numeric, control, info);
6357
(*current_liboctave_error_handler)
6358
("SparseComplexMatrix::solve solve failed");
6360
UMFPACK_ZNAME (report_status) (control, status);
6368
UMFPACK_ZNAME (report_info) (control, info);
6370
UMFPACK_ZNAME (free_numeric) (&Numeric);
6373
mattype.mark_as_rectangular ();
6376
(*current_liboctave_error_handler) ("UMFPACK not installed");
6379
else if (typ != MatrixType::Hermitian)
6380
(*current_liboctave_error_handler) ("incorrect matrix type");
6387
SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseComplexMatrix& b,
6388
octave_idx_type& err, double& rcond,
6389
solve_singularity_handler sing_handler,
6390
bool calc_cond) const
6392
SparseComplexMatrix retval;
6394
octave_idx_type nr = rows ();
6395
octave_idx_type nc = cols ();
6398
if (nr != nc || nr != b.rows ())
6399
(*current_liboctave_error_handler)
6400
("matrix dimension mismatch solution of linear equations");
6401
else if (nr == 0 || b.cols () == 0)
6402
retval = SparseComplexMatrix (nc, b.cols ());
6405
// Print spparms("spumoni") info if requested
6406
volatile int typ = mattype.type ();
6409
if (typ == MatrixType::Hermitian)
6412
cholmod_common Common;
6413
cholmod_common *cm = &Common;
6415
// Setup initial parameters
6416
CHOLMOD_NAME(start) (cm);
6417
cm->prefer_zomplex = false;
6419
double spu = octave_sparse_params::get_key ("spumoni");
6423
cm->print_function = 0;
6427
cm->print = static_cast<int> (spu) + 2;
6428
cm->print_function =&SparseCholPrint;
6431
cm->error_handler = &SparseCholError;
6432
cm->complex_divide = CHOLMOD_NAME(divcomplex);
6433
cm->hypotenuse = CHOLMOD_NAME(hypot);
6435
cm->final_ll = true;
6437
cholmod_sparse Astore;
6438
cholmod_sparse *A = &Astore;
6449
#ifdef USE_64_BIT_IDX_T
6450
A->itype = CHOLMOD_LONG;
6452
A->itype = CHOLMOD_INT;
6454
A->dtype = CHOLMOD_DOUBLE;
6456
A->xtype = CHOLMOD_COMPLEX;
6463
cholmod_sparse Bstore;
6464
cholmod_sparse *B = &Bstore;
6465
B->nrow = b.rows ();
6466
B->ncol = b.cols ();
6469
B->nzmax = b.nnz ();
6473
#ifdef USE_64_BIT_IDX_T
6474
B->itype = CHOLMOD_LONG;
6476
B->itype = CHOLMOD_INT;
6478
B->dtype = CHOLMOD_DOUBLE;
6480
B->xtype = CHOLMOD_COMPLEX;
6482
if (b.rows () < 1 || b.cols () < 1)
6488
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6489
L = CHOLMOD_NAME(analyze) (A, cm);
6490
CHOLMOD_NAME(factorize) (A, L, cm);
6492
rcond = CHOLMOD_NAME(rcond)(L, cm);
6495
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6499
// Either its indefinite or singular. Try UMFPACK
6500
mattype.mark_as_unsymmetric ();
6501
typ = MatrixType::Full;
6505
volatile double rcond_plus_one = rcond + 1.0;
6507
if (rcond_plus_one == 1.0 || xisnan (rcond))
6513
sing_handler (rcond);
6514
mattype.mark_as_rectangular ();
6517
(*current_liboctave_error_handler)
6518
("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
6525
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6526
X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6527
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6529
retval = SparseComplexMatrix
6530
(static_cast<octave_idx_type>(X->nrow),
6531
static_cast<octave_idx_type>(X->ncol),
6532
static_cast<octave_idx_type>(X->nzmax));
6533
for (octave_idx_type j = 0;
6534
j <= static_cast<octave_idx_type>(X->ncol); j++)
6535
retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6536
for (octave_idx_type j = 0;
6537
j < static_cast<octave_idx_type>(X->nzmax); j++)
6539
retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j];
6540
retval.xdata (j) = static_cast<Complex *>(X->x)[j];
6543
BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6544
CHOLMOD_NAME(free_sparse) (&X, cm);
6545
CHOLMOD_NAME(free_factor) (&L, cm);
6546
CHOLMOD_NAME(finish) (cm);
6547
static char tmp[] = " ";
6548
CHOLMOD_NAME(print_common) (tmp, cm);
6549
END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6552
(*current_liboctave_warning_handler)
6553
("CHOLMOD not installed");
6555
mattype.mark_as_unsymmetric ();
6556
typ = MatrixType::Full;
6560
if (typ == MatrixType::Full)
6563
Matrix Control, Info;
6564
void *Numeric = factorize (err, rcond, Control, Info,
6565
sing_handler, calc_cond);
6569
octave_idx_type b_nr = b.rows ();
6570
octave_idx_type b_nc = b.cols ();
6572
double *control = Control.fortran_vec ();
6573
double *info = Info.fortran_vec ();
6574
const octave_idx_type *Ap = cidx ();
6575
const octave_idx_type *Ai = ridx ();
6576
const Complex *Ax = data ();
6578
OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr);
6580
// Take a first guess that the number of non-zero terms
6581
// will be as many as in b
6582
octave_idx_type x_nz = b.nnz ();
6583
octave_idx_type ii = 0;
6584
retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
6586
OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
6588
retval.xcidx (0) = 0;
6589
for (octave_idx_type j = 0; j < b_nc; j++)
6591
for (octave_idx_type i = 0; i < b_nr; i++)
6594
status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
6596
reinterpret_cast<const double *> (Ax),
6598
reinterpret_cast<double *> (Xx),
6600
reinterpret_cast<double *> (Bx),
6601
0, Numeric, control, info);
6605
(*current_liboctave_error_handler)
6606
("SparseComplexMatrix::solve solve failed");
6608
UMFPACK_ZNAME (report_status) (control, status);
6615
for (octave_idx_type i = 0; i < b_nr; i++)
6617
Complex tmp = Xx[i];
6622
// Resize the sparse matrix
6623
octave_idx_type sz = x_nz * (b_nc - j) / b_nc;
6624
sz = (sz > 10 ? sz : 10) + x_nz;
6625
retval.change_capacity (sz);
6628
retval.xdata (ii) = tmp;
6629
retval.xridx (ii++) = i;
6632
retval.xcidx (j+1) = ii;
6635
retval.maybe_compress ();
6637
rcond = Info (UMFPACK_RCOND);
6638
volatile double rcond_plus_one = rcond + 1.0;
6640
if (status == UMFPACK_WARNING_singular_matrix ||
6641
rcond_plus_one == 1.0 || xisnan (rcond))
6646
sing_handler (rcond);
6648
(*current_liboctave_error_handler)
6649
("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
6654
UMFPACK_ZNAME (report_info) (control, info);
6656
UMFPACK_ZNAME (free_numeric) (&Numeric);
6659
mattype.mark_as_rectangular ();
6662
(*current_liboctave_error_handler) ("UMFPACK not installed");
6665
else if (typ != MatrixType::Hermitian)
6666
(*current_liboctave_error_handler) ("incorrect matrix type");
6673
SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b) const
6675
octave_idx_type info;
6677
return solve (mattype, b, info, rcond, 0);
6681
SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
6682
octave_idx_type& info) const
6685
return solve (mattype, b, info, rcond, 0);
6689
SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
6690
octave_idx_type& info, double& rcond) const
6692
return solve (mattype, b, info, rcond, 0);
6696
SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b,
6697
octave_idx_type& err, double& rcond,
6698
solve_singularity_handler sing_handler,
6699
bool singular_fallback) const
6701
ComplexMatrix retval;
6702
int typ = mattype.type (false);
6704
if (typ == MatrixType::Unknown)
6705
typ = mattype.type (*this);
6707
if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
6708
retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6709
else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6710
retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6711
else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6712
retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6713
else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6714
retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6715
else if (typ == MatrixType::Tridiagonal ||
6716
typ == MatrixType::Tridiagonal_Hermitian)
6717
retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6718
else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6719
retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6720
else if (typ != MatrixType::Rectangular)
6722
(*current_liboctave_error_handler) ("unknown matrix type");
6723
return ComplexMatrix ();
6726
if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6730
retval = qrsolve (*this, b, err);
6732
retval = dmsolve<ComplexMatrix, SparseComplexMatrix, Matrix>
6741
SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b) const
6743
octave_idx_type info;
6745
return solve (mattype, b, info, rcond, 0);
6749
SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
6750
octave_idx_type& info) const
6753
return solve (mattype, b, info, rcond, 0);
6757
SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
6758
octave_idx_type& info, double& rcond) const
6760
return solve (mattype, b, info, rcond, 0);
6764
SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
6765
octave_idx_type& err, double& rcond,
6766
solve_singularity_handler sing_handler,
6767
bool singular_fallback) const
6769
SparseComplexMatrix retval;
6770
int typ = mattype.type (false);
6772
if (typ == MatrixType::Unknown)
6773
typ = mattype.type (*this);
6775
if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
6776
retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6777
else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6778
retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6779
else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6780
retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6781
else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6782
retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6783
else if (typ == MatrixType::Tridiagonal ||
6784
typ == MatrixType::Tridiagonal_Hermitian)
6785
retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6786
else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6787
retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6788
else if (typ != MatrixType::Rectangular)
6790
(*current_liboctave_error_handler) ("unknown matrix type");
6791
return SparseComplexMatrix ();
6794
if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6798
retval = qrsolve (*this, b, err);
6800
retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, SparseMatrix>
6809
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b) const
6811
octave_idx_type info;
6813
return solve (mattype, b, info, rcond, 0);
6817
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
6818
octave_idx_type& info) const
6821
return solve (mattype, b, info, rcond, 0);
6825
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
6826
octave_idx_type& info, double& rcond) const
6828
return solve (mattype, b, info, rcond, 0);
6832
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
6833
octave_idx_type& err, double& rcond,
6834
solve_singularity_handler sing_handler,
6835
bool singular_fallback) const
6837
ComplexMatrix retval;
6838
int typ = mattype.type (false);
6840
if (typ == MatrixType::Unknown)
6841
typ = mattype.type (*this);
6843
if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
6844
retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6845
else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6846
retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6847
else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6848
retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6849
else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6850
retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6851
else if (typ == MatrixType::Tridiagonal ||
6852
typ == MatrixType::Tridiagonal_Hermitian)
6853
retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6854
else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6855
retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6856
else if (typ != MatrixType::Rectangular)
6858
(*current_liboctave_error_handler) ("unknown matrix type");
6859
return ComplexMatrix ();
6862
if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6866
retval = qrsolve (*this, b, err);
6868
retval = dmsolve<ComplexMatrix, SparseComplexMatrix, ComplexMatrix>
6877
SparseComplexMatrix::solve (MatrixType &mattype,
6878
const SparseComplexMatrix& b) const
6880
octave_idx_type info;
6882
return solve (mattype, b, info, rcond, 0);
6886
SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
6887
octave_idx_type& info) const
6890
return solve (mattype, b, info, rcond, 0);
6894
SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
6895
octave_idx_type& info, double& rcond) const
6897
return solve (mattype, b, info, rcond, 0);
6901
SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
6902
octave_idx_type& err, double& rcond,
6903
solve_singularity_handler sing_handler,
6904
bool singular_fallback) const
6906
SparseComplexMatrix retval;
6907
int typ = mattype.type (false);
6909
if (typ == MatrixType::Unknown)
6910
typ = mattype.type (*this);
6912
if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal)
6913
retval = dsolve (mattype, b, err, rcond, sing_handler, false);
6914
else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
6915
retval = utsolve (mattype, b, err, rcond, sing_handler, false);
6916
else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
6917
retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
6918
else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
6919
retval = bsolve (mattype, b, err, rcond, sing_handler, false);
6920
else if (typ == MatrixType::Tridiagonal ||
6921
typ == MatrixType::Tridiagonal_Hermitian)
6922
retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6923
else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6924
retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6925
else if (typ != MatrixType::Rectangular)
6927
(*current_liboctave_error_handler) ("unknown matrix type");
6928
return SparseComplexMatrix ();
6931
if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6935
retval = qrsolve (*this, b, err);
6937
retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
6938
SparseComplexMatrix> (*this, b, err);
6946
SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b) const
6948
octave_idx_type info; double rcond;
6949
return solve (mattype, b, info, rcond);
6953
SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
6954
octave_idx_type& info) const
6957
return solve (mattype, b, info, rcond);
6961
SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
6962
octave_idx_type& info, double& rcond) const
6964
return solve (mattype, b, info, rcond, 0);
6968
SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
6969
octave_idx_type& info, double& rcond,
6970
solve_singularity_handler sing_handler) const
6973
return solve (mattype, tmp, info, rcond,
6974
sing_handler).column (static_cast<octave_idx_type> (0));
6978
SparseComplexMatrix::solve (MatrixType &mattype,
6979
const ComplexColumnVector& b) const
6981
octave_idx_type info;
6983
return solve (mattype, b, info, rcond, 0);
6987
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
6988
octave_idx_type& info) const
6991
return solve (mattype, b, info, rcond, 0);
6995
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
6996
octave_idx_type& info, double& rcond) const
6998
return solve (mattype, b, info, rcond, 0);
7002
SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
7003
octave_idx_type& info, double& rcond,
7004
solve_singularity_handler sing_handler) const
7006
ComplexMatrix tmp (b);
7007
return solve (mattype, tmp, info, rcond,
7008
sing_handler).column (static_cast<octave_idx_type> (0));
7012
SparseComplexMatrix::solve (const Matrix& b) const
7014
octave_idx_type info;
7016
return solve (b, info, rcond, 0);
7020
SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
7023
return solve (b, info, rcond, 0);
7027
SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info,
7028
double& rcond) const
7030
return solve (b, info, rcond, 0);
7034
SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& err,
7036
solve_singularity_handler sing_handler) const
7038
MatrixType mattype (*this);
7039
return solve (mattype, b, err, rcond, sing_handler);
7043
SparseComplexMatrix::solve (const SparseMatrix& b) const
7045
octave_idx_type info;
7047
return solve (b, info, rcond, 0);
7051
SparseComplexMatrix::solve (const SparseMatrix& b,
7052
octave_idx_type& info) const
7055
return solve (b, info, rcond, 0);
7059
SparseComplexMatrix::solve (const SparseMatrix& b,
7060
octave_idx_type& info, double& rcond) const
7062
return solve (b, info, rcond, 0);
7066
SparseComplexMatrix::solve (const SparseMatrix& b,
7067
octave_idx_type& err, double& rcond,
7068
solve_singularity_handler sing_handler) const
7070
MatrixType mattype (*this);
7071
return solve (mattype, b, err, rcond, sing_handler);
7075
SparseComplexMatrix::solve (const ComplexMatrix& b,
7076
octave_idx_type& info) const
7079
return solve (b, info, rcond, 0);
7083
SparseComplexMatrix::solve (const ComplexMatrix& b,
7084
octave_idx_type& info, double& rcond) const
7086
return solve (b, info, rcond, 0);
7090
SparseComplexMatrix::solve (const ComplexMatrix& b,
7091
octave_idx_type& err, double& rcond,
7092
solve_singularity_handler sing_handler) const
7094
MatrixType mattype (*this);
7095
return solve (mattype, b, err, rcond, sing_handler);
7099
SparseComplexMatrix::solve (const SparseComplexMatrix& b) const
7101
octave_idx_type info;
7103
return solve (b, info, rcond, 0);
7107
SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7108
octave_idx_type& info) const
7111
return solve (b, info, rcond, 0);
7115
SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7116
octave_idx_type& info, double& rcond) const
7118
return solve (b, info, rcond, 0);
7122
SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7123
octave_idx_type& err, double& rcond,
7124
solve_singularity_handler sing_handler) const
7126
MatrixType mattype (*this);
7127
return solve (mattype, b, err, rcond, sing_handler);
7131
SparseComplexMatrix::solve (const ColumnVector& b) const
7133
octave_idx_type info; double rcond;
7134
return solve (b, info, rcond);
7138
SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const
7141
return solve (b, info, rcond);
7145
SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
7146
double& rcond) const
7148
return solve (b, info, rcond, 0);
7152
SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
7154
solve_singularity_handler sing_handler) const
7157
return solve (tmp, info, rcond,
7158
sing_handler).column (static_cast<octave_idx_type> (0));
7162
SparseComplexMatrix::solve (const ComplexColumnVector& b) const
7164
octave_idx_type info;
7166
return solve (b, info, rcond, 0);
7170
SparseComplexMatrix::solve (const ComplexColumnVector& b,
7171
octave_idx_type& info) const
7174
return solve (b, info, rcond, 0);
7178
SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
7179
double& rcond) const
7181
return solve (b, info, rcond, 0);
7185
SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
7187
solve_singularity_handler sing_handler) const
7189
ComplexMatrix tmp (b);
7190
return solve (tmp, info, rcond,
7191
sing_handler).column (static_cast<octave_idx_type> (0));
7196
SparseComplexMatrix::operator ! (void) const
7198
if (any_element_is_nan ())
7199
gripe_nan_to_logical_conversion ();
7201
octave_idx_type nr = rows ();
7202
octave_idx_type nc = cols ();
7203
octave_idx_type nz1 = nnz ();
7204
octave_idx_type nz2 = nr*nc - nz1;
7206
SparseBoolMatrix r (nr, nc, nz2);
7208
octave_idx_type ii = 0;
7209
octave_idx_type jj = 0;
7211
for (octave_idx_type i = 0; i < nc; i++)
7213
for (octave_idx_type j = 0; j < nr; j++)
7215
if (jj < cidx (i+1) && ridx (jj) == j)
7230
SparseComplexMatrix::squeeze (void) const
7232
return MSparse<Complex>::squeeze ();
7236
SparseComplexMatrix::reshape (const dim_vector& new_dims) const
7238
return MSparse<Complex>::reshape (new_dims);
7242
SparseComplexMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const
7244
return MSparse<Complex>::permute (vec, inv);
7248
SparseComplexMatrix::ipermute (const Array<octave_idx_type>& vec) const
7250
return MSparse<Complex>::ipermute (vec);
7256
SparseComplexMatrix::any_element_is_nan (void) const
7258
octave_idx_type nel = nnz ();
7260
for (octave_idx_type i = 0; i < nel; i++)
7262
Complex val = data (i);
7271
SparseComplexMatrix::any_element_is_inf_or_nan (void) const
7273
octave_idx_type nel = nnz ();
7275
for (octave_idx_type i = 0; i < nel; i++)
7277
Complex val = data (i);
7278
if (xisinf (val) || xisnan (val))
7285
// Return true if no elements have imaginary components.
7288
SparseComplexMatrix::all_elements_are_real (void) const
7290
return mx_inline_all_real (nnz (), data ());
7293
// Return nonzero if any element of CM has a non-integer real or
7294
// imaginary part. Also extract the largest and smallest (real or
7295
// imaginary) values and return them in MAX_VAL and MIN_VAL.
7298
SparseComplexMatrix::all_integers (double& max_val, double& min_val) const
7300
octave_idx_type nel = nnz ();
7305
max_val = std::real (data (0));
7306
min_val = std::real (data (0));
7308
for (octave_idx_type i = 0; i < nel; i++)
7310
Complex val = data (i);
7312
double r_val = std::real (val);
7313
double i_val = std::imag (val);
7315
if (r_val > max_val)
7318
if (i_val > max_val)
7321
if (r_val < min_val)
7324
if (i_val < min_val)
7327
if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
7335
SparseComplexMatrix::too_large_for_float (void) const
7337
return test_any (xtoo_large_for_float);
7340
// FIXME: Do these really belong here? Maybe they should be in a base class?
7343
SparseComplexMatrix::all (int dim) const
7345
SPARSE_ALL_OP (dim);
7349
SparseComplexMatrix::any (int dim) const
7351
SPARSE_ANY_OP (dim);
7355
SparseComplexMatrix::cumprod (int dim) const
7357
SPARSE_CUMPROD (SparseComplexMatrix, Complex, cumprod);
7361
SparseComplexMatrix::cumsum (int dim) const
7363
SPARSE_CUMSUM (SparseComplexMatrix, Complex, cumsum);
7367
SparseComplexMatrix::prod (int dim) const
7369
if ((rows () == 1 && dim == -1) || dim == 1)
7370
return transpose (). prod (0). transpose ();
7373
SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, *=,
7374
(cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0);
7379
SparseComplexMatrix::sum (int dim) const
7381
SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, +=, 0.0, 0.0);
7385
SparseComplexMatrix::sumsq (int dim) const
7388
Complex d = data (i); \
7389
tmp[ridx (i)] += d * conj (d)
7392
Complex d = data (i); \
7393
tmp[j] += d * conj (d)
7395
SPARSE_BASE_REDUCTION_OP (SparseComplexMatrix, Complex, ROW_EXPR,
7396
COL_EXPR, 0.0, 0.0);
7402
SparseMatrix SparseComplexMatrix::abs (void) const
7404
octave_idx_type nz = nnz ();
7405
octave_idx_type nc = cols ();
7407
SparseMatrix retval (rows (), nc, nz);
7409
for (octave_idx_type i = 0; i < nc + 1; i++)
7410
retval.cidx (i) = cidx (i);
7412
for (octave_idx_type i = 0; i < nz; i++)
7414
retval.data (i) = std::abs (data (i));
7415
retval.ridx (i) = ridx (i);
7422
SparseComplexMatrix::diag (octave_idx_type k) const
7424
return MSparse<Complex>::diag (k);
7428
operator << (std::ostream& os, const SparseComplexMatrix& a)
7430
octave_idx_type nc = a.cols ();
7432
// add one to the printed indices to go from
7433
// zero-based to one-based arrays
7434
for (octave_idx_type j = 0; j < nc; j++)
7437
for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7439
os << a.ridx (i) + 1 << " " << j + 1 << " ";
7440
octave_write_complex (os, a.data (i));
7449
operator >> (std::istream& is, SparseComplexMatrix& a)
7451
typedef SparseComplexMatrix::element_type elt_type;
7453
return read_sparse_matrix<elt_type> (is, a, octave_read_value<Complex>);
7457
operator * (const SparseComplexMatrix& m, const SparseMatrix& a)
7459
SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, double);
7463
operator * (const SparseMatrix& m, const SparseComplexMatrix& a)
7465
SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
7469
operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a)
7471
SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex);
7475
operator * (const ComplexMatrix& m, const SparseMatrix& a)
7477
FULL_SPARSE_MUL (ComplexMatrix, double, Complex (0.,0.));
7481
operator * (const Matrix& m, const SparseComplexMatrix& a)
7483
FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
7487
operator * (const ComplexMatrix& m, const SparseComplexMatrix& a)
7489
FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.));
7493
mul_trans (const ComplexMatrix& m, const SparseComplexMatrix& a)
7495
FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), );
7499
mul_herm (const ComplexMatrix& m, const SparseComplexMatrix& a)
7501
FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), conj);
7505
operator * (const SparseComplexMatrix& m, const Matrix& a)
7507
SPARSE_FULL_MUL (ComplexMatrix, double, Complex (0.,0.));
7511
operator * (const SparseMatrix& m, const ComplexMatrix& a)
7513
SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
7517
operator * (const SparseComplexMatrix& m, const ComplexMatrix& a)
7519
SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.));
7523
trans_mul (const SparseComplexMatrix& m, const ComplexMatrix& a)
7525
SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), );
7529
herm_mul (const SparseComplexMatrix& m, const ComplexMatrix& a)
7531
SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), conj);
7534
// diag * sparse and sparse * diag
7536
operator * (const DiagMatrix& d, const SparseComplexMatrix& a)
7538
return do_mul_dm_sm<SparseComplexMatrix> (d, a);
7541
operator * (const SparseComplexMatrix& a, const DiagMatrix& d)
7543
return do_mul_sm_dm<SparseComplexMatrix> (a, d);
7547
operator * (const ComplexDiagMatrix& d, const SparseMatrix& a)
7549
return do_mul_dm_sm<SparseComplexMatrix> (d, a);
7552
operator * (const SparseMatrix& a, const ComplexDiagMatrix& d)
7554
return do_mul_sm_dm<SparseComplexMatrix> (a, d);
7558
operator * (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
7560
return do_mul_dm_sm<SparseComplexMatrix> (d, a);
7563
operator * (const SparseComplexMatrix& a, const ComplexDiagMatrix& d)
7565
return do_mul_sm_dm<SparseComplexMatrix> (a, d);
7569
operator + (const ComplexDiagMatrix& d, const SparseMatrix& a)
7571
return do_add_dm_sm<SparseComplexMatrix> (d, a);
7574
operator + (const DiagMatrix& d, const SparseComplexMatrix& a)
7576
return do_add_dm_sm<SparseComplexMatrix> (d, a);
7579
operator + (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
7581
return do_add_dm_sm<SparseComplexMatrix> (d, a);
7584
operator + (const SparseMatrix& a, const ComplexDiagMatrix& d)
7586
return do_add_sm_dm<SparseComplexMatrix> (a, d);
7589
operator + (const SparseComplexMatrix& a, const DiagMatrix& d)
7591
return do_add_sm_dm<SparseComplexMatrix> (a, d);
7594
operator + (const SparseComplexMatrix&a, const ComplexDiagMatrix& d)
7596
return do_add_sm_dm<SparseComplexMatrix> (a, d);
7600
operator - (const ComplexDiagMatrix& d, const SparseMatrix& a)
7602
return do_sub_dm_sm<SparseComplexMatrix> (d, a);
7605
operator - (const DiagMatrix& d, const SparseComplexMatrix& a)
7607
return do_sub_dm_sm<SparseComplexMatrix> (d, a);
7610
operator - (const ComplexDiagMatrix& d, const SparseComplexMatrix& a)
7612
return do_sub_dm_sm<SparseComplexMatrix> (d, a);
7615
operator - (const SparseMatrix& a, const ComplexDiagMatrix& d)
7617
return do_sub_sm_dm<SparseComplexMatrix> (a, d);
7620
operator - (const SparseComplexMatrix& a, const DiagMatrix& d)
7622
return do_sub_sm_dm<SparseComplexMatrix> (a, d);
7625
operator - (const SparseComplexMatrix&a, const ComplexDiagMatrix& d)
7627
return do_sub_sm_dm<SparseComplexMatrix> (a, d);
7630
// perm * sparse and sparse * perm
7633
operator * (const PermMatrix& p, const SparseComplexMatrix& a)
7635
return octinternal_do_mul_pm_sm (p, a);
7639
operator * (const SparseComplexMatrix& a, const PermMatrix& p)
7641
return octinternal_do_mul_sm_pm (a, p);
7644
// FIXME: it would be nice to share code among the min/max functions below.
7646
#define EMPTY_RETURN_CHECK(T) \
7647
if (nr == 0 || nc == 0) \
7651
min (const Complex& c, const SparseComplexMatrix& m)
7653
SparseComplexMatrix result;
7655
octave_idx_type nr = m.rows ();
7656
octave_idx_type nc = m.columns ();
7658
EMPTY_RETURN_CHECK (SparseComplexMatrix);
7661
return SparseComplexMatrix (nr, nc);
7664
result = SparseComplexMatrix (m);
7666
for (octave_idx_type j = 0; j < nc; j++)
7667
for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7668
result.data (i) = xmin (c, m.data (i));
7675
min (const SparseComplexMatrix& m, const Complex& c)
7681
min (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
7683
SparseComplexMatrix r;
7685
if ((a.rows () == b.rows ()) && (a.cols () == b.cols ()))
7687
octave_idx_type a_nr = a.rows ();
7688
octave_idx_type a_nc = a.cols ();
7690
octave_idx_type b_nr = b.rows ();
7691
octave_idx_type b_nc = b.cols ();
7693
if (a_nr == 0 || b_nc == 0 || a.nnz () == 0 || b.nnz () == 0)
7694
return SparseComplexMatrix (a_nr, a_nc);
7696
if (a_nr != b_nr || a_nc != b_nc)
7697
gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
7700
r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7702
octave_idx_type jx = 0;
7704
for (octave_idx_type i = 0 ; i < a_nc ; i++)
7706
octave_idx_type ja = a.cidx (i);
7707
octave_idx_type ja_max = a.cidx (i+1);
7708
bool ja_lt_max= ja < ja_max;
7710
octave_idx_type jb = b.cidx (i);
7711
octave_idx_type jb_max = b.cidx (i+1);
7712
bool jb_lt_max = jb < jb_max;
7714
while (ja_lt_max || jb_lt_max )
7717
if ((! jb_lt_max) ||
7718
(ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7720
Complex tmp = xmin (a.data (ja), 0.);
7723
r.ridx (jx) = a.ridx (ja);
7728
ja_lt_max= ja < ja_max;
7730
else if (( !ja_lt_max ) ||
7731
(jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
7733
Complex tmp = xmin (0., b.data (jb));
7736
r.ridx (jx) = b.ridx (jb);
7741
jb_lt_max= jb < jb_max;
7745
Complex tmp = xmin (a.data (ja), b.data (jb));
7749
r.ridx (jx) = a.ridx (ja);
7753
ja_lt_max= ja < ja_max;
7755
jb_lt_max= jb < jb_max;
7761
r.maybe_compress ();
7765
(*current_liboctave_error_handler) ("matrix size mismatch");
7771
max (const Complex& c, const SparseComplexMatrix& m)
7773
SparseComplexMatrix result;
7775
octave_idx_type nr = m.rows ();
7776
octave_idx_type nc = m.columns ();
7778
EMPTY_RETURN_CHECK (SparseComplexMatrix);
7780
// Count the number of non-zero elements
7781
if (xmax (c, 0.) != 0.)
7783
result = SparseComplexMatrix (nr, nc, c);
7784
for (octave_idx_type j = 0; j < nc; j++)
7785
for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
7786
result.xdata (m.ridx (i) + j * nr) = xmax (c, m.data (i));
7789
result = SparseComplexMatrix (m);
7795
max (const SparseComplexMatrix& m, const Complex& c)
7801
max (const SparseComplexMatrix& a, const SparseComplexMatrix& b)
7803
SparseComplexMatrix r;
7805
if ((a.rows () == b.rows ()) && (a.cols () == b.cols ()))
7807
octave_idx_type a_nr = a.rows ();
7808
octave_idx_type a_nc = a.cols ();
7810
octave_idx_type b_nr = b.rows ();
7811
octave_idx_type b_nc = b.cols ();
7813
if (a_nr == 0 || b_nc == 0)
7814
return SparseComplexMatrix (a_nr, a_nc);
7816
return SparseComplexMatrix (b);
7818
return SparseComplexMatrix (a);
7820
if (a_nr != b_nr || a_nc != b_nc)
7821
gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
7824
r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
7826
octave_idx_type jx = 0;
7828
for (octave_idx_type i = 0 ; i < a_nc ; i++)
7830
octave_idx_type ja = a.cidx (i);
7831
octave_idx_type ja_max = a.cidx (i+1);
7832
bool ja_lt_max= ja < ja_max;
7834
octave_idx_type jb = b.cidx (i);
7835
octave_idx_type jb_max = b.cidx (i+1);
7836
bool jb_lt_max = jb < jb_max;
7838
while (ja_lt_max || jb_lt_max )
7841
if ((! jb_lt_max) ||
7842
(ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
7844
Complex tmp = xmax (a.data (ja), 0.);
7847
r.ridx (jx) = a.ridx (ja);
7852
ja_lt_max= ja < ja_max;
7854
else if (( !ja_lt_max ) ||
7855
(jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
7857
Complex tmp = xmax (0., b.data (jb));
7860
r.ridx (jx) = b.ridx (jb);
7865
jb_lt_max= jb < jb_max;
7869
Complex tmp = xmax (a.data (ja), b.data (jb));
7873
r.ridx (jx) = a.ridx (ja);
7877
ja_lt_max= ja < ja_max;
7879
jb_lt_max= jb < jb_max;
7885
r.maybe_compress ();
7889
(*current_liboctave_error_handler) ("matrix size mismatch");
7894
SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex,
7896
SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0)
7898
SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix,
7900
SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0)
7902
SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix,
7904
SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0)