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
solve_singularity_warning(double rcond)47 solve_singularity_warning (double rcond)
48 {
49   octave::warn_singular_matrix (rcond);
50 }
51 
52 template <typename T1, typename T2>
53 bool
mx_leftdiv_conform(const T1 & a,const T2 & b)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 
73 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseMatrix);
74 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseComplexMatrix);
75 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseMatrix);
76 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix);
77 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, Matrix);
78 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, ComplexMatrix);
79 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, Matrix);
80 INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, ComplexMatrix);
81 INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseMatrix);
82 INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseComplexMatrix);
83 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseMatrix);
84 INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseComplexMatrix);
85 
86 template <typename T1, typename T2>
87 bool
mx_div_conform(const T1 & a,const T2 & b)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 
107 INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseMatrix);
108 INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseComplexMatrix);
109 INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseMatrix);
110 INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix);
111 INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseMatrix);
112 INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseComplexMatrix);
113 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseMatrix);
114 INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseComplexMatrix);
115 INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, DiagMatrix);
116 INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, ComplexDiagMatrix);
117 INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, DiagMatrix);
118 INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, ComplexDiagMatrix);
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
xdiv(const Matrix & a,const SparseMatrix & b,MatrixType & typ)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,
147                               solve_singularity_warning);
148 
149   typ = btyp.transpose ();
150   return result.transpose ();
151 }
152 
153 // -*- 2 -*-
154 ComplexMatrix
xdiv(const Matrix & a,const SparseComplexMatrix & b,MatrixType & typ)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 -*-
174 ComplexMatrix
xdiv(const ComplexMatrix & a,const SparseMatrix & b,MatrixType & typ)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 -*-
194 ComplexMatrix
xdiv(const ComplexMatrix & a,const SparseComplexMatrix & b,MatrixType & typ)195 xdiv (const ComplexMatrix& a, const SparseComplexMatrix& b, MatrixType& typ)
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 -*-
214 SparseMatrix
xdiv(const SparseMatrix & a,const SparseMatrix & b,MatrixType & typ)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,
227                                     solve_singularity_warning);
228 
229   typ = btyp.transpose ();
230   return result.transpose ();
231 }
232 
233 // -*- 6 -*-
234 SparseComplexMatrix
xdiv(const SparseMatrix & a,const SparseComplexMatrix & b,MatrixType & typ)235 xdiv (const SparseMatrix& a, const SparseComplexMatrix& b, MatrixType& typ)
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 -*-
254 SparseComplexMatrix
xdiv(const SparseComplexMatrix & a,const SparseMatrix & b,MatrixType & typ)255 xdiv (const SparseComplexMatrix& a, const SparseMatrix& b, MatrixType& typ)
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 -*-
274 SparseComplexMatrix
xdiv(const SparseComplexMatrix & a,const SparseComplexMatrix & b,MatrixType & typ)275 xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
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 -*-
336 SparseMatrix
xdiv(const SparseMatrix & a,const DiagMatrix & b,MatrixType &)337 xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType &)
338 {
339   return do_rightdiv_sm_dm<SparseMatrix> (a, b);
340 }
341 
342 // -*- 10 -*-
343 SparseComplexMatrix
xdiv(const SparseMatrix & a,const ComplexDiagMatrix & b,MatrixType &)344 xdiv (const SparseMatrix& a, const ComplexDiagMatrix& b, MatrixType &)
345 {
346   return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
347 }
348 
349 // -*- 11 -*-
350 SparseComplexMatrix
xdiv(const SparseComplexMatrix & a,const DiagMatrix & b,MatrixType &)351 xdiv (const SparseComplexMatrix& a, const DiagMatrix& b, MatrixType &)
352 {
353   return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b);
354 }
355 
356 // -*- 12 -*-
357 SparseComplexMatrix
xdiv(const SparseComplexMatrix & a,const ComplexDiagMatrix & b,MatrixType &)358 xdiv (const SparseComplexMatrix& a, const ComplexDiagMatrix& b, MatrixType &)
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
x_el_div(double a,const SparseMatrix & b)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 
396 ComplexMatrix
x_el_div(double a,const SparseComplexMatrix & b)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 
402   ComplexMatrix result (nr, nc, Complex (octave::numeric_limits<double>::NaN (),
403                                          octave::numeric_limits<double>::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 
415 ComplexMatrix
x_el_div(const Complex a,const SparseMatrix & b)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 
433 ComplexMatrix
x_el_div(const Complex a,const SparseComplexMatrix & b)434 x_el_div (const Complex a, const SparseComplexMatrix& b)
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
xleftdiv(const SparseMatrix & a,const Matrix & b,MatrixType & typ)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 -*-
477 ComplexMatrix
xleftdiv(const SparseMatrix & a,const ComplexMatrix & b,MatrixType & typ)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 -*-
489 SparseMatrix
xleftdiv(const SparseMatrix & a,const SparseMatrix & b,MatrixType & typ)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 -*-
501 SparseComplexMatrix
xleftdiv(const SparseMatrix & a,const SparseComplexMatrix & b,MatrixType & typ)502 xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b, MatrixType& typ)
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 -*-
513 ComplexMatrix
xleftdiv(const SparseComplexMatrix & a,const Matrix & b,MatrixType & typ)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 -*-
525 ComplexMatrix
xleftdiv(const SparseComplexMatrix & a,const ComplexMatrix & b,MatrixType & typ)526 xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b, MatrixType& typ)
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 -*-
537 SparseComplexMatrix
xleftdiv(const SparseComplexMatrix & a,const SparseMatrix & b,MatrixType & typ)538 xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b, MatrixType& typ)
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 -*-
549 SparseComplexMatrix
xleftdiv(const SparseComplexMatrix & a,const SparseComplexMatrix & b,MatrixType & typ)550 xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
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 -*-
609 SparseMatrix
xleftdiv(const DiagMatrix & d,const SparseMatrix & a,MatrixType &)610 xleftdiv (const DiagMatrix& d, const SparseMatrix& a,  MatrixType&)
611 {
612   return do_leftdiv_dm_sm<SparseMatrix> (d, a);
613 }
614 
615 // -*- 10 -*-
616 SparseComplexMatrix
xleftdiv(const DiagMatrix & d,const SparseComplexMatrix & a,MatrixType &)617 xleftdiv (const DiagMatrix& d, const SparseComplexMatrix& a,  MatrixType&)
618 {
619   return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
620 }
621 
622 // -*- 11 -*-
623 SparseComplexMatrix
xleftdiv(const ComplexDiagMatrix & d,const SparseMatrix & a,MatrixType &)624 xleftdiv (const ComplexDiagMatrix& d, const SparseMatrix& a,  MatrixType&)
625 {
626   return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
627 }
628 
629 // -*- 12 -*-
630 SparseComplexMatrix
xleftdiv(const ComplexDiagMatrix & d,const SparseComplexMatrix & a,MatrixType &)631 xleftdiv (const ComplexDiagMatrix& d, const SparseComplexMatrix& a,
632           MatrixType&)
633 {
634   return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a);
635 }
636