GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sparse-xdiv.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26#if defined (HAVE_CONFIG_H)
27# include "config.h"
28#endif
29
30#include <cassert>
31
32#include "Array-util.h"
33#include "lo-array-errwarn.h"
34#include "oct-cmplx.h"
35#include "quit.h"
36#include "error.h"
37#include "lo-ieee.h"
38
39#include "dSparse.h"
40#include "dDiagMatrix.h"
41#include "CSparse.h"
42#include "CDiagMatrix.h"
43#include "oct-spparms.h"
44#include "sparse-xdiv.h"
45
46OCTAVE_NAMESPACE_BEGIN
47
48static void
50{
52}
53
54template <typename T1, typename T2>
55bool
56mx_leftdiv_conform (const T1& a, const T2& b)
57{
58 octave_idx_type a_nr = a.rows ();
59 octave_idx_type b_nr = b.rows ();
60
61 if (a_nr != b_nr)
62 {
63 octave_idx_type a_nc = a.cols ();
64 octave_idx_type b_nc = b.cols ();
65
66 octave::err_nonconformant (R"(operator \)", a_nr, a_nc, b_nr, b_nc);
67 }
68
69 return true;
70}
71
72#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
73 template bool mx_leftdiv_conform (const T1&, const T2&)
74
87
88template <typename T1, typename T2>
89bool
90mx_div_conform (const T1& a, const T2& b)
91{
92 octave_idx_type a_nc = a.cols ();
93 octave_idx_type b_nc = b.cols ();
94
95 if (a_nc != b_nc)
96 {
97 octave_idx_type a_nr = a.rows ();
98 octave_idx_type b_nr = b.rows ();
99
100 octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
101 }
102
103 return true;
104}
105
106#define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \
107 template bool mx_div_conform (const T1&, const T2&)
108
121
122// Right division functions. X / Y = X * inv (Y) = (inv (Y') * X')'
123//
124// Y / X: m cm sm scm
125// +-- +---+----+----+----+
126// sparse matrix | 1 | 3 | 5 | 7 |
127// +---+----+----+----+
128// sparse complex_matrix | 2 | 4 | 6 | 8 |
129// +---+----+----+----+
130// diagonal matrix | 9 | 11 |
131// +----+----+
132// complex diag. matrix | 10 | 12 |
133// +----+----+
134
135// -*- 1 -*-
136Matrix
137xdiv (const Matrix& a, const SparseMatrix& b, MatrixType& typ)
138{
139 if (! mx_div_conform (a, b))
140 return Matrix ();
141
142 Matrix atmp = a.transpose ();
143 SparseMatrix btmp = b.transpose ();
144 MatrixType btyp = typ.transpose ();
145
146 octave_idx_type info;
147 double rcond = 0.0;
148 Matrix result = btmp.solve (btyp, atmp, info, rcond,
150
151 typ = btyp.transpose ();
152 return result.transpose ();
153}
154
155// -*- 2 -*-
157xdiv (const Matrix& a, const SparseComplexMatrix& b, MatrixType& typ)
158{
159 if (! mx_div_conform (a, b))
160 return ComplexMatrix ();
161
162 Matrix atmp = a.transpose ();
163 SparseComplexMatrix btmp = b.hermitian ();
164 MatrixType btyp = typ.transpose ();
165
166 octave_idx_type info;
167 double rcond = 0.0;
168 ComplexMatrix result
169 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
170
171 typ = btyp.transpose ();
172 return result.hermitian ();
173}
174
175// -*- 3 -*-
177xdiv (const ComplexMatrix& a, const SparseMatrix& b, MatrixType& typ)
178{
179 if (! mx_div_conform (a, b))
180 return ComplexMatrix ();
181
182 ComplexMatrix atmp = a.hermitian ();
183 SparseMatrix btmp = b.transpose ();
184 MatrixType btyp = typ.transpose ();
185
186 octave_idx_type info;
187 double rcond = 0.0;
188 ComplexMatrix result
189 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
190
191 typ = btyp.transpose ();
192 return result.hermitian ();
193}
194
195// -*- 4 -*-
198{
199 if (! mx_div_conform (a, b))
200 return ComplexMatrix ();
201
202 ComplexMatrix atmp = a.hermitian ();
203 SparseComplexMatrix btmp = b.hermitian ();
204 MatrixType btyp = typ.transpose ();
205
206 octave_idx_type info;
207 double rcond = 0.0;
208 ComplexMatrix result
209 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
210
211 typ = btyp.transpose ();
212 return result.hermitian ();
213}
214
215// -*- 5 -*-
217xdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
218{
219 if (! mx_div_conform (a, b))
220 return SparseMatrix ();
221
222 SparseMatrix atmp = a.transpose ();
223 SparseMatrix btmp = b.transpose ();
224 MatrixType btyp = typ.transpose ();
225
226 octave_idx_type info;
227 double rcond = 0.0;
228 SparseMatrix result = btmp.solve (btyp, atmp, info, rcond,
230
231 typ = btyp.transpose ();
232 return result.transpose ();
233}
234
235// -*- 6 -*-
238{
239 if (! mx_div_conform (a, b))
240 return SparseComplexMatrix ();
241
242 SparseMatrix atmp = a.transpose ();
243 SparseComplexMatrix btmp = b.hermitian ();
244 MatrixType btyp = typ.transpose ();
245
246 octave_idx_type info;
247 double rcond = 0.0;
249 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
250
251 typ = btyp.transpose ();
252 return result.hermitian ();
253}
254
255// -*- 7 -*-
258{
259 if (! mx_div_conform (a, b))
260 return SparseComplexMatrix ();
261
262 SparseComplexMatrix atmp = a.hermitian ();
263 SparseMatrix btmp = b.transpose ();
264 MatrixType btyp = typ.transpose ();
265
266 octave_idx_type info;
267 double rcond = 0.0;
269 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
270
271 typ = btyp.transpose ();
272 return result.hermitian ();
273}
274
275// -*- 8 -*-
278 MatrixType& typ)
279{
280 if (! mx_div_conform (a, b))
281 return SparseComplexMatrix ();
282
283 SparseComplexMatrix atmp = a.hermitian ();
284 SparseComplexMatrix btmp = b.hermitian ();
285 MatrixType btyp = typ.transpose ();
286
287 octave_idx_type info;
288 double rcond = 0.0;
290 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
291
292 typ = btyp.transpose ();
293 return result.hermitian ();
294}
295
296template <typename RT, typename SM, typename DM>
297RT do_rightdiv_sm_dm (const SM& a, const DM& d)
298{
299 const octave_idx_type d_nr = d.rows ();
300
301 const octave_idx_type a_nr = a.rows ();
302 const octave_idx_type a_nc = a.cols ();
303
304 using std::min;
305 const octave_idx_type nc = min (d_nr, a_nc);
306
307 if (! mx_div_conform (a, d))
308 return RT ();
309
310 const octave_idx_type nz = a.nnz ();
311 RT r (a_nr, nc, nz);
312
313 typedef typename DM::element_type DM_elt_type;
314 const DM_elt_type zero = DM_elt_type ();
315
316 octave_idx_type k_result = 0;
317 for (octave_idx_type j = 0; j < nc; ++j)
318 {
319 octave_quit ();
320 const DM_elt_type s = d.dgelem (j);
321 const octave_idx_type colend = a.cidx (j+1);
322 r.xcidx (j) = k_result;
323 if (s != zero)
324 for (octave_idx_type k = a.cidx (j); k < colend; ++k)
325 {
326 r.xdata (k_result) = a.data (k) / s;
327 r.xridx (k_result) = a.ridx (k);
328 ++k_result;
329 }
330 }
331 r.xcidx (nc) = k_result;
332
333 r.maybe_compress (true);
334 return r;
335}
336
337// -*- 9 -*-
340{
341 return do_rightdiv_sm_dm<SparseMatrix> (a, b);
342}
343
344// -*- 10 -*-
347{
348 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
349}
350
351// -*- 11 -*-
354{
355 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
356}
357
358// -*- 12 -*-
361{
362 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
363}
364
365// Funny element by element division operations.
366//
367// op2 \ op1: s cs
368// +-- +---+----+
369// matrix | 1 | 3 |
370// +---+----+
371// complex_matrix | 2 | 4 |
372// +---+----+
373
374Matrix
375elem_xdiv (double a, const SparseMatrix& b)
376{
377 octave_idx_type nr = b.rows ();
378 octave_idx_type nc = b.cols ();
379
380 Matrix result;
381 if (a == 0.)
382 result = Matrix (nr, nc, octave::numeric_limits<double>::NaN ());
383 else if (a > 0.)
384 result = Matrix (nr, nc, octave::numeric_limits<double>::Inf ());
385 else
386 result = Matrix (nr, nc, -octave::numeric_limits<double>::Inf ());
387
388 for (octave_idx_type j = 0; j < nc; j++)
389 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
390 {
391 octave_quit ();
392 result.elem (b.ridx (i), j) = a / b.data (i);
393 }
394
395 return result;
396}
397
399elem_xdiv (double a, const SparseComplexMatrix& b)
400{
401 octave_idx_type nr = b.rows ();
402 octave_idx_type nc = b.cols ();
403
406
407 for (octave_idx_type j = 0; j < nc; j++)
408 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
409 {
410 octave_quit ();
411 result.elem (b.ridx (i), j) = a / b.data (i);
412 }
413
414 return result;
415}
416
418elem_xdiv (const Complex& a, const SparseMatrix& b)
419{
420 octave_idx_type nr = b.rows ();
421 octave_idx_type nc = b.cols ();
422
423 ComplexMatrix result (nr, nc, (a / 0.0));
424
425 for (octave_idx_type j = 0; j < nc; j++)
426 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
427 {
428 octave_quit ();
429 result.elem (b.ridx (i), j) = a / b.data (i);
430 }
431
432 return result;
433}
434
437{
438 octave_idx_type nr = b.rows ();
439 octave_idx_type nc = b.cols ();
440
441 ComplexMatrix result (nr, nc, (a / 0.0));
442
443 for (octave_idx_type j = 0; j < nc; j++)
444 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
445 {
446 octave_quit ();
447 result.elem (b.ridx (i), j) = a / b.data (i);
448 }
449
450 return result;
451}
452
453// Left division functions. X \ Y = inv (X) * Y
454//
455// Y \ X : sm scm dm dcm
456// +-- +---+----+
457// matrix | 1 | 5 |
458// +---+----+
459// complex_matrix | 2 | 6 |
460// +---+----+----+----+
461// sparse matrix | 3 | 7 | 9 | 11 |
462// +---+----+----+----+
463// sparse complex_matrix | 4 | 8 | 10 | 12 |
464// +---+----+----+----+
465
466// -*- 1 -*-
467Matrix
468xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ)
469{
470 if (! mx_leftdiv_conform (a, b))
471 return Matrix ();
472
473 octave_idx_type info;
474 double rcond = 0.0;
475 return a.solve (typ, b, info, rcond, solve_singularity_warning);
476}
477
478// -*- 2 -*-
481{
482 if (! mx_leftdiv_conform (a, b))
483 return ComplexMatrix ();
484
485 octave_idx_type info;
486 double rcond = 0.0;
487 return a.solve (typ, b, info, rcond, solve_singularity_warning);
488}
489
490// -*- 3 -*-
493{
494 if (! mx_leftdiv_conform (a, b))
495 return SparseMatrix ();
496
497 octave_idx_type info;
498 double rcond = 0.0;
499 return a.solve (typ, b, info, rcond, solve_singularity_warning);
500}
501
502// -*- 4 -*-
505{
506 if (! mx_leftdiv_conform (a, b))
507 return SparseComplexMatrix ();
508
509 octave_idx_type info;
510 double rcond = 0.0;
511 return a.solve (typ, b, info, rcond, solve_singularity_warning);
512}
513
514// -*- 5 -*-
517{
518 if (! mx_leftdiv_conform (a, b))
519 return ComplexMatrix ();
520
521 octave_idx_type info;
522 double rcond = 0.0;
523 return a.solve (typ, b, info, rcond, solve_singularity_warning);
524}
525
526// -*- 6 -*-
529{
530 if (! mx_leftdiv_conform (a, b))
531 return ComplexMatrix ();
532
533 octave_idx_type info;
534 double rcond = 0.0;
535 return a.solve (typ, b, info, rcond, solve_singularity_warning);
536}
537
538// -*- 7 -*-
541{
542 if (! mx_leftdiv_conform (a, b))
543 return SparseComplexMatrix ();
544
545 octave_idx_type info;
546 double rcond = 0.0;
547 return a.solve (typ, b, info, rcond, solve_singularity_warning);
548}
549
550// -*- 8 -*-
553 MatrixType& typ)
554{
555 if (! mx_leftdiv_conform (a, b))
556 return SparseComplexMatrix ();
557
558 octave_idx_type info;
559 double rcond = 0.0;
560 return a.solve (typ, b, info, rcond, solve_singularity_warning);
561}
562
563template <typename RT, typename DM, typename SM>
564RT do_leftdiv_dm_sm (const DM& d, const SM& a)
565{
566 const octave_idx_type a_nr = a.rows ();
567 const octave_idx_type a_nc = a.cols ();
568
569 const octave_idx_type d_nc = d.cols ();
570
571 using std::min;
572 const octave_idx_type nr = min (d_nc, a_nr);
573
574 if (! mx_leftdiv_conform (d, a))
575 return RT ();
576
577 const octave_idx_type nz = a.nnz ();
578 RT r (nr, a_nc, nz);
579
580 typedef typename DM::element_type DM_elt_type;
581 const DM_elt_type zero = DM_elt_type ();
582
583 octave_idx_type k_result = 0;
584 for (octave_idx_type j = 0; j < a_nc; ++j)
585 {
586 octave_quit ();
587 const octave_idx_type colend = a.cidx (j+1);
588 r.xcidx (j) = k_result;
589 for (octave_idx_type k = a.cidx (j); k < colend; ++k)
590 {
591 const octave_idx_type i = a.ridx (k);
592 if (i < nr)
593 {
594 const DM_elt_type s = d.dgelem (i);
595 if (s != zero)
596 {
597 r.xdata (k_result) = a.data (k) / s;
598 r.xridx (k_result) = i;
599 ++k_result;
600 }
601 }
602 }
603 }
604 r.xcidx (a_nc) = k_result;
605
606 r.maybe_compress (true);
607 return r;
608}
609
610// -*- 9 -*-
613{
614 return do_leftdiv_dm_sm<SparseMatrix> (d, a);
615}
616
617// -*- 10 -*-
620{
621 return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
622}
623
624// -*- 11 -*-
627{
628 return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
629}
630
631// -*- 12 -*-
634 MatrixType&)
635{
636 return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
637}
638
639OCTAVE_NAMESPACE_END
#define Inf
Definition: Faddeeva.cc:260
#define NaN
Definition: Faddeeva.cc:261
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:534
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:170
OCTAVE_API MatrixType transpose(void) const
Definition: MatrixType.cc:973
Definition: dMatrix.h:42
Matrix transpose(void) const
Definition: dMatrix.h:140
OCTAVE_API SparseComplexMatrix hermitian(void) const
Definition: CSparse.cc:606
OCTAVE_API ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition: CSparse.cc:6702
SparseMatrix transpose(void) const
Definition: dSparse.h:124
OCTAVE_API Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6692
octave_idx_type rows(void) const
Definition: Sparse.h:351
T * data(void)
Definition: Sparse.h:574
octave_idx_type * ridx(void)
Definition: Sparse.h:583
octave_idx_type * cidx(void)
Definition: Sparse.h:596
octave_idx_type cols(void) const
Definition: Sparse.h:352
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
class OCTAVE_API SparseMatrix
Definition: mx-fwd.h:55
class OCTAVE_API SparseComplexMatrix
Definition: mx-fwd.h:56
void warn_singular_matrix(double rcond)
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
std::complex< double > Complex
Definition: oct-cmplx.h:33
static OCTAVE_NAMESPACE_BEGIN void solve_singularity_warning(double rcond)
Definition: sparse-xdiv.cc:49
Matrix elem_xdiv(double a, const SparseMatrix &b)
Definition: sparse-xdiv.cc:375
Matrix xdiv(const Matrix &a, const SparseMatrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:137
bool mx_leftdiv_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:56
bool mx_div_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:90
RT do_rightdiv_sm_dm(const SM &a, const DM &d)
Definition: sparse-xdiv.cc:297
#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:72
#define INSTANTIATE_MX_DIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:106
Matrix xleftdiv(const SparseMatrix &a, const Matrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:468
RT do_leftdiv_dm_sm(const DM &d, const SM &a)
Definition: sparse-xdiv.cc:564