GNU Octave  6.2.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-2021 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 
46 static void
48 {
50 }
51 
52 template <typename T1, typename T2>
53 bool
54 mx_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 
86 template <typename T1, typename T2>
87 bool
88 mx_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 -*-
134 Matrix
135 xdiv (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,
148 
149  typ = btyp.transpose ();
150  return result.transpose ();
151 }
152 
153 // -*- 2 -*-
155 xdiv (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 -*-
175 xdiv (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 -*-
215 xdiv (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,
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;
246  SparseComplexMatrix result
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;
266  SparseComplexMatrix result
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;
287  SparseComplexMatrix result
288  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
289 
290  typ = btyp.transpose ();
291  return result.hermitian ();
292 }
293 
294 template <typename RT, typename SM, typename DM>
295 RT do_rightdiv_sm_dm (const SM& a, const DM& d)
296 {
297  const octave_idx_type d_nr = d.rows ();
298 
299  const octave_idx_type a_nr = a.rows ();
300  const octave_idx_type a_nc = a.cols ();
301 
302  using std::min;
303  const octave_idx_type nc = min (d_nr, a_nc);
304 
305  if (! mx_div_conform (a, d))
306  return RT ();
307 
308  const octave_idx_type nz = a.nnz ();
309  RT r (a_nr, nc, nz);
310 
311  typedef typename DM::element_type DM_elt_type;
312  const DM_elt_type zero = DM_elt_type ();
313 
314  octave_idx_type k_result = 0;
315  for (octave_idx_type j = 0; j < nc; ++j)
316  {
317  octave_quit ();
318  const DM_elt_type s = d.dgelem (j);
319  const octave_idx_type colend = a.cidx (j+1);
320  r.xcidx (j) = k_result;
321  if (s != zero)
322  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
323  {
324  r.xdata (k_result) = a.data (k) / s;
325  r.xridx (k_result) = a.ridx (k);
326  ++k_result;
327  }
328  }
329  r.xcidx (nc) = k_result;
330 
331  r.maybe_compress (true);
332  return r;
333 }
334 
335 // -*- 9 -*-
337 xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType &)
338 {
339  return do_rightdiv_sm_dm<SparseMatrix> (a, b);
340 }
341 
342 // -*- 10 -*-
345 {
346  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
347 }
348 
349 // -*- 11 -*-
352 {
353  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
354 }
355 
356 // -*- 12 -*-
359 {
360  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
361 }
362 
363 // Funny element by element division operations.
364 //
365 // op2 \ op1: s cs
366 // +-- +---+----+
367 // matrix | 1 | 3 |
368 // +---+----+
369 // complex_matrix | 2 | 4 |
370 // +---+----+
371 
372 Matrix
373 x_el_div (double a, const SparseMatrix& b)
374 {
375  octave_idx_type nr = b.rows ();
376  octave_idx_type nc = b.cols ();
377 
378  Matrix result;
379  if (a == 0.)
380  result = Matrix (nr, nc, octave::numeric_limits<double>::NaN ());
381  else if (a > 0.)
382  result = Matrix (nr, nc, octave::numeric_limits<double>::Inf ());
383  else
384  result = Matrix (nr, nc, -octave::numeric_limits<double>::Inf ());
385 
386  for (octave_idx_type j = 0; j < nc; j++)
387  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
388  {
389  octave_quit ();
390  result.elem (b.ridx (i), j) = a / b.data (i);
391  }
392 
393  return result;
394 }
395 
397 x_el_div (double a, const SparseComplexMatrix& b)
398 {
399  octave_idx_type nr = b.rows ();
400  octave_idx_type nc = b.cols ();
401 
404 
405  for (octave_idx_type j = 0; j < nc; j++)
406  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
407  {
408  octave_quit ();
409  result.elem (b.ridx (i), j) = a / b.data (i);
410  }
411 
412  return result;
413 }
414 
416 x_el_div (const Complex a, const SparseMatrix& b)
417 {
418  octave_idx_type nr = b.rows ();
419  octave_idx_type nc = b.cols ();
420 
421  ComplexMatrix result (nr, nc, (a / 0.0));
422 
423  for (octave_idx_type j = 0; j < nc; j++)
424  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
425  {
426  octave_quit ();
427  result.elem (b.ridx (i), j) = a / b.data (i);
428  }
429 
430  return result;
431 }
432 
435 {
436  octave_idx_type nr = b.rows ();
437  octave_idx_type nc = b.cols ();
438 
439  ComplexMatrix result (nr, nc, (a / 0.0));
440 
441  for (octave_idx_type j = 0; j < nc; j++)
442  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
443  {
444  octave_quit ();
445  result.elem (b.ridx (i), j) = a / b.data (i);
446  }
447 
448  return result;
449 }
450 
451 // Left division functions. X \ Y = inv (X) * Y
452 //
453 // Y \ X : sm scm dm dcm
454 // +-- +---+----+
455 // matrix | 1 | 5 |
456 // +---+----+
457 // complex_matrix | 2 | 6 |
458 // +---+----+----+----+
459 // sparse matrix | 3 | 7 | 9 | 11 |
460 // +---+----+----+----+
461 // sparse complex_matrix | 4 | 8 | 10 | 12 |
462 // +---+----+----+----+
463 
464 // -*- 1 -*-
465 Matrix
466 xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ)
467 {
468  if (! mx_leftdiv_conform (a, b))
469  return Matrix ();
470 
471  octave_idx_type info;
472  double rcond = 0.0;
473  return a.solve (typ, b, info, rcond, solve_singularity_warning);
474 }
475 
476 // -*- 2 -*-
478 xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, MatrixType& typ)
479 {
480  if (! mx_leftdiv_conform (a, b))
481  return ComplexMatrix ();
482 
483  octave_idx_type info;
484  double rcond = 0.0;
485  return a.solve (typ, b, info, rcond, solve_singularity_warning);
486 }
487 
488 // -*- 3 -*-
490 xleftdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
491 {
492  if (! mx_leftdiv_conform (a, b))
493  return SparseMatrix ();
494 
495  octave_idx_type info;
496  double rcond = 0.0;
497  return a.solve (typ, b, info, rcond, solve_singularity_warning);
498 }
499 
500 // -*- 4 -*-
503 {
504  if (! mx_leftdiv_conform (a, b))
505  return SparseComplexMatrix ();
506 
507  octave_idx_type info;
508  double rcond = 0.0;
509  return a.solve (typ, b, info, rcond, solve_singularity_warning);
510 }
511 
512 // -*- 5 -*-
514 xleftdiv (const SparseComplexMatrix& a, const Matrix& b, MatrixType& typ)
515 {
516  if (! mx_leftdiv_conform (a, b))
517  return ComplexMatrix ();
518 
519  octave_idx_type info;
520  double rcond = 0.0;
521  return a.solve (typ, b, info, rcond, solve_singularity_warning);
522 }
523 
524 // -*- 6 -*-
527 {
528  if (! mx_leftdiv_conform (a, b))
529  return ComplexMatrix ();
530 
531  octave_idx_type info;
532  double rcond = 0.0;
533  return a.solve (typ, b, info, rcond, solve_singularity_warning);
534 }
535 
536 // -*- 7 -*-
539 {
540  if (! mx_leftdiv_conform (a, b))
541  return SparseComplexMatrix ();
542 
543  octave_idx_type info;
544  double rcond = 0.0;
545  return a.solve (typ, b, info, rcond, solve_singularity_warning);
546 }
547 
548 // -*- 8 -*-
551  MatrixType& typ)
552 {
553  if (! mx_leftdiv_conform (a, b))
554  return SparseComplexMatrix ();
555 
556  octave_idx_type info;
557  double rcond = 0.0;
558  return a.solve (typ, b, info, rcond, solve_singularity_warning);
559 }
560 
561 template <typename RT, typename DM, typename SM>
562 RT do_leftdiv_dm_sm (const DM& d, const SM& a)
563 {
564  const octave_idx_type a_nr = a.rows ();
565  const octave_idx_type a_nc = a.cols ();
566 
567  const octave_idx_type d_nc = d.cols ();
568 
569  using std::min;
570  const octave_idx_type nr = min (d_nc, a_nr);
571 
572  if (! mx_leftdiv_conform (d, a))
573  return RT ();
574 
575  const octave_idx_type nz = a.nnz ();
576  RT r (nr, a_nc, nz);
577 
578  typedef typename DM::element_type DM_elt_type;
579  const DM_elt_type zero = DM_elt_type ();
580 
581  octave_idx_type k_result = 0;
582  for (octave_idx_type j = 0; j < a_nc; ++j)
583  {
584  octave_quit ();
585  const octave_idx_type colend = a.cidx (j+1);
586  r.xcidx (j) = k_result;
587  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
588  {
589  const octave_idx_type i = a.ridx (k);
590  if (i < nr)
591  {
592  const DM_elt_type s = d.dgelem (i);
593  if (s != zero)
594  {
595  r.xdata (k_result) = a.data (k) / s;
596  r.xridx (k_result) = i;
597  ++k_result;
598  }
599  }
600  }
601  }
602  r.xcidx (a_nc) = k_result;
603 
604  r.maybe_compress (true);
605  return r;
606 }
607 
608 // -*- 9 -*-
611 {
612  return do_leftdiv_dm_sm<SparseMatrix> (d, a);
613 }
614 
615 // -*- 10 -*-
618 {
619  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
620 }
621 
622 // -*- 11 -*-
625 {
626  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
627 }
628 
629 // -*- 12 -*-
632  MatrixType&)
633 {
634  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
635 }
#define Inf
Definition: Faddeeva.cc:247
#define NaN
Definition: Faddeeva.cc:248
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:499
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:169
MatrixType transpose(void) const
Definition: MatrixType.cc:966
Definition: dMatrix.h:42
Matrix transpose(void) const
Definition: dMatrix.h:135
SparseComplexMatrix hermitian(void) const
Definition: CSparse.cc:606
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition: CSparse.cc:6734
SparseMatrix transpose(void) const
Definition: dSparse.h:128
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6725
octave_idx_type cols(void) const
Definition: Sparse.h:251
T * data(void)
Definition: Sparse.h:470
octave_idx_type * cidx(void)
Definition: Sparse.h:492
octave_idx_type rows(void) const
Definition: Sparse.h:250
octave_idx_type * ridx(void)
Definition: Sparse.h:479
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
T * r
Definition: mx-inlines.cc:773
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
Matrix x_el_div(double a, const SparseMatrix &b)
Definition: sparse-xdiv.cc:373
Matrix xdiv(const Matrix &a, const SparseMatrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:135
bool mx_leftdiv_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:54
bool mx_div_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:88
RT do_rightdiv_sm_dm(const SM &a, const DM &d)
Definition: sparse-xdiv.cc:295
static void solve_singularity_warning(double rcond)
Definition: sparse-xdiv.cc:47
#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:70
#define INSTANTIATE_MX_DIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:104
Matrix xleftdiv(const SparseMatrix &a, const Matrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:466
RT do_leftdiv_dm_sm(const DM &d, const SM &a)
Definition: sparse-xdiv.cc:562