1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1993-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 "CMatrix.h"
34 #include "dMatrix.h"
35 #include "CNDArray.h"
36 #include "dNDArray.h"
37 #include "fCMatrix.h"
38 #include "fMatrix.h"
39 #include "fCNDArray.h"
40 #include "fNDArray.h"
41 #include "oct-cmplx.h"
42 #include "dDiagMatrix.h"
43 #include "fDiagMatrix.h"
44 #include "CDiagMatrix.h"
45 #include "fCDiagMatrix.h"
46 #include "lo-array-errwarn.h"
47 #include "quit.h"
48
49 #include "error.h"
50 #include "xdiv.h"
51
52 static void
solve_singularity_warning(double rcond)53 solve_singularity_warning (double rcond)
54 {
55 octave::warn_singular_matrix (rcond);
56 }
57
58 template <typename T1, typename T2>
59 bool
mx_leftdiv_conform(const T1 & a,const T2 & b,blas_trans_type blas_trans)60 mx_leftdiv_conform (const T1& a, const T2& b, blas_trans_type blas_trans)
61 {
62 octave_idx_type a_nr = (blas_trans == blas_no_trans ? a.rows () : a.cols ());
63 octave_idx_type b_nr = b.rows ();
64
65 if (a_nr != b_nr)
66 {
67 octave_idx_type a_nc = (blas_trans == blas_no_trans ? a.cols ()
68 : a.rows ());
69 octave_idx_type b_nc = b.cols ();
70
71 octave::err_nonconformant (R"(operator \)", a_nr, a_nc, b_nr, b_nc);
72 }
73
74 return true;
75 }
76
77 #define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \
78 template bool mx_leftdiv_conform (const T1&, const T2&, blas_trans_type)
79
80 INSTANTIATE_MX_LEFTDIV_CONFORM (Matrix, Matrix);
81 INSTANTIATE_MX_LEFTDIV_CONFORM (Matrix, ComplexMatrix);
82 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexMatrix, Matrix);
83 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexMatrix, ComplexMatrix);
84
85 template <typename T1, typename T2>
86 bool
mx_div_conform(const T1 & a,const T2 & b)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 octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc);
98 }
99
100 return true;
101 }
102
103 #define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \
104 template bool mx_div_conform (const T1&, const T2&)
105
106 INSTANTIATE_MX_DIV_CONFORM (Matrix, Matrix);
107 INSTANTIATE_MX_DIV_CONFORM (Matrix, ComplexMatrix);
108 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, Matrix);
109 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, ComplexMatrix);
110
111 // Right division functions.
112 //
113 // op2 / op1: m cm
114 // +-- +---+----+
115 // matrix | 1 | 3 |
116 // +---+----+
117 // complex_matrix | 2 | 4 |
118 // +---+----+
119
120 // -*- 1 -*-
121 Matrix
xdiv(const Matrix & a,const Matrix & b,MatrixType & typ)122 xdiv (const Matrix& a, const Matrix& b, MatrixType& typ)
123 {
124 if (! mx_div_conform (a, b))
125 return Matrix ();
126
127 octave_idx_type info;
128 double rcond = 0.0;
129
130 Matrix result
131 = b.solve (typ, a.transpose (), info, rcond,
132 solve_singularity_warning, true, blas_trans);
133
134 return result.transpose ();
135 }
136
137 // -*- 2 -*-
138 ComplexMatrix
xdiv(const Matrix & a,const ComplexMatrix & b,MatrixType & typ)139 xdiv (const Matrix& a, const ComplexMatrix& b, MatrixType& typ)
140 {
141 if (! mx_div_conform (a, b))
142 return ComplexMatrix ();
143
144 octave_idx_type info;
145 double rcond = 0.0;
146
147 ComplexMatrix result
148 = b.solve (typ, a.transpose (), info, rcond,
149 solve_singularity_warning, true, blas_trans);
150
151 return result.transpose ();
152 }
153
154 // -*- 3 -*-
155 ComplexMatrix
xdiv(const ComplexMatrix & a,const Matrix & b,MatrixType & typ)156 xdiv (const ComplexMatrix& a, const Matrix& b, MatrixType& typ)
157 {
158 if (! mx_div_conform (a, b))
159 return ComplexMatrix ();
160
161 octave_idx_type info;
162 double rcond = 0.0;
163
164 ComplexMatrix result
165 = b.solve (typ, a.transpose (), info, rcond,
166 solve_singularity_warning, true, blas_trans);
167
168 return result.transpose ();
169 }
170
171 // -*- 4 -*-
172 ComplexMatrix
xdiv(const ComplexMatrix & a,const ComplexMatrix & b,MatrixType & typ)173 xdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType& typ)
174 {
175 if (! mx_div_conform (a, b))
176 return ComplexMatrix ();
177
178 octave_idx_type info;
179 double rcond = 0.0;
180
181 ComplexMatrix result
182 = b.solve (typ, a.transpose (), info, rcond,
183 solve_singularity_warning, true, blas_trans);
184
185 return result.transpose ();
186 }
187
188 // Funny element by element division operations.
189 //
190 // op2 \ op1: s cs
191 // +-- +---+----+
192 // matrix | 1 | 3 |
193 // +---+----+
194 // complex_matrix | 2 | 4 |
195 // +---+----+
196
197 Matrix
x_el_div(double a,const Matrix & b)198 x_el_div (double a, const Matrix& b)
199 {
200 octave_idx_type nr = b.rows ();
201 octave_idx_type nc = b.columns ();
202
203 Matrix result (nr, nc);
204
205 for (octave_idx_type j = 0; j < nc; j++)
206 for (octave_idx_type i = 0; i < nr; i++)
207 {
208 octave_quit ();
209 result (i, j) = a / b (i, j);
210 }
211
212 return result;
213 }
214
215 ComplexMatrix
x_el_div(double a,const ComplexMatrix & b)216 x_el_div (double a, const ComplexMatrix& b)
217 {
218 octave_idx_type nr = b.rows ();
219 octave_idx_type nc = b.columns ();
220
221 ComplexMatrix result (nr, nc);
222
223 for (octave_idx_type j = 0; j < nc; j++)
224 for (octave_idx_type i = 0; i < nr; i++)
225 {
226 octave_quit ();
227 result (i, j) = a / b (i, j);
228 }
229
230 return result;
231 }
232
233 ComplexMatrix
x_el_div(const Complex a,const Matrix & b)234 x_el_div (const Complex a, const Matrix& b)
235 {
236 octave_idx_type nr = b.rows ();
237 octave_idx_type nc = b.columns ();
238
239 ComplexMatrix result (nr, nc);
240
241 for (octave_idx_type j = 0; j < nc; j++)
242 for (octave_idx_type i = 0; i < nr; i++)
243 {
244 octave_quit ();
245 result (i, j) = a / b (i, j);
246 }
247
248 return result;
249 }
250
251 ComplexMatrix
x_el_div(const Complex a,const ComplexMatrix & b)252 x_el_div (const Complex a, const ComplexMatrix& b)
253 {
254 octave_idx_type nr = b.rows ();
255 octave_idx_type nc = b.columns ();
256
257 ComplexMatrix result (nr, nc);
258
259 for (octave_idx_type j = 0; j < nc; j++)
260 for (octave_idx_type i = 0; i < nr; i++)
261 {
262 octave_quit ();
263 result (i, j) = a / b (i, j);
264 }
265
266 return result;
267 }
268
269 // Funny element by element division operations.
270 //
271 // op2 \ op1: s cs
272 // +-- +---+----+
273 // N-D array | 1 | 3 |
274 // +---+----+
275 // complex N-D array | 2 | 4 |
276 // +---+----+
277
278 NDArray
x_el_div(double a,const NDArray & b)279 x_el_div (double a, const NDArray& b)
280 {
281 NDArray result (b.dims ());
282
283 for (octave_idx_type i = 0; i < b.numel (); i++)
284 {
285 octave_quit ();
286 result (i) = a / b (i);
287 }
288
289 return result;
290 }
291
292 ComplexNDArray
x_el_div(double a,const ComplexNDArray & b)293 x_el_div (double a, const ComplexNDArray& b)
294 {
295 ComplexNDArray result (b.dims ());
296
297 for (octave_idx_type i = 0; i < b.numel (); i++)
298 {
299 octave_quit ();
300 result (i) = a / b (i);
301 }
302
303 return result;
304 }
305
306 ComplexNDArray
x_el_div(const Complex a,const NDArray & b)307 x_el_div (const Complex a, const NDArray& b)
308 {
309 ComplexNDArray result (b.dims ());
310
311 for (octave_idx_type i = 0; i < b.numel (); i++)
312 {
313 octave_quit ();
314 result (i) = a / b (i);
315 }
316
317 return result;
318 }
319
320 ComplexNDArray
x_el_div(const Complex a,const ComplexNDArray & b)321 x_el_div (const Complex a, const ComplexNDArray& b)
322 {
323 ComplexNDArray result (b.dims ());
324
325 for (octave_idx_type i = 0; i < b.numel (); i++)
326 {
327 octave_quit ();
328 result (i) = a / b (i);
329 }
330
331 return result;
332 }
333
334 // Left division functions.
335 //
336 // op2 \ op1: m cm
337 // +-- +---+----+
338 // matrix | 1 | 3 |
339 // +---+----+
340 // complex_matrix | 2 | 4 |
341 // +---+----+
342
343 // -*- 1 -*-
344 Matrix
xleftdiv(const Matrix & a,const Matrix & b,MatrixType & typ,blas_trans_type transt)345 xleftdiv (const Matrix& a, const Matrix& b, MatrixType& typ,
346 blas_trans_type transt)
347 {
348 if (! mx_leftdiv_conform (a, b, transt))
349 return Matrix ();
350
351 octave_idx_type info;
352 double rcond = 0.0;
353 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
354 }
355
356 // -*- 2 -*-
357 ComplexMatrix
xleftdiv(const Matrix & a,const ComplexMatrix & b,MatrixType & typ,blas_trans_type transt)358 xleftdiv (const Matrix& a, const ComplexMatrix& b, MatrixType& typ,
359 blas_trans_type transt)
360 {
361 if (! mx_leftdiv_conform (a, b, transt))
362 return ComplexMatrix ();
363
364 octave_idx_type info;
365 double rcond = 0.0;
366
367 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
368 }
369
370 // -*- 3 -*-
371 ComplexMatrix
xleftdiv(const ComplexMatrix & a,const Matrix & b,MatrixType & typ,blas_trans_type transt)372 xleftdiv (const ComplexMatrix& a, const Matrix& b, MatrixType& typ,
373 blas_trans_type transt)
374 {
375 if (! mx_leftdiv_conform (a, b, transt))
376 return ComplexMatrix ();
377
378 octave_idx_type info;
379 double rcond = 0.0;
380 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
381 }
382
383 // -*- 4 -*-
384 ComplexMatrix
xleftdiv(const ComplexMatrix & a,const ComplexMatrix & b,MatrixType & typ,blas_trans_type transt)385 xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType& typ,
386 blas_trans_type transt)
387 {
388 if (! mx_leftdiv_conform (a, b, transt))
389 return ComplexMatrix ();
390
391 octave_idx_type info;
392 double rcond = 0.0;
393 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
394 }
395
396 static void
solve_singularity_warning(float rcond)397 solve_singularity_warning (float rcond)
398 {
399 octave::warn_singular_matrix (rcond);
400 }
401
402 INSTANTIATE_MX_LEFTDIV_CONFORM (FloatMatrix, FloatMatrix);
403 INSTANTIATE_MX_LEFTDIV_CONFORM (FloatMatrix, FloatComplexMatrix);
404 INSTANTIATE_MX_LEFTDIV_CONFORM (FloatComplexMatrix, FloatMatrix);
405 INSTANTIATE_MX_LEFTDIV_CONFORM (FloatComplexMatrix, FloatComplexMatrix);
406
407 INSTANTIATE_MX_DIV_CONFORM (FloatMatrix, FloatMatrix);
408 INSTANTIATE_MX_DIV_CONFORM (FloatMatrix, FloatComplexMatrix);
409 INSTANTIATE_MX_DIV_CONFORM (FloatComplexMatrix, FloatMatrix);
410 INSTANTIATE_MX_DIV_CONFORM (FloatComplexMatrix, FloatComplexMatrix);
411
412 // Right division functions.
413 //
414 // op2 / op1: m cm
415 // +-- +---+----+
416 // matrix | 1 | 3 |
417 // +---+----+
418 // complex_matrix | 2 | 4 |
419 // +---+----+
420
421 // -*- 1 -*-
422 FloatMatrix
xdiv(const FloatMatrix & a,const FloatMatrix & b,MatrixType & typ)423 xdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType& typ)
424 {
425 if (! mx_div_conform (a, b))
426 return FloatMatrix ();
427
428 octave_idx_type info;
429 float rcond = 0.0;
430
431 FloatMatrix result
432 = b.solve (typ, a.transpose (), info, rcond,
433 solve_singularity_warning, true, blas_trans);
434
435 return result.transpose ();
436 }
437
438 // -*- 2 -*-
439 FloatComplexMatrix
xdiv(const FloatMatrix & a,const FloatComplexMatrix & b,MatrixType & typ)440 xdiv (const FloatMatrix& a, const FloatComplexMatrix& b, MatrixType& typ)
441 {
442 if (! mx_div_conform (a, b))
443 return FloatComplexMatrix ();
444
445 octave_idx_type info;
446 float rcond = 0.0;
447
448 FloatComplexMatrix result
449 = b.solve (typ, a.transpose (), info, rcond,
450 solve_singularity_warning, true, blas_trans);
451
452 return result.transpose ();
453 }
454
455 // -*- 3 -*-
456 FloatComplexMatrix
xdiv(const FloatComplexMatrix & a,const FloatMatrix & b,MatrixType & typ)457 xdiv (const FloatComplexMatrix& a, const FloatMatrix& b, MatrixType& typ)
458 {
459 if (! mx_div_conform (a, b))
460 return FloatComplexMatrix ();
461
462 octave_idx_type info;
463 float rcond = 0.0;
464
465 FloatComplexMatrix result
466 = b.solve (typ, a.transpose (), info, rcond,
467 solve_singularity_warning, true, blas_trans);
468
469 return result.transpose ();
470 }
471
472 // -*- 4 -*-
473 FloatComplexMatrix
xdiv(const FloatComplexMatrix & a,const FloatComplexMatrix & b,MatrixType & typ)474 xdiv (const FloatComplexMatrix& a, const FloatComplexMatrix& b, MatrixType& typ)
475 {
476 if (! mx_div_conform (a, b))
477 return FloatComplexMatrix ();
478
479 octave_idx_type info;
480 float rcond = 0.0;
481
482 FloatComplexMatrix result
483 = b.solve (typ, a.transpose (), info, rcond,
484 solve_singularity_warning, true, blas_trans);
485
486 return result.transpose ();
487 }
488
489 // Funny element by element division operations.
490 //
491 // op2 \ op1: s cs
492 // +-- +---+----+
493 // matrix | 1 | 3 |
494 // +---+----+
495 // complex_matrix | 2 | 4 |
496 // +---+----+
497
498 FloatMatrix
x_el_div(float a,const FloatMatrix & b)499 x_el_div (float a, const FloatMatrix& b)
500 {
501 octave_idx_type nr = b.rows ();
502 octave_idx_type nc = b.columns ();
503
504 FloatMatrix result (nr, nc);
505
506 for (octave_idx_type j = 0; j < nc; j++)
507 for (octave_idx_type i = 0; i < nr; i++)
508 {
509 octave_quit ();
510 result (i, j) = a / b (i, j);
511 }
512
513 return result;
514 }
515
516 FloatComplexMatrix
x_el_div(float a,const FloatComplexMatrix & b)517 x_el_div (float a, const FloatComplexMatrix& b)
518 {
519 octave_idx_type nr = b.rows ();
520 octave_idx_type nc = b.columns ();
521
522 FloatComplexMatrix result (nr, nc);
523
524 for (octave_idx_type j = 0; j < nc; j++)
525 for (octave_idx_type i = 0; i < nr; i++)
526 {
527 octave_quit ();
528 result (i, j) = a / b (i, j);
529 }
530
531 return result;
532 }
533
534 FloatComplexMatrix
x_el_div(const FloatComplex a,const FloatMatrix & b)535 x_el_div (const FloatComplex a, const FloatMatrix& b)
536 {
537 octave_idx_type nr = b.rows ();
538 octave_idx_type nc = b.columns ();
539
540 FloatComplexMatrix result (nr, nc);
541
542 for (octave_idx_type j = 0; j < nc; j++)
543 for (octave_idx_type i = 0; i < nr; i++)
544 {
545 octave_quit ();
546 result (i, j) = a / b (i, j);
547 }
548
549 return result;
550 }
551
552 FloatComplexMatrix
x_el_div(const FloatComplex a,const FloatComplexMatrix & b)553 x_el_div (const FloatComplex a, const FloatComplexMatrix& b)
554 {
555 octave_idx_type nr = b.rows ();
556 octave_idx_type nc = b.columns ();
557
558 FloatComplexMatrix result (nr, nc);
559
560 for (octave_idx_type j = 0; j < nc; j++)
561 for (octave_idx_type i = 0; i < nr; i++)
562 {
563 octave_quit ();
564 result (i, j) = a / b (i, j);
565 }
566
567 return result;
568 }
569
570 // Funny element by element division operations.
571 //
572 // op2 \ op1: s cs
573 // +-- +---+----+
574 // N-D array | 1 | 3 |
575 // +---+----+
576 // complex N-D array | 2 | 4 |
577 // +---+----+
578
579 FloatNDArray
x_el_div(float a,const FloatNDArray & b)580 x_el_div (float a, const FloatNDArray& b)
581 {
582 FloatNDArray result (b.dims ());
583
584 for (octave_idx_type i = 0; i < b.numel (); i++)
585 {
586 octave_quit ();
587 result (i) = a / b (i);
588 }
589
590 return result;
591 }
592
593 FloatComplexNDArray
x_el_div(float a,const FloatComplexNDArray & b)594 x_el_div (float a, const FloatComplexNDArray& b)
595 {
596 FloatComplexNDArray result (b.dims ());
597
598 for (octave_idx_type i = 0; i < b.numel (); i++)
599 {
600 octave_quit ();
601 result (i) = a / b (i);
602 }
603
604 return result;
605 }
606
607 FloatComplexNDArray
x_el_div(const FloatComplex a,const FloatNDArray & b)608 x_el_div (const FloatComplex a, const FloatNDArray& b)
609 {
610 FloatComplexNDArray result (b.dims ());
611
612 for (octave_idx_type i = 0; i < b.numel (); i++)
613 {
614 octave_quit ();
615 result (i) = a / b (i);
616 }
617
618 return result;
619 }
620
621 FloatComplexNDArray
x_el_div(const FloatComplex a,const FloatComplexNDArray & b)622 x_el_div (const FloatComplex a, const FloatComplexNDArray& b)
623 {
624 FloatComplexNDArray result (b.dims ());
625
626 for (octave_idx_type i = 0; i < b.numel (); i++)
627 {
628 octave_quit ();
629 result (i) = a / b (i);
630 }
631
632 return result;
633 }
634
635 // Left division functions.
636 //
637 // op2 \ op1: m cm
638 // +-- +---+----+
639 // matrix | 1 | 3 |
640 // +---+----+
641 // complex_matrix | 2 | 4 |
642 // +---+----+
643
644 // -*- 1 -*-
645 FloatMatrix
xleftdiv(const FloatMatrix & a,const FloatMatrix & b,MatrixType & typ,blas_trans_type transt)646 xleftdiv (const FloatMatrix& a, const FloatMatrix& b, MatrixType& typ,
647 blas_trans_type transt)
648 {
649 if (! mx_leftdiv_conform (a, b, transt))
650 return FloatMatrix ();
651
652 octave_idx_type info;
653 float rcond = 0.0;
654 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
655 }
656
657 // -*- 2 -*-
658 FloatComplexMatrix
xleftdiv(const FloatMatrix & a,const FloatComplexMatrix & b,MatrixType & typ,blas_trans_type transt)659 xleftdiv (const FloatMatrix& a, const FloatComplexMatrix& b, MatrixType& typ,
660 blas_trans_type transt)
661 {
662 if (! mx_leftdiv_conform (a, b, transt))
663 return FloatComplexMatrix ();
664
665 octave_idx_type info;
666 float rcond = 0.0;
667
668 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
669 }
670
671 // -*- 3 -*-
672 FloatComplexMatrix
xleftdiv(const FloatComplexMatrix & a,const FloatMatrix & b,MatrixType & typ,blas_trans_type transt)673 xleftdiv (const FloatComplexMatrix& a, const FloatMatrix& b, MatrixType& typ,
674 blas_trans_type transt)
675 {
676 if (! mx_leftdiv_conform (a, b, transt))
677 return FloatComplexMatrix ();
678
679 octave_idx_type info;
680 float rcond = 0.0;
681 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
682 }
683
684 // -*- 4 -*-
685 FloatComplexMatrix
xleftdiv(const FloatComplexMatrix & a,const FloatComplexMatrix & b,MatrixType & typ,blas_trans_type transt)686 xleftdiv (const FloatComplexMatrix& a, const FloatComplexMatrix& b,
687 MatrixType& typ, blas_trans_type transt)
688 {
689 if (! mx_leftdiv_conform (a, b, transt))
690 return FloatComplexMatrix ();
691
692 octave_idx_type info;
693 float rcond = 0.0;
694 return a.solve (typ, b, info, rcond, solve_singularity_warning, true, transt);
695 }
696
697 // Diagonal matrix division.
698
699 template <typename MT, typename DMT>
700 MT
701 mdm_div_impl (const MT& a, const DMT& d)
702 {
703 if (! mx_div_conform (a, d))
704 return MT ();
705
706 octave_idx_type m = a.rows ();
707 octave_idx_type n = d.rows ();
708 octave_idx_type l = d.length ();
709 MT x (m, n);
710 typedef typename DMT::element_type S;
711 typedef typename MT::element_type T;
712 const T *aa = a.data ();
713 const S *dd = d.data ();
714 T *xx = x.fortran_vec ();
715
716 for (octave_idx_type j = 0; j < l; j++)
717 {
718 const S del = dd[j];
719 if (del != S ())
720 for (octave_idx_type i = 0; i < m; i++)
721 xx[i] = aa[i] / del;
722 else
723 for (octave_idx_type i = 0; i < m; i++)
724 xx[i] = T ();
725 aa += m; xx += m;
726 }
727
728 for (octave_idx_type i = l*m; i < n*m; i++)
729 xx[i] = T ();
730
731 return x;
732 }
733
734 // Right division functions.
735 //
736 // op2 / op1: dm cdm
737 // +-- +---+----+
738 // matrix | 1 | |
739 // +---+----+
740 // complex_matrix | 2 | 3 |
741 // +---+----+
742
743 // -*- 1 -*-
744 Matrix
xdiv(const Matrix & a,const DiagMatrix & b)745 xdiv (const Matrix& a, const DiagMatrix& b)
746 { return mdm_div_impl (a, b); }
747
748 // -*- 2 -*-
749 ComplexMatrix
xdiv(const ComplexMatrix & a,const DiagMatrix & b)750 xdiv (const ComplexMatrix& a, const DiagMatrix& b)
751 { return mdm_div_impl (a, b); }
752
753 // -*- 3 -*-
754 ComplexMatrix
xdiv(const ComplexMatrix & a,const ComplexDiagMatrix & b)755 xdiv (const ComplexMatrix& a, const ComplexDiagMatrix& b)
756 { return mdm_div_impl (a, b); }
757
758 // Right division functions, float type.
759 //
760 // op2 / op1: dm cdm
761 // +-- +---+----+
762 // matrix | 1 | |
763 // +---+----+
764 // complex_matrix | 2 | 3 |
765 // +---+----+
766
767 // -*- 1 -*-
768 FloatMatrix
xdiv(const FloatMatrix & a,const FloatDiagMatrix & b)769 xdiv (const FloatMatrix& a, const FloatDiagMatrix& b)
770 { return mdm_div_impl (a, b); }
771
772 // -*- 2 -*-
773 FloatComplexMatrix
xdiv(const FloatComplexMatrix & a,const FloatDiagMatrix & b)774 xdiv (const FloatComplexMatrix& a, const FloatDiagMatrix& b)
775 { return mdm_div_impl (a, b); }
776
777 // -*- 3 -*-
778 FloatComplexMatrix
xdiv(const FloatComplexMatrix & a,const FloatComplexDiagMatrix & b)779 xdiv (const FloatComplexMatrix& a, const FloatComplexDiagMatrix& b)
780 { return mdm_div_impl (a, b); }
781
782 template <typename MT, typename DMT>
783 MT
784 dmm_leftdiv_impl (const DMT& d, const MT& a)
785 {
786 if (! mx_leftdiv_conform (d, a, blas_no_trans))
787 return MT ();
788
789 octave_idx_type m = d.cols ();
790 octave_idx_type n = a.cols ();
791 octave_idx_type k = a.rows ();
792 octave_idx_type l = d.length ();
793 MT x (m, n);
794 typedef typename DMT::element_type S;
795 typedef typename MT::element_type T;
796 const T *aa = a.data ();
797 const S *dd = d.data ();
798 T *xx = x.fortran_vec ();
799
800 for (octave_idx_type j = 0; j < n; j++)
801 {
802 for (octave_idx_type i = 0; i < l; i++)
803 xx[i] = (dd[i] != S () ? aa[i] / dd[i] : T ());
804 for (octave_idx_type i = l; i < m; i++)
805 xx[i] = T ();
806 aa += k; xx += m;
807 }
808
809 return x;
810 }
811
812 // Left division functions.
813 //
814 // op2 \ op1: m cm
815 // +---+----+
816 // diag_matrix | 1 | 2 |
817 // +---+----+
818 // complex_diag_matrix | | 3 |
819 // +---+----+
820
821 // -*- 1 -*-
822 Matrix
xleftdiv(const DiagMatrix & a,const Matrix & b)823 xleftdiv (const DiagMatrix& a, const Matrix& b)
824 { return dmm_leftdiv_impl (a, b); }
825
826 // -*- 2 -*-
827 ComplexMatrix
xleftdiv(const DiagMatrix & a,const ComplexMatrix & b)828 xleftdiv (const DiagMatrix& a, const ComplexMatrix& b)
829 { return dmm_leftdiv_impl (a, b); }
830
831 // -*- 3 -*-
832 ComplexMatrix
xleftdiv(const ComplexDiagMatrix & a,const ComplexMatrix & b)833 xleftdiv (const ComplexDiagMatrix& a, const ComplexMatrix& b)
834 { return dmm_leftdiv_impl (a, b); }
835
836 // Left division functions, float type.
837 //
838 // op2 \ op1: m cm
839 // +---+----+
840 // diag_matrix | 1 | 2 |
841 // +---+----+
842 // complex_diag_matrix | | 3 |
843 // +---+----+
844
845 // -*- 1 -*-
846 FloatMatrix
xleftdiv(const FloatDiagMatrix & a,const FloatMatrix & b)847 xleftdiv (const FloatDiagMatrix& a, const FloatMatrix& b)
848 { return dmm_leftdiv_impl (a, b); }
849
850 // -*- 2 -*-
851 FloatComplexMatrix
xleftdiv(const FloatDiagMatrix & a,const FloatComplexMatrix & b)852 xleftdiv (const FloatDiagMatrix& a, const FloatComplexMatrix& b)
853 { return dmm_leftdiv_impl (a, b); }
854
855 // -*- 3 -*-
856 FloatComplexMatrix
xleftdiv(const FloatComplexDiagMatrix & a,const FloatComplexMatrix & b)857 xleftdiv (const FloatComplexDiagMatrix& a, const FloatComplexMatrix& b)
858 { return dmm_leftdiv_impl (a, b); }
859
860 // Diagonal by diagonal matrix division.
861
862 template <typename MT, typename DMT>
863 MT
864 dmdm_div_impl (const MT& a, const DMT& d)
865 {
866 if (! mx_div_conform (a, d))
867 return MT ();
868
869 octave_idx_type m = a.rows ();
870 octave_idx_type n = d.rows ();
871 octave_idx_type k = d.cols ();
872 octave_idx_type l = std::min (m, n);
873 octave_idx_type lk = std::min (l, k);
874 MT x (m, n);
875 typedef typename DMT::element_type S;
876 typedef typename MT::element_type T;
877 const T *aa = a.data ();
878 const S *dd = d.data ();
879 T *xx = x.fortran_vec ();
880
881 for (octave_idx_type i = 0; i < lk; i++)
882 xx[i] = (dd[i] != S () ? aa[i] / dd[i] : T ());
883 for (octave_idx_type i = lk; i < l; i++)
884 xx[i] = T ();
885
886 return x;
887 }
888
889 // Right division functions.
890 //
891 // op2 / op1: dm cdm
892 // +-- +---+----+
893 // diag_matrix | 1 | |
894 // +---+----+
895 // complex_diag_matrix | 2 | 3 |
896 // +---+----+
897
898 // -*- 1 -*-
899 DiagMatrix
xdiv(const DiagMatrix & a,const DiagMatrix & b)900 xdiv (const DiagMatrix& a, const DiagMatrix& b)
901 { return dmdm_div_impl (a, b); }
902
903 // -*- 2 -*-
904 ComplexDiagMatrix
xdiv(const ComplexDiagMatrix & a,const DiagMatrix & b)905 xdiv (const ComplexDiagMatrix& a, const DiagMatrix& b)
906 { return dmdm_div_impl (a, b); }
907
908 // -*- 3 -*-
909 ComplexDiagMatrix
xdiv(const ComplexDiagMatrix & a,const ComplexDiagMatrix & b)910 xdiv (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
911 { return dmdm_div_impl (a, b); }
912
913 // Right division functions, float type.
914 //
915 // op2 / op1: dm cdm
916 // +-- +---+----+
917 // diag_matrix | 1 | |
918 // +---+----+
919 // complex_diag_matrix | 2 | 3 |
920 // +---+----+
921
922 // -*- 1 -*-
923 FloatDiagMatrix
xdiv(const FloatDiagMatrix & a,const FloatDiagMatrix & b)924 xdiv (const FloatDiagMatrix& a, const FloatDiagMatrix& b)
925 { return dmdm_div_impl (a, b); }
926
927 // -*- 2 -*-
928 FloatComplexDiagMatrix
xdiv(const FloatComplexDiagMatrix & a,const FloatDiagMatrix & b)929 xdiv (const FloatComplexDiagMatrix& a, const FloatDiagMatrix& b)
930 { return dmdm_div_impl (a, b); }
931
932 // -*- 3 -*-
933 FloatComplexDiagMatrix
xdiv(const FloatComplexDiagMatrix & a,const FloatComplexDiagMatrix & b)934 xdiv (const FloatComplexDiagMatrix& a, const FloatComplexDiagMatrix& b)
935 { return dmdm_div_impl (a, b); }
936
937 template <typename MT, typename DMT>
938 MT
939 dmdm_leftdiv_impl (const DMT& d, const MT& a)
940 {
941 if (! mx_leftdiv_conform (d, a, blas_no_trans))
942 return MT ();
943
944 octave_idx_type m = d.cols ();
945 octave_idx_type n = a.cols ();
946 octave_idx_type k = d.rows ();
947 octave_idx_type l = std::min (m, n);
948 octave_idx_type lk = std::min (l, k);
949 MT x (m, n);
950 typedef typename DMT::element_type S;
951 typedef typename MT::element_type T;
952 const T *aa = a.data ();
953 const S *dd = d.data ();
954 T *xx = x.fortran_vec ();
955
956 for (octave_idx_type i = 0; i < lk; i++)
957 xx[i] = (dd[i] != S () ? aa[i] / dd[i] : T ());
958 for (octave_idx_type i = lk; i < l; i++)
959 xx[i] = T ();
960
961 return x;
962 }
963
964 // Left division functions.
965 //
966 // op2 \ op1: dm cdm
967 // +---+----+
968 // diag_matrix | 1 | 2 |
969 // +---+----+
970 // complex_diag_matrix | | 3 |
971 // +---+----+
972
973 // -*- 1 -*-
974 DiagMatrix
xleftdiv(const DiagMatrix & a,const DiagMatrix & b)975 xleftdiv (const DiagMatrix& a, const DiagMatrix& b)
976 { return dmdm_leftdiv_impl (a, b); }
977
978 // -*- 2 -*-
979 ComplexDiagMatrix
xleftdiv(const DiagMatrix & a,const ComplexDiagMatrix & b)980 xleftdiv (const DiagMatrix& a, const ComplexDiagMatrix& b)
981 { return dmdm_leftdiv_impl (a, b); }
982
983 // -*- 3 -*-
984 ComplexDiagMatrix
xleftdiv(const ComplexDiagMatrix & a,const ComplexDiagMatrix & b)985 xleftdiv (const ComplexDiagMatrix& a, const ComplexDiagMatrix& b)
986 { return dmdm_leftdiv_impl (a, b); }
987
988 // Left division functions, float type.
989 //
990 // op2 \ op1: dm cdm
991 // +---+----+
992 // diag_matrix | 1 | 2 |
993 // +---+----+
994 // complex_diag_matrix | | 3 |
995 // +---+----+
996
997 // -*- 1 -*-
998 FloatDiagMatrix
xleftdiv(const FloatDiagMatrix & a,const FloatDiagMatrix & b)999 xleftdiv (const FloatDiagMatrix& a, const FloatDiagMatrix& b)
1000 { return dmdm_leftdiv_impl (a, b); }
1001
1002 // -*- 2 -*-
1003 FloatComplexDiagMatrix
xleftdiv(const FloatDiagMatrix & a,const FloatComplexDiagMatrix & b)1004 xleftdiv (const FloatDiagMatrix& a, const FloatComplexDiagMatrix& b)
1005 { return dmdm_leftdiv_impl (a, b); }
1006
1007 // -*- 3 -*-
1008 FloatComplexDiagMatrix
xleftdiv(const FloatComplexDiagMatrix & a,const FloatComplexDiagMatrix & b)1009 xleftdiv (const FloatComplexDiagMatrix& a, const FloatComplexDiagMatrix& b)
1010 { return dmdm_leftdiv_impl (a, b); }
1011