GNU Octave  9.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-2024 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
49 solve_singularity_warning (double rcond)
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,
149  solve_singularity_warning);
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,
229  solve_singularity_warning);
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
298 do_rightdiv_sm_dm (const SM& a, const DM& d)
299 {
300  const octave_idx_type d_nr = d.rows ();
301 
302  const octave_idx_type a_nr = a.rows ();
303  const octave_idx_type a_nc = a.cols ();
304 
305  using std::min;
306  const octave_idx_type nc = min (d_nr, a_nc);
307 
308  if (! mx_div_conform (a, d))
309  return RT ();
310 
311  const octave_idx_type nz = a.nnz ();
312  RT r (a_nr, nc, nz);
313 
314  typedef typename DM::element_type DM_elt_type;
315  const DM_elt_type zero = DM_elt_type ();
316 
317  octave_idx_type k_result = 0;
318  for (octave_idx_type j = 0; j < nc; ++j)
319  {
320  octave_quit ();
321  const DM_elt_type s = d.dgelem (j);
322  const octave_idx_type colend = a.cidx (j+1);
323  r.xcidx (j) = k_result;
324  if (s != zero)
325  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
326  {
327  r.xdata (k_result) = a.data (k) / s;
328  r.xridx (k_result) = a.ridx (k);
329  ++k_result;
330  }
331  }
332  r.xcidx (nc) = k_result;
333 
334  r.maybe_compress (true);
335  return r;
336 }
337 
338 // -*- 9 -*-
340 xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType&)
341 {
342  return do_rightdiv_sm_dm<SparseMatrix> (a, b);
343 }
344 
345 // -*- 10 -*-
348 {
349  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
350 }
351 
352 // -*- 11 -*-
355 {
356  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
357 }
358 
359 // -*- 12 -*-
362 {
363  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
364 }
365 
366 // Funny element by element division operations.
367 //
368 // op2 \ op1: s cs
369 // +-- +---+----+
370 // matrix | 1 | 3 |
371 // +---+----+
372 // complex_matrix | 2 | 4 |
373 // +---+----+
374 
375 Matrix
376 elem_xdiv (double a, const SparseMatrix& b)
377 {
378  octave_idx_type nr = b.rows ();
379  octave_idx_type nc = b.cols ();
380 
381  Matrix result;
382  if (a == 0.)
383  result = Matrix (nr, nc, octave::numeric_limits<double>::NaN ());
384  else if (a > 0.)
385  result = Matrix (nr, nc, octave::numeric_limits<double>::Inf ());
386  else
387  result = Matrix (nr, nc, -octave::numeric_limits<double>::Inf ());
388 
389  for (octave_idx_type j = 0; j < nc; j++)
390  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
391  {
392  octave_quit ();
393  result.elem (b.ridx (i), j) = a / b.data (i);
394  }
395 
396  return result;
397 }
398 
400 elem_xdiv (double a, const SparseComplexMatrix& b)
401 {
402  octave_idx_type nr = b.rows ();
403  octave_idx_type nc = b.cols ();
404 
407 
408  for (octave_idx_type j = 0; j < nc; j++)
409  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
410  {
411  octave_quit ();
412  result.elem (b.ridx (i), j) = a / b.data (i);
413  }
414 
415  return result;
416 }
417 
419 elem_xdiv (const Complex& a, const SparseMatrix& b)
420 {
421  octave_idx_type nr = b.rows ();
422  octave_idx_type nc = b.cols ();
423 
424  ComplexMatrix result (nr, nc, (a / 0.0));
425 
426  for (octave_idx_type j = 0; j < nc; j++)
427  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
428  {
429  octave_quit ();
430  result.elem (b.ridx (i), j) = a / b.data (i);
431  }
432 
433  return result;
434 }
435 
438 {
439  octave_idx_type nr = b.rows ();
440  octave_idx_type nc = b.cols ();
441 
442  ComplexMatrix result (nr, nc, (a / 0.0));
443 
444  for (octave_idx_type j = 0; j < nc; j++)
445  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
446  {
447  octave_quit ();
448  result.elem (b.ridx (i), j) = a / b.data (i);
449  }
450 
451  return result;
452 }
453 
454 // Left division functions. X \ Y = inv (X) * Y
455 //
456 // Y \ X : sm scm dm dcm
457 // +-- +---+----+
458 // matrix | 1 | 5 |
459 // +---+----+
460 // complex_matrix | 2 | 6 |
461 // +---+----+----+----+
462 // sparse matrix | 3 | 7 | 9 | 11 |
463 // +---+----+----+----+
464 // sparse complex_matrix | 4 | 8 | 10 | 12 |
465 // +---+----+----+----+
466 
467 // -*- 1 -*-
468 Matrix
469 xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ)
470 {
471  if (! mx_leftdiv_conform (a, b))
472  return Matrix ();
473 
474  octave_idx_type info;
475  double rcond = 0.0;
476  return a.solve (typ, b, info, rcond, solve_singularity_warning);
477 }
478 
479 // -*- 2 -*-
481 xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, MatrixType& typ)
482 {
483  if (! mx_leftdiv_conform (a, b))
484  return ComplexMatrix ();
485 
486  octave_idx_type info;
487  double rcond = 0.0;
488  return a.solve (typ, b, info, rcond, solve_singularity_warning);
489 }
490 
491 // -*- 3 -*-
493 xleftdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
494 {
495  if (! mx_leftdiv_conform (a, b))
496  return SparseMatrix ();
497 
498  octave_idx_type info;
499  double rcond = 0.0;
500  return a.solve (typ, b, info, rcond, solve_singularity_warning);
501 }
502 
503 // -*- 4 -*-
506 {
507  if (! mx_leftdiv_conform (a, b))
508  return SparseComplexMatrix ();
509 
510  octave_idx_type info;
511  double rcond = 0.0;
512  return a.solve (typ, b, info, rcond, solve_singularity_warning);
513 }
514 
515 // -*- 5 -*-
517 xleftdiv (const SparseComplexMatrix& a, const Matrix& b, MatrixType& typ)
518 {
519  if (! mx_leftdiv_conform (a, b))
520  return ComplexMatrix ();
521 
522  octave_idx_type info;
523  double rcond = 0.0;
524  return a.solve (typ, b, info, rcond, solve_singularity_warning);
525 }
526 
527 // -*- 6 -*-
530 {
531  if (! mx_leftdiv_conform (a, b))
532  return ComplexMatrix ();
533 
534  octave_idx_type info;
535  double rcond = 0.0;
536  return a.solve (typ, b, info, rcond, solve_singularity_warning);
537 }
538 
539 // -*- 7 -*-
542 {
543  if (! mx_leftdiv_conform (a, b))
544  return SparseComplexMatrix ();
545 
546  octave_idx_type info;
547  double rcond = 0.0;
548  return a.solve (typ, b, info, rcond, solve_singularity_warning);
549 }
550 
551 // -*- 8 -*-
554  MatrixType& typ)
555 {
556  if (! mx_leftdiv_conform (a, b))
557  return SparseComplexMatrix ();
558 
559  octave_idx_type info;
560  double rcond = 0.0;
561  return a.solve (typ, b, info, rcond, solve_singularity_warning);
562 }
563 
564 template <typename RT, typename DM, typename SM>
565 RT
566 do_leftdiv_dm_sm (const DM& d, const SM& a)
567 {
568  const octave_idx_type a_nr = a.rows ();
569  const octave_idx_type a_nc = a.cols ();
570 
571  const octave_idx_type d_nc = d.cols ();
572 
573  using std::min;
574  const octave_idx_type nr = min (d_nc, a_nr);
575 
576  if (! mx_leftdiv_conform (d, a))
577  return RT ();
578 
579  const octave_idx_type nz = a.nnz ();
580  RT r (nr, a_nc, nz);
581 
582  typedef typename DM::element_type DM_elt_type;
583  const DM_elt_type zero = DM_elt_type ();
584 
585  octave_idx_type k_result = 0;
586  for (octave_idx_type j = 0; j < a_nc; ++j)
587  {
588  octave_quit ();
589  const octave_idx_type colend = a.cidx (j+1);
590  r.xcidx (j) = k_result;
591  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
592  {
593  const octave_idx_type i = a.ridx (k);
594  if (i < nr)
595  {
596  const DM_elt_type s = d.dgelem (i);
597  if (s != zero)
598  {
599  r.xdata (k_result) = a.data (k) / s;
600  r.xridx (k_result) = i;
601  ++k_result;
602  }
603  }
604  }
605  }
606  r.xcidx (a_nc) = k_result;
607 
608  r.maybe_compress (true);
609  return r;
610 }
611 
612 // -*- 9 -*-
615 {
616  return do_leftdiv_dm_sm<SparseMatrix> (d, a);
617 }
618 
619 // -*- 10 -*-
622 {
623  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
624 }
625 
626 // -*- 11 -*-
629 {
630  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
631 }
632 
633 // -*- 12 -*-
636  MatrixType&)
637 {
638  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
639 }
640 
641 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
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:562
ComplexMatrix hermitian() const
Definition: CMatrix.h:170
MatrixType transpose() const
Definition: MatrixType.cc:973
Definition: dMatrix.h:42
Matrix transpose() const
Definition: dMatrix.h:140
SparseComplexMatrix hermitian() const
Definition: CSparse.cc:606
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition: CSparse.cc:6706
SparseMatrix transpose() const
Definition: dSparse.h:126
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6697
octave_idx_type cols() const
Definition: Sparse.h:352
octave_idx_type * cidx()
Definition: Sparse.h:596
T * data()
Definition: Sparse.h:574
octave_idx_type * ridx()
Definition: Sparse.h:583
octave_idx_type rows() const
Definition: Sparse.h:351
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
T * r
Definition: mx-inlines.cc:781
std::complex< double > Complex
Definition: oct-cmplx.h:33
Matrix elem_xdiv(double a, const SparseMatrix &b)
Definition: sparse-xdiv.cc:376
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:298
#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:469
RT do_leftdiv_dm_sm(const DM &d, const SM &a)
Definition: sparse-xdiv.cc:566