GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sparse-xdiv.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2015 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 the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 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 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <cassert>
29 
30 #include "Array-util.h"
31 #include "lo-array-gripes.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 {
47  gripe_singular_matrix (rcond);
48 }
49 
50 template <class T1, class 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  gripe_nonconformant ("operator \\", a_nr, a_nc, b_nr, b_nc);
63  return false;
64  }
65 
66  return true;
67 }
68 
69 #define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
70  template bool mx_leftdiv_conform (const T1&, const T2&)
71 
84 
85 template <class T1, class T2>
86 bool
87 mx_div_conform (const T1& a, const T2& b)
88 {
89  octave_idx_type a_nc = a.cols ();
90  octave_idx_type b_nc = b.cols ();
91 
92  if (a_nc != b_nc)
93  {
94  octave_idx_type a_nr = a.rows ();
95  octave_idx_type b_nr = b.rows ();
96 
97  gripe_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
98  return false;
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_NaN);
381  else if (a > 0.)
382  result = Matrix (nr, nc, octave_Inf);
383  else
384  result = Matrix (nr, nc, -octave_Inf);
385 
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 
398 x_el_div (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_NaN, octave_NaN));
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 }
octave_idx_type cols(void) const
Definition: Sparse.h:264
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
T * data(void)
Definition: Sparse.h:509
octave_idx_type rows(void) const
Definition: Sparse.h:263
static void solve_singularity_warning(double rcond)
Definition: sparse-xdiv.cc:45
ComplexMatrix solve(MatrixType &typ, const Matrix &b) const
Definition: CSparse.cc:6595
SparseMatrix transpose(void) const
Definition: dSparse.h:133
octave_idx_type * cidx(void)
Definition: Sparse.h:531
bool mx_leftdiv_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:52
T & elem(octave_idx_type n)
Definition: Array.h:380
Matrix solve(MatrixType &typ, const Matrix &b) const
Definition: dSparse.cc:6735
F77_RET_T const double const double double * d
ComplexMatrix hermitian(void) const
Definition: CMatrix.h:149
#define octave_Inf
Definition: lo-ieee.h:31
Matrix x_el_div(double a, const SparseMatrix &b)
Definition: sparse-xdiv.cc:373
MatrixType transpose(void) const
Definition: MatrixType.cc:1274
#define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:69
Matrix transpose(void) const
Definition: dMatrix.h:114
Matrix xleftdiv(const SparseMatrix &a, const Matrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:466
void gripe_singular_matrix(double rcond)
SparseComplexMatrix hermitian(void) const
Definition: CSparse.cc:679
Definition: dMatrix.h:35
bool mx_div_conform(const T1 &a, const T2 &b)
Definition: sparse-xdiv.cc:87
octave_idx_type * ridx(void)
Definition: Sparse.h:518
Matrix xdiv(const Matrix &a, const SparseMatrix &b, MatrixType &typ)
Definition: sparse-xdiv.cc:135
#define octave_NaN
Definition: lo-ieee.h:37
#define INSTANTIATE_MX_DIV_CONFORM(T1, T2)
Definition: sparse-xdiv.cc:104
std::complex< double > Complex
Definition: oct-cmplx.h:29
RT do_leftdiv_dm_sm(const DM &d, const SM &a)
Definition: sparse-xdiv.cc:562
RT do_rightdiv_sm_dm(const SM &a, const DM &d)
Definition: sparse-xdiv.cc:295
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210