GNU Octave  4.4.1
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) 2004-2018 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
5 
6 This file is part of Octave.
7 
8 Octave is free software: you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <cassert>
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 
44 static void
46 {
48 }
49 
50 template <typename T1, typename T2>
51 bool
52 mx_leftdiv_conform (const T1& a, const T2& b)
53 {
54  octave_idx_type a_nr = a.rows ();
55  octave_idx_type b_nr = b.rows ();
56 
57  if (a_nr != b_nr)
58  {
59  octave_idx_type a_nc = a.cols ();
60  octave_idx_type b_nc = b.cols ();
61 
62  octave::err_nonconformant (R"(operator \)", a_nr, a_nc, b_nr, b_nc);
63  }
64 
65  return true;
66 }
67 
68 #define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
69  template bool mx_leftdiv_conform (const T1&, const T2&)
70 
83 
84 template <typename T1, typename T2>
85 bool
86 mx_div_conform (const T1& a, const T2& b)
87 {
88  octave_idx_type a_nc = a.cols ();
89  octave_idx_type b_nc = b.cols ();
90 
91  if (a_nc != b_nc)
92  {
93  octave_idx_type a_nr = a.rows ();
94  octave_idx_type b_nr = b.rows ();
95 
96  octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
97  }
98 
99  return true;
100 }
101 
102 #define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \
103  template bool mx_div_conform (const T1&, const T2&)
104 
117 
118 // Right division functions. X / Y = X * inv (Y) = (inv (Y') * X')'
119 //
120 // Y / X: m cm sm scm
121 // +-- +---+----+----+----+
122 // sparse matrix | 1 | 3 | 5 | 7 |
123 // +---+----+----+----+
124 // sparse complex_matrix | 2 | 4 | 6 | 8 |
125 // +---+----+----+----+
126 // diagonal matrix | 9 | 11 |
127 // +----+----+
128 // complex diag. matrix | 10 | 12 |
129 // +----+----+
130 
131 // -*- 1 -*-
132 Matrix
133 xdiv (const Matrix& a, const SparseMatrix& b, MatrixType& typ)
134 {
135  if (! mx_div_conform (a, b))
136  return Matrix ();
137 
138  Matrix atmp = a.transpose ();
139  SparseMatrix btmp = b.transpose ();
140  MatrixType btyp = typ.transpose ();
141 
142  octave_idx_type info;
143  double rcond = 0.0;
144  Matrix result = btmp.solve (btyp, atmp, info, rcond,
146 
147  typ = btyp.transpose ();
148  return result.transpose ();
149 }
150 
151 // -*- 2 -*-
153 xdiv (const Matrix& a, const SparseComplexMatrix& b, MatrixType& typ)
154 {
155  if (! mx_div_conform (a, b))
156  return ComplexMatrix ();
157 
158  Matrix atmp = a.transpose ();
159  SparseComplexMatrix btmp = b.hermitian ();
160  MatrixType btyp = typ.transpose ();
161 
162  octave_idx_type info;
163  double rcond = 0.0;
165  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
166 
167  typ = btyp.transpose ();
168  return result.hermitian ();
169 }
170 
171 // -*- 3 -*-
173 xdiv (const ComplexMatrix& a, const SparseMatrix& b, MatrixType& typ)
174 {
175  if (! mx_div_conform (a, b))
176  return ComplexMatrix ();
177 
178  ComplexMatrix atmp = a.hermitian ();
179  SparseMatrix btmp = b.transpose ();
180  MatrixType btyp = typ.transpose ();
181 
182  octave_idx_type info;
183  double rcond = 0.0;
185  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
186 
187  typ = btyp.transpose ();
188  return result.hermitian ();
189 }
190 
191 // -*- 4 -*-
194 {
195  if (! mx_div_conform (a, b))
196  return ComplexMatrix ();
197 
198  ComplexMatrix atmp = a.hermitian ();
199  SparseComplexMatrix btmp = b.hermitian ();
200  MatrixType btyp = typ.transpose ();
201 
202  octave_idx_type info;
203  double rcond = 0.0;
205  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
206 
207  typ = btyp.transpose ();
208  return result.hermitian ();
209 }
210 
211 // -*- 5 -*-
213 xdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ)
214 {
215  if (! mx_div_conform (a, b))
216  return SparseMatrix ();
217 
218  SparseMatrix atmp = a.transpose ();
219  SparseMatrix btmp = b.transpose ();
220  MatrixType btyp = typ.transpose ();
221 
222  octave_idx_type info;
223  double rcond = 0.0;
224  SparseMatrix result = btmp.solve (btyp, atmp, info, rcond,
226 
227  typ = btyp.transpose ();
228  return result.transpose ();
229 }
230 
231 // -*- 6 -*-
234 {
235  if (! mx_div_conform (a, b))
236  return SparseComplexMatrix ();
237 
238  SparseMatrix atmp = a.transpose ();
239  SparseComplexMatrix btmp = b.hermitian ();
240  MatrixType btyp = typ.transpose ();
241 
242  octave_idx_type info;
243  double rcond = 0.0;
245  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
246 
247  typ = btyp.transpose ();
248  return result.hermitian ();
249 }
250 
251 // -*- 7 -*-
254 {
255  if (! mx_div_conform (a, b))
256  return SparseComplexMatrix ();
257 
258  SparseComplexMatrix atmp = a.hermitian ();
259  SparseMatrix btmp = b.transpose ();
260  MatrixType btyp = typ.transpose ();
261 
262  octave_idx_type info;
263  double rcond = 0.0;
265  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
266 
267  typ = btyp.transpose ();
268  return result.hermitian ();
269 }
270 
271 // -*- 8 -*-
274  MatrixType& typ)
275 {
276  if (! mx_div_conform (a, b))
277  return SparseComplexMatrix ();
278 
279  SparseComplexMatrix atmp = a.hermitian ();
280  SparseComplexMatrix btmp = b.hermitian ();
281  MatrixType btyp = typ.transpose ();
282 
283  octave_idx_type info;
284  double rcond = 0.0;
286  = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
287 
288  typ = btyp.transpose ();
289  return result.hermitian ();
290 }
291 
292 template <typename RT, typename SM, typename DM>
293 RT do_rightdiv_sm_dm (const SM& a, const DM& d)
294 {
295  const octave_idx_type d_nr = d.rows ();
296 
297  const octave_idx_type a_nr = a.rows ();
298  const octave_idx_type a_nc = a.cols ();
299 
300  using std::min;
301  const octave_idx_type nc = min (d_nr, a_nc);
302 
303  if (! mx_div_conform (a, d))
304  return RT ();
305 
306  const octave_idx_type nz = a.nnz ();
307  RT r (a_nr, nc, nz);
308 
309  typedef typename DM::element_type DM_elt_type;
310  const DM_elt_type zero = DM_elt_type ();
311 
312  octave_idx_type k_result = 0;
313  for (octave_idx_type j = 0; j < nc; ++j)
314  {
315  octave_quit ();
316  const DM_elt_type s = d.dgelem (j);
317  const octave_idx_type colend = a.cidx (j+1);
318  r.xcidx (j) = k_result;
319  if (s != zero)
320  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
321  {
322  r.xdata (k_result) = a.data (k) / s;
323  r.xridx (k_result) = a.ridx (k);
324  ++k_result;
325  }
326  }
327  r.xcidx (nc) = k_result;
328 
329  r.maybe_compress (true);
330  return r;
331 }
332 
333 // -*- 9 -*-
336 {
337  return do_rightdiv_sm_dm<SparseMatrix> (a, b);
338 }
339 
340 // -*- 10 -*-
343 {
344  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
345 }
346 
347 // -*- 11 -*-
350 {
351  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
352 }
353 
354 // -*- 12 -*-
357 {
358  return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
359 }
360 
361 // Funny element by element division operations.
362 //
363 // op2 \ op1: s cs
364 // +-- +---+----+
365 // matrix | 1 | 3 |
366 // +---+----+
367 // complex_matrix | 2 | 4 |
368 // +---+----+
369 
370 Matrix
371 x_el_div (double a, const SparseMatrix& b)
372 {
373  octave_idx_type nr = b.rows ();
374  octave_idx_type nc = b.cols ();
375 
376  Matrix result;
377  if (a == 0.)
379  else if (a > 0.)
381  else
383 
384  for (octave_idx_type j = 0; j < nc; j++)
385  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
386  {
387  octave_quit ();
388  result.elem (b.ridx (i), j) = a / b.data (i);
389  }
390 
391  return result;
392 }
393 
395 x_el_div (double a, const SparseComplexMatrix& b)
396 {
397  octave_idx_type nr = b.rows ();
398  octave_idx_type nc = b.cols ();
399 
402 
403  for (octave_idx_type j = 0; j < nc; j++)
404  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
405  {
406  octave_quit ();
407  result.elem (b.ridx (i), j) = a / b.data (i);
408  }
409 
410  return result;
411 }
412 
414 x_el_div (const Complex a, const SparseMatrix& b)
415 {
416  octave_idx_type nr = b.rows ();
417  octave_idx_type nc = b.cols ();
418 
419  ComplexMatrix result (nr, nc, (a / 0.0));
420 
421  for (octave_idx_type j = 0; j < nc; j++)
422  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
423  {
424  octave_quit ();
425  result.elem (b.ridx (i), j) = a / b.data (i);
426  }
427 
428  return result;
429 }
430 
433 {
434  octave_idx_type nr = b.rows ();
435  octave_idx_type nc = b.cols ();
436 
437  ComplexMatrix result (nr, nc, (a / 0.0));
438 
439  for (octave_idx_type j = 0; j < nc; j++)
440  for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
441  {
442  octave_quit ();
443  result.elem (b.ridx (i), j) = a / b.data (i);
444  }
445 
446  return result;
447 }
448 
449 // Left division functions. X \ Y = inv (X) * Y
450 //
451 // Y \ X : sm scm dm dcm
452 // +-- +---+----+
453 // matrix | 1 | 5 |
454 // +---+----+
455 // complex_matrix | 2 | 6 |
456 // +---+----+----+----+
457 // sparse matrix | 3 | 7 | 9 | 11 |
458 // +---+----+----+----+
459 // sparse complex_matrix | 4 | 8 | 10 | 12 |
460 // +---+----+----+----+
461 
462 // -*- 1 -*-
463 Matrix
464 xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ)
465 {
466  if (! mx_leftdiv_conform (a, b))
467  return Matrix ();
468 
469  octave_idx_type info;
470  double rcond = 0.0;
471  return a.solve (typ, b, info, rcond, solve_singularity_warning);
472 }
473 
474 // -*- 2 -*-
477 {
478  if (! mx_leftdiv_conform (a, b))
479  return ComplexMatrix ();
480 
481  octave_idx_type info;
482  double rcond = 0.0;
483  return a.solve (typ, b, info, rcond, solve_singularity_warning);
484 }
485 
486 // -*- 3 -*-
489 {
490  if (! mx_leftdiv_conform (a, b))
491  return SparseMatrix ();
492 
493  octave_idx_type info;
494  double rcond = 0.0;
495  return a.solve (typ, b, info, rcond, solve_singularity_warning);
496 }
497 
498 // -*- 4 -*-
501 {
502  if (! mx_leftdiv_conform (a, b))
503  return SparseComplexMatrix ();
504 
505  octave_idx_type info;
506  double rcond = 0.0;
507  return a.solve (typ, b, info, rcond, solve_singularity_warning);
508 }
509 
510 // -*- 5 -*-
513 {
514  if (! mx_leftdiv_conform (a, b))
515  return ComplexMatrix ();
516 
517  octave_idx_type info;
518  double rcond = 0.0;
519  return a.solve (typ, b, info, rcond, solve_singularity_warning);
520 }
521 
522 // -*- 6 -*-
525 {
526  if (! mx_leftdiv_conform (a, b))
527  return ComplexMatrix ();
528 
529  octave_idx_type info;
530  double rcond = 0.0;
531  return a.solve (typ, b, info, rcond, solve_singularity_warning);
532 }
533 
534 // -*- 7 -*-
537 {
538  if (! mx_leftdiv_conform (a, b))
539  return SparseComplexMatrix ();
540 
541  octave_idx_type info;
542  double rcond = 0.0;
543  return a.solve (typ, b, info, rcond, solve_singularity_warning);
544 }
545 
546 // -*- 8 -*-
549  MatrixType& typ)
550 {
551  if (! mx_leftdiv_conform (a, b))
552  return SparseComplexMatrix ();
553 
554  octave_idx_type info;
555  double rcond = 0.0;
556  return a.solve (typ, b, info, rcond, solve_singularity_warning);
557 }
558 
559 template <typename RT, typename DM, typename SM>
560 RT do_leftdiv_dm_sm (const DM& d, const SM& a)
561 {
562  const octave_idx_type a_nr = a.rows ();
563  const octave_idx_type a_nc = a.cols ();
564 
565  const octave_idx_type d_nc = d.cols ();
566 
567  using std::min;
568  const octave_idx_type nr = min (d_nc, a_nr);
569 
570  if (! mx_leftdiv_conform (d, a))
571  return RT ();
572 
573  const octave_idx_type nz = a.nnz ();
574  RT r (nr, a_nc, nz);
575 
576  typedef typename DM::element_type DM_elt_type;
577  const DM_elt_type zero = DM_elt_type ();
578 
579  octave_idx_type k_result = 0;
580  for (octave_idx_type j = 0; j < a_nc; ++j)
581  {
582  octave_quit ();
583  const octave_idx_type colend = a.cidx (j+1);
584  r.xcidx (j) = k_result;
585  for (octave_idx_type k = a.cidx (j); k < colend; ++k)
586  {
587  const octave_idx_type i = a.ridx (k);
588  if (i < nr)
589  {
590  const DM_elt_type s = d.dgelem (i);
591  if (s != zero)
592  {
593  r.xdata (k_result) = a.data (k) / s;
594  r.xridx (k_result) = i;
595  ++k_result;
596  }
597  }
598  }
599  }
600  r.xcidx (a_nc) = k_result;
601 
602  r.maybe_compress (true);
603  return r;
604 }
605 
606 // -*- 9 -*-
609 {
610  return do_leftdiv_dm_sm<SparseMatrix> (d, a);
611 }
612 
613 // -*- 10 -*-
616 {
617  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
618 }
619 
620 // -*- 11 -*-
623 {
624  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
625 }
626 
627 // -*- 12 -*-
630  MatrixType&)
631 {
632  return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
633 }
static void solve_singularity_warning(double rcond)
Definition: sparse-xdiv.cc:45
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by zero($0/0$)
for large enough k
Definition: lu.cc:617
octave_idx_type b_nr
Definition: sylvester.cc:76
bool mx_leftdiv_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:52
s
Definition: file-io.cc:2729
octave_idx_type a_nc
Definition: sylvester.cc:74
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
octave_idx_type a_nr
Definition: sylvester.cc:73
Matrix x_el_div(double a, const SparseMatrix &b)
Definition: sparse-xdiv.cc:371
#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:68
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
Matrix xleftdiv(const SparseMatrix &a, const Matrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:464
Definition: dMatrix.h:36
With real return the complex result
Definition: data.cc:3260
bool mx_div_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:86
#define Inf
Definition: Faddeeva.cc:247
Matrix xdiv(const Matrix &a, const SparseMatrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:133
octave_idx_type b_nc
Definition: sylvester.cc:77
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol NaN(Not a Number). NaN is the result of operations which do not produce a well defined 0 result. Common operations which produce a NaN are arithmetic with infinity ex($\infty - \infty$)
ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition: CSparse.cc:6735
b
Definition: cellfun.cc:400
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6725
#define INSTANTIATE_MX_DIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:102
for i
Definition: data.cc:5264
MatrixType transpose(void) const
Definition: MatrixType.cc:961
std::complex< double > Complex
Definition: oct-cmplx.h:31
RT do_leftdiv_dm_sm(const DM &d, const SM &a)
Definition: sparse-xdiv.cc:560
RT do_rightdiv_sm_dm(const SM &a, const DM &d)
Definition: sparse-xdiv.cc:293
void warn_singular_matrix(double rcond)
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204