GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
sparse-xdiv.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1998-2025 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 "Array-util.h"
31#include "lo-array-errwarn.h"
32#include "oct-cmplx.h"
33#include "quit.h"
34#include "error.h"
35#include "lo-ieee.h"
36
37#include "dSparse.h"
38#include "dDiagMatrix.h"
39#include "CSparse.h"
40#include "CDiagMatrix.h"
41#include "oct-spparms.h"
42#include "sparse-xdiv.h"
43
45
46static void
47solve_singularity_warning (double rcond)
48{
49 octave::warn_singular_matrix (rcond);
50}
51
52template <typename T1, typename T2>
53bool
54mx_leftdiv_conform (const T1& a, const T2& b)
55{
56 octave_idx_type a_nr = a.rows ();
57 octave_idx_type b_nr = b.rows ();
58
59 if (a_nr != b_nr)
60 {
61 octave_idx_type a_nc = a.cols ();
62 octave_idx_type b_nc = b.cols ();
63
64 octave::err_nonconformant (R"(operator \)", a_nr, a_nc, b_nr, b_nc);
65 }
66
67 return true;
68}
69
70#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
71 template bool mx_leftdiv_conform (const T1&, const T2&)
72
85
86template <typename T1, typename T2>
87bool
88mx_div_conform (const T1& a, const T2& b)
89{
90 octave_idx_type a_nc = a.cols ();
91 octave_idx_type b_nc = b.cols ();
92
93 if (a_nc != b_nc)
94 {
95 octave_idx_type a_nr = a.rows ();
96 octave_idx_type b_nr = b.rows ();
97
98 octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
99 }
100
101 return true;
102}
103
104#define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \
105 template bool mx_div_conform (const T1&, const T2&)
106
119
120// Right division functions. X / Y = X * inv (Y) = (inv (Y') * X')'
121//
122// Y / X: m cm sm scm
123// +-- +---+----+----+----+
124// sparse matrix | 1 | 3 | 5 | 7 |
125// +---+----+----+----+
126// sparse complex_matrix | 2 | 4 | 6 | 8 |
127// +---+----+----+----+
128// diagonal matrix | 9 | 11 |
129// +----+----+
130// complex diag. matrix | 10 | 12 |
131// +----+----+
132
133// -*- 1 -*-
134Matrix
135xdiv (const Matrix& a, const SparseMatrix& b, MatrixType& typ)
136{
137 if (! mx_div_conform (a, b))
138 return Matrix ();
139
140 Matrix atmp = a.transpose ();
141 SparseMatrix btmp = b.transpose ();
142 MatrixType btyp = typ.transpose ();
143
144 octave_idx_type info;
145 double rcond = 0.0;
146 Matrix result = btmp.solve (btyp, atmp, info, rcond,
147 solve_singularity_warning);
148
149 typ = btyp.transpose ();
150 return result.transpose ();
151}
152
153// -*- 2 -*-
155xdiv (const Matrix& a, const SparseComplexMatrix& b, MatrixType& typ)
156{
157 if (! mx_div_conform (a, b))
158 return ComplexMatrix ();
159
160 Matrix atmp = a.transpose ();
161 SparseComplexMatrix btmp = b.hermitian ();
162 MatrixType btyp = typ.transpose ();
163
164 octave_idx_type info;
165 double rcond = 0.0;
166 ComplexMatrix result
167 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
168
169 typ = btyp.transpose ();
170 return result.hermitian ();
171}
172
173// -*- 3 -*-
175xdiv (const ComplexMatrix& a, const SparseMatrix& b, MatrixType& typ)
176{
177 if (! mx_div_conform (a, b))
178 return ComplexMatrix ();
179
180 ComplexMatrix atmp = a.hermitian ();
181 SparseMatrix btmp = b.transpose ();
182 MatrixType btyp = typ.transpose ();
183
184 octave_idx_type info;
185 double rcond = 0.0;
186 ComplexMatrix result
187 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
188
189 typ = btyp.transpose ();
190 return result.hermitian ();
191}
192
193// -*- 4 -*-
196{
197 if (! mx_div_conform (a, b))
198 return ComplexMatrix ();
199
200 ComplexMatrix atmp = a.hermitian ();
201 SparseComplexMatrix btmp = b.hermitian ();
202 MatrixType btyp = typ.transpose ();
203
204 octave_idx_type info;
205 double rcond = 0.0;
206 ComplexMatrix result
207 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
208
209 typ = btyp.transpose ();
210 return result.hermitian ();
211}
212
213// -*- 5 -*-
215xdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
216{
217 if (! mx_div_conform (a, b))
218 return SparseMatrix ();
219
220 SparseMatrix atmp = a.transpose ();
221 SparseMatrix btmp = b.transpose ();
222 MatrixType btyp = typ.transpose ();
223
224 octave_idx_type info;
225 double rcond = 0.0;
226 SparseMatrix result = btmp.solve (btyp, atmp, info, rcond,
227 solve_singularity_warning);
228
229 typ = btyp.transpose ();
230 return result.transpose ();
231}
232
233// -*- 6 -*-
236{
237 if (! mx_div_conform (a, b))
238 return SparseComplexMatrix ();
239
240 SparseMatrix atmp = a.transpose ();
241 SparseComplexMatrix btmp = b.hermitian ();
242 MatrixType btyp = typ.transpose ();
243
244 octave_idx_type info;
245 double rcond = 0.0;
247 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
248
249 typ = btyp.transpose ();
250 return result.hermitian ();
251}
252
253// -*- 7 -*-
256{
257 if (! mx_div_conform (a, b))
258 return SparseComplexMatrix ();
259
260 SparseComplexMatrix atmp = a.hermitian ();
261 SparseMatrix btmp = b.transpose ();
262 MatrixType btyp = typ.transpose ();
263
264 octave_idx_type info;
265 double rcond = 0.0;
267 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
268
269 typ = btyp.transpose ();
270 return result.hermitian ();
271}
272
273// -*- 8 -*-
276 MatrixType& typ)
277{
278 if (! mx_div_conform (a, b))
279 return SparseComplexMatrix ();
280
281 SparseComplexMatrix atmp = a.hermitian ();
282 SparseComplexMatrix btmp = b.hermitian ();
283 MatrixType btyp = typ.transpose ();
284
285 octave_idx_type info;
286 double rcond = 0.0;
288 = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
289
290 typ = btyp.transpose ();
291 return result.hermitian ();
292}
293
294template <typename RT, typename SM, typename DM>
295RT
296do_rightdiv_sm_dm (const SM& a, const DM& d)
297{
298 const octave_idx_type d_nr = d.rows ();
299
300 const octave_idx_type a_nr = a.rows ();
301 const octave_idx_type a_nc = a.cols ();
302
303 using std::min;
304 const octave_idx_type nc = min (d_nr, a_nc);
305
306 if (! mx_div_conform (a, d))
307 return RT ();
308
309 const octave_idx_type nz = a.nnz ();
310 RT r (a_nr, nc, nz);
311
312 typedef typename DM::element_type DM_elt_type;
313 const DM_elt_type zero = DM_elt_type ();
314
315 octave_idx_type k_result = 0;
316 for (octave_idx_type j = 0; j < nc; ++j)
317 {
318 octave_quit ();
319 const DM_elt_type s = d.dgelem (j);
320 const octave_idx_type colend = a.cidx (j+1);
321 r.xcidx (j) = k_result;
322 if (s != zero)
323 for (octave_idx_type k = a.cidx (j); k < colend; ++k)
324 {
325 r.xdata (k_result) = a.data (k) / s;
326 r.xridx (k_result) = a.ridx (k);
327 ++k_result;
328 }
329 }
330 r.xcidx (nc) = k_result;
331
332 r.maybe_compress (true);
333 return r;
334}
335
336// -*- 9 -*-
339{
340 return do_rightdiv_sm_dm<SparseMatrix> (a, b);
341}
342
343// -*- 10 -*-
346{
347 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
348}
349
350// -*- 11 -*-
353{
354 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
355}
356
357// -*- 12 -*-
360{
361 return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
362}
363
364// Funny element by element division operations.
365//
366// op2 \ op1: s cs
367// +-- +---+----+
368// matrix | 1 | 3 |
369// +---+----+
370// complex_matrix | 2 | 4 |
371// +---+----+
372
373Matrix
374elem_xdiv (double a, const SparseMatrix& b)
375{
376 octave_idx_type nr = b.rows ();
377 octave_idx_type nc = b.cols ();
378
379 Matrix result;
380 if (a == 0.)
381 result = Matrix (nr, nc, octave::numeric_limits<double>::NaN ());
382 else if (a > 0.)
383 result = Matrix (nr, nc, octave::numeric_limits<double>::Inf ());
384 else
385 result = Matrix (nr, nc, -octave::numeric_limits<double>::Inf ());
386
387 for (octave_idx_type j = 0; j < nc; j++)
388 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
389 {
390 octave_quit ();
391 result.elem (b.ridx (i), j) = a / b.data (i);
392 }
393
394 return result;
395}
396
398elem_xdiv (double a, const SparseComplexMatrix& b)
399{
400 octave_idx_type nr = b.rows ();
401 octave_idx_type nc = b.cols ();
402
403 ComplexMatrix result (nr, nc, Complex (octave::numeric_limits<double>::NaN (),
404 octave::numeric_limits<double>::NaN ()));
405
406 for (octave_idx_type j = 0; j < nc; j++)
407 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
408 {
409 octave_quit ();
410 result.elem (b.ridx (i), j) = a / b.data (i);
411 }
412
413 return result;
414}
415
417elem_xdiv (const Complex& a, const SparseMatrix& b)
418{
419 octave_idx_type nr = b.rows ();
420 octave_idx_type nc = b.cols ();
421
422 ComplexMatrix result (nr, nc, (a / 0.0));
423
424 for (octave_idx_type j = 0; j < nc; j++)
425 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
426 {
427 octave_quit ();
428 result.elem (b.ridx (i), j) = a / b.data (i);
429 }
430
431 return result;
432}
433
436{
437 octave_idx_type nr = b.rows ();
438 octave_idx_type nc = b.cols ();
439
440 ComplexMatrix result (nr, nc, (a / 0.0));
441
442 for (octave_idx_type j = 0; j < nc; j++)
443 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
444 {
445 octave_quit ();
446 result.elem (b.ridx (i), j) = a / b.data (i);
447 }
448
449 return result;
450}
451
452// Left division functions. X \ Y = inv (X) * Y
453//
454// Y \ X : sm scm dm dcm
455// +-- +---+----+
456// matrix | 1 | 5 |
457// +---+----+
458// complex_matrix | 2 | 6 |
459// +---+----+----+----+
460// sparse matrix | 3 | 7 | 9 | 11 |
461// +---+----+----+----+
462// sparse complex_matrix | 4 | 8 | 10 | 12 |
463// +---+----+----+----+
464
465// -*- 1 -*-
466Matrix
467xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ)
468{
469 if (! mx_leftdiv_conform (a, b))
470 return Matrix ();
471
472 octave_idx_type info;
473 double rcond = 0.0;
474 return a.solve (typ, b, info, rcond, solve_singularity_warning);
475}
476
477// -*- 2 -*-
480{
481 if (! mx_leftdiv_conform (a, b))
482 return ComplexMatrix ();
483
484 octave_idx_type info;
485 double rcond = 0.0;
486 return a.solve (typ, b, info, rcond, solve_singularity_warning);
487}
488
489// -*- 3 -*-
492{
493 if (! mx_leftdiv_conform (a, b))
494 return SparseMatrix ();
495
496 octave_idx_type info;
497 double rcond = 0.0;
498 return a.solve (typ, b, info, rcond, solve_singularity_warning);
499}
500
501// -*- 4 -*-
504{
505 if (! mx_leftdiv_conform (a, b))
506 return SparseComplexMatrix ();
507
508 octave_idx_type info;
509 double rcond = 0.0;
510 return a.solve (typ, b, info, rcond, solve_singularity_warning);
511}
512
513// -*- 5 -*-
516{
517 if (! mx_leftdiv_conform (a, b))
518 return ComplexMatrix ();
519
520 octave_idx_type info;
521 double rcond = 0.0;
522 return a.solve (typ, b, info, rcond, solve_singularity_warning);
523}
524
525// -*- 6 -*-
528{
529 if (! mx_leftdiv_conform (a, b))
530 return ComplexMatrix ();
531
532 octave_idx_type info;
533 double rcond = 0.0;
534 return a.solve (typ, b, info, rcond, solve_singularity_warning);
535}
536
537// -*- 7 -*-
540{
541 if (! mx_leftdiv_conform (a, b))
542 return SparseComplexMatrix ();
543
544 octave_idx_type info;
545 double rcond = 0.0;
546 return a.solve (typ, b, info, rcond, solve_singularity_warning);
547}
548
549// -*- 8 -*-
552 MatrixType& typ)
553{
554 if (! mx_leftdiv_conform (a, b))
555 return SparseComplexMatrix ();
556
557 octave_idx_type info;
558 double rcond = 0.0;
559 return a.solve (typ, b, info, rcond, solve_singularity_warning);
560}
561
562template <typename RT, typename DM, typename SM>
563RT
564do_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_END_NAMESPACE(octave)
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:563
ComplexMatrix hermitian() const
Definition CMatrix.h:168
MatrixType transpose() const
Matrix transpose() const
Definition dMatrix.h:138
SparseComplexMatrix hermitian() const
Definition CSparse.cc:604
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition CSparse.cc:6669
SparseMatrix transpose() const
Definition dSparse.h:125
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition dSparse.cc:6661
octave_idx_type cols() const
Definition Sparse.h:349
octave_idx_type * cidx()
Definition Sparse.h:593
T * data()
Definition Sparse.h:571
octave_idx_type * ridx()
Definition Sparse.h:580
octave_idx_type rows() const
Definition Sparse.h:348
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
std::complex< double > Complex
Definition oct-cmplx.h:33
Matrix elem_xdiv(double a, const SparseMatrix &b)
Matrix xdiv(const Matrix &a, const SparseMatrix &b, MatrixType &typ)
bool mx_leftdiv_conform(const T1 &a, const T2 &b)
bool mx_div_conform(const T1 &a, const T2 &b)
RT do_rightdiv_sm_dm(const SM &a, const DM &d)
#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2)
#define INSTANTIATE_MX_DIV_CONFORM(T1, T2)
Matrix xleftdiv(const SparseMatrix &a, const Matrix &b, MatrixType &typ)
RT do_leftdiv_dm_sm(const DM &d, const SM &a)