GNU Octave  8.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-2023 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 
47 
48 static void
50 {
52 }
53 
54 template <typename T1, typename T2>
55 bool
56 mx_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 
88 template <typename T1, typename T2>
89 bool
90 mx_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 -*-
136 Matrix
137 xdiv (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 -*-
157 xdiv (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 -*-
177 xdiv (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 -*-
217 xdiv (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;
248  SparseComplexMatrix result
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;
268  SparseComplexMatrix result
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;
289  SparseComplexMatrix result
290  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
291 
292  typ = btyp.transpose ();
293  return result.hermitian ();
294 }
295 
296 template <typename RT, typename SM, typename DM>
297 RT 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 -*-
339 xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType&)
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 
374 Matrix
375 elem_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 
399 elem_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 
418 elem_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 -*-
467 Matrix
468 xleftdiv (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 -*-
480 xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, MatrixType& typ)
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 -*-
492 xleftdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
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 -*-
516 xleftdiv (const SparseComplexMatrix& a, const Matrix& b, MatrixType& typ)
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 
563 template <typename RT, typename DM, typename SM>
564 RT 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 
OCTAVE_END_NAMESPACE(octave)
#define Inf
Definition: Faddeeva.cc:260
#define NaN
Definition: Faddeeva.cc:261
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
OCTARRAY_OVERRIDABLE_FUNC_API T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
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
octave_idx_type * cidx(void)
Definition: Sparse.h:596
octave_idx_type * ridx(void)
Definition: Sparse.h:583
T * data(void)
Definition: Sparse.h:574
octave_idx_type cols(void) const
Definition: Sparse.h:352
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void warn_singular_matrix(double rcond)
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
T * r
Definition: mx-inlines.cc:773
std::complex< double > Complex
Definition: oct-cmplx.h:33
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
static void solve_singularity_warning(double rcond)
Definition: sparse-xdiv.cc:49
#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