1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-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 <algorithm>
31 #include <complex>
32 #include <istream>
33 #include <limits>
34 #include <ostream>
35
36 #include "Array-util.h"
37 #include "CDiagMatrix.h"
38 #include "CMatrix.h"
39 #include "CNDArray.h"
40 #include "CRowVector.h"
41 #include "DET.h"
42 #include "boolMatrix.h"
43 #include "chMatrix.h"
44 #include "chol.h"
45 #include "dDiagMatrix.h"
46 #include "dMatrix.h"
47 #include "dRowVector.h"
48 #include "lo-blas-proto.h"
49 #include "lo-error.h"
50 #include "lo-ieee.h"
51 #include "lo-lapack-proto.h"
52 #include "lo-mappers.h"
53 #include "lo-utils.h"
54 #include "mx-cm-dm.h"
55 #include "mx-cm-s.h"
56 #include "mx-dm-cm.h"
57 #include "mx-inlines.cc"
58 #include "mx-op-defs.h"
59 #include "oct-cmplx.h"
60 #include "oct-fftw.h"
61 #include "oct-locbuf.h"
62 #include "oct-norm.h"
63 #include "schur.h"
64 #include "svd.h"
65
66 static const Complex Complex_NaN_result (octave::numeric_limits<double>::NaN (),
67 octave::numeric_limits<double>::NaN ());
68
69 // Complex Matrix class
70
ComplexMatrix(const Matrix & a)71 ComplexMatrix::ComplexMatrix (const Matrix& a)
72 : ComplexNDArray (a)
73 { }
74
ComplexMatrix(const RowVector & rv)75 ComplexMatrix::ComplexMatrix (const RowVector& rv)
76 : ComplexNDArray (rv)
77 { }
78
ComplexMatrix(const ColumnVector & cv)79 ComplexMatrix::ComplexMatrix (const ColumnVector& cv)
80 : ComplexNDArray (cv)
81 { }
82
ComplexMatrix(const DiagMatrix & a)83 ComplexMatrix::ComplexMatrix (const DiagMatrix& a)
84 : ComplexNDArray (a.dims (), 0.0)
85 {
86 for (octave_idx_type i = 0; i < a.length (); i++)
87 elem (i, i) = a.elem (i, i);
88 }
89
ComplexMatrix(const MDiagArray2<double> & a)90 ComplexMatrix::ComplexMatrix (const MDiagArray2<double>& a)
91 : ComplexNDArray (a.dims (), 0.0)
92 {
93 for (octave_idx_type i = 0; i < a.length (); i++)
94 elem (i, i) = a.elem (i, i);
95 }
96
ComplexMatrix(const DiagArray2<double> & a)97 ComplexMatrix::ComplexMatrix (const DiagArray2<double>& a)
98 : ComplexNDArray (a.dims (), 0.0)
99 {
100 for (octave_idx_type i = 0; i < a.length (); i++)
101 elem (i, i) = a.elem (i, i);
102 }
103
ComplexMatrix(const ComplexRowVector & rv)104 ComplexMatrix::ComplexMatrix (const ComplexRowVector& rv)
105 : ComplexNDArray (rv)
106 { }
107
ComplexMatrix(const ComplexColumnVector & cv)108 ComplexMatrix::ComplexMatrix (const ComplexColumnVector& cv)
109 : ComplexNDArray (cv)
110 { }
111
ComplexMatrix(const ComplexDiagMatrix & a)112 ComplexMatrix::ComplexMatrix (const ComplexDiagMatrix& a)
113 : ComplexNDArray (a.dims (), 0.0)
114 {
115 for (octave_idx_type i = 0; i < a.length (); i++)
116 elem (i, i) = a.elem (i, i);
117 }
118
ComplexMatrix(const MDiagArray2<Complex> & a)119 ComplexMatrix::ComplexMatrix (const MDiagArray2<Complex>& a)
120 : ComplexNDArray (a.dims (), 0.0)
121 {
122 for (octave_idx_type i = 0; i < a.length (); i++)
123 elem (i, i) = a.elem (i, i);
124 }
125
ComplexMatrix(const DiagArray2<Complex> & a)126 ComplexMatrix::ComplexMatrix (const DiagArray2<Complex>& a)
127 : ComplexNDArray (a.dims (), 0.0)
128 {
129 for (octave_idx_type i = 0; i < a.length (); i++)
130 elem (i, i) = a.elem (i, i);
131 }
132
133 // FIXME: could we use a templated mixed-type copy function here?
134
ComplexMatrix(const boolMatrix & a)135 ComplexMatrix::ComplexMatrix (const boolMatrix& a)
136 : ComplexNDArray (a)
137 { }
138
ComplexMatrix(const charMatrix & a)139 ComplexMatrix::ComplexMatrix (const charMatrix& a)
140 : ComplexNDArray (a.dims (), 0.0)
141 {
142 for (octave_idx_type i = 0; i < a.rows (); i++)
143 for (octave_idx_type j = 0; j < a.cols (); j++)
144 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
145 }
146
ComplexMatrix(const Matrix & re,const Matrix & im)147 ComplexMatrix::ComplexMatrix (const Matrix& re, const Matrix& im)
148 : ComplexNDArray (re.dims ())
149 {
150 if (im.rows () != rows () || im.cols () != cols ())
151 (*current_liboctave_error_handler) ("complex: internal error");
152
153 octave_idx_type nel = numel ();
154 for (octave_idx_type i = 0; i < nel; i++)
155 xelem (i) = Complex (re(i), im(i));
156 }
157
158 bool
operator ==(const ComplexMatrix & a) const159 ComplexMatrix::operator == (const ComplexMatrix& a) const
160 {
161 if (rows () != a.rows () || cols () != a.cols ())
162 return false;
163
164 return mx_inline_equal (numel (), data (), a.data ());
165 }
166
167 bool
operator !=(const ComplexMatrix & a) const168 ComplexMatrix::operator != (const ComplexMatrix& a) const
169 {
170 return !(*this == a);
171 }
172
173 bool
ishermitian(void) const174 ComplexMatrix::ishermitian (void) const
175 {
176 octave_idx_type nr = rows ();
177 octave_idx_type nc = cols ();
178
179 if (issquare () && nr > 0)
180 {
181 for (octave_idx_type i = 0; i < nr; i++)
182 for (octave_idx_type j = i; j < nc; j++)
183 if (elem (i, j) != conj (elem (j, i)))
184 return false;
185
186 return true;
187 }
188
189 return false;
190 }
191
192 // destructive insert/delete/reorder operations
193
194 ComplexMatrix&
insert(const Matrix & a,octave_idx_type r,octave_idx_type c)195 ComplexMatrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c)
196 {
197 octave_idx_type a_nr = a.rows ();
198 octave_idx_type a_nc = a.cols ();
199
200 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
201 (*current_liboctave_error_handler) ("range error for insert");
202
203 if (a_nr >0 && a_nc > 0)
204 {
205 make_unique ();
206
207 for (octave_idx_type j = 0; j < a_nc; j++)
208 for (octave_idx_type i = 0; i < a_nr; i++)
209 xelem (r+i, c+j) = a.elem (i, j);
210 }
211
212 return *this;
213 }
214
215 ComplexMatrix&
insert(const RowVector & a,octave_idx_type r,octave_idx_type c)216 ComplexMatrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c)
217 {
218 octave_idx_type a_len = a.numel ();
219
220 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
221 (*current_liboctave_error_handler) ("range error for insert");
222
223 if (a_len > 0)
224 {
225 make_unique ();
226
227 for (octave_idx_type i = 0; i < a_len; i++)
228 xelem (r, c+i) = a.elem (i);
229 }
230
231 return *this;
232 }
233
234 ComplexMatrix&
insert(const ColumnVector & a,octave_idx_type r,octave_idx_type c)235 ComplexMatrix::insert (const ColumnVector& a,
236 octave_idx_type r, octave_idx_type c)
237 {
238 octave_idx_type a_len = a.numel ();
239
240 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
241 (*current_liboctave_error_handler) ("range error for insert");
242
243 if (a_len > 0)
244 {
245 make_unique ();
246
247 for (octave_idx_type i = 0; i < a_len; i++)
248 xelem (r+i, c) = a.elem (i);
249 }
250
251 return *this;
252 }
253
254 ComplexMatrix&
insert(const DiagMatrix & a,octave_idx_type r,octave_idx_type c)255 ComplexMatrix::insert (const DiagMatrix& a,
256 octave_idx_type r, octave_idx_type c)
257 {
258 octave_idx_type a_nr = a.rows ();
259 octave_idx_type a_nc = a.cols ();
260
261 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
262 (*current_liboctave_error_handler) ("range error for insert");
263
264 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
265
266 octave_idx_type a_len = a.length ();
267
268 if (a_len > 0)
269 {
270 make_unique ();
271
272 for (octave_idx_type i = 0; i < a_len; i++)
273 xelem (r+i, c+i) = a.elem (i, i);
274 }
275
276 return *this;
277 }
278
279 ComplexMatrix&
insert(const ComplexMatrix & a,octave_idx_type r,octave_idx_type c)280 ComplexMatrix::insert (const ComplexMatrix& a,
281 octave_idx_type r, octave_idx_type c)
282 {
283 ComplexNDArray::insert (a, r, c);
284 return *this;
285 }
286
287 ComplexMatrix&
insert(const ComplexRowVector & a,octave_idx_type r,octave_idx_type c)288 ComplexMatrix::insert (const ComplexRowVector& a,
289 octave_idx_type r, octave_idx_type c)
290 {
291 octave_idx_type a_len = a.numel ();
292 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
293 (*current_liboctave_error_handler) ("range error for insert");
294
295 for (octave_idx_type i = 0; i < a_len; i++)
296 elem (r, c+i) = a.elem (i);
297
298 return *this;
299 }
300
301 ComplexMatrix&
insert(const ComplexColumnVector & a,octave_idx_type r,octave_idx_type c)302 ComplexMatrix::insert (const ComplexColumnVector& a,
303 octave_idx_type r, octave_idx_type c)
304 {
305 octave_idx_type a_len = a.numel ();
306
307 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
308 (*current_liboctave_error_handler) ("range error for insert");
309
310 if (a_len > 0)
311 {
312 make_unique ();
313
314 for (octave_idx_type i = 0; i < a_len; i++)
315 xelem (r+i, c) = a.elem (i);
316 }
317
318 return *this;
319 }
320
321 ComplexMatrix&
insert(const ComplexDiagMatrix & a,octave_idx_type r,octave_idx_type c)322 ComplexMatrix::insert (const ComplexDiagMatrix& a,
323 octave_idx_type r, octave_idx_type c)
324 {
325 octave_idx_type a_nr = a.rows ();
326 octave_idx_type a_nc = a.cols ();
327
328 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
329 (*current_liboctave_error_handler) ("range error for insert");
330
331 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
332
333 octave_idx_type a_len = a.length ();
334
335 if (a_len > 0)
336 {
337 make_unique ();
338
339 for (octave_idx_type i = 0; i < a_len; i++)
340 xelem (r+i, c+i) = a.elem (i, i);
341 }
342
343 return *this;
344 }
345
346 ComplexMatrix&
fill(double val)347 ComplexMatrix::fill (double val)
348 {
349 octave_idx_type nr = rows ();
350 octave_idx_type nc = cols ();
351
352 if (nr > 0 && nc > 0)
353 {
354 make_unique ();
355
356 for (octave_idx_type j = 0; j < nc; j++)
357 for (octave_idx_type i = 0; i < nr; i++)
358 xelem (i, j) = val;
359 }
360
361 return *this;
362 }
363
364 ComplexMatrix&
fill(const Complex & val)365 ComplexMatrix::fill (const Complex& val)
366 {
367 octave_idx_type nr = rows ();
368 octave_idx_type nc = cols ();
369
370 if (nr > 0 && nc > 0)
371 {
372 make_unique ();
373
374 for (octave_idx_type j = 0; j < nc; j++)
375 for (octave_idx_type i = 0; i < nr; i++)
376 xelem (i, j) = val;
377 }
378
379 return *this;
380 }
381
382 ComplexMatrix&
fill(double val,octave_idx_type r1,octave_idx_type c1,octave_idx_type r2,octave_idx_type c2)383 ComplexMatrix::fill (double val, octave_idx_type r1, octave_idx_type c1,
384 octave_idx_type r2, octave_idx_type c2)
385 {
386 octave_idx_type nr = rows ();
387 octave_idx_type nc = cols ();
388
389 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
390 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
391 (*current_liboctave_error_handler) ("range error for fill");
392
393 if (r1 > r2) { std::swap (r1, r2); }
394 if (c1 > c2) { std::swap (c1, c2); }
395
396 if (r2 >= r1 && c2 >= c1)
397 {
398 make_unique ();
399
400 for (octave_idx_type j = c1; j <= c2; j++)
401 for (octave_idx_type i = r1; i <= r2; i++)
402 xelem (i, j) = val;
403 }
404
405 return *this;
406 }
407
408 ComplexMatrix&
fill(const Complex & val,octave_idx_type r1,octave_idx_type c1,octave_idx_type r2,octave_idx_type c2)409 ComplexMatrix::fill (const Complex& val, octave_idx_type r1, octave_idx_type c1,
410 octave_idx_type r2, octave_idx_type c2)
411 {
412 octave_idx_type nr = rows ();
413 octave_idx_type nc = cols ();
414
415 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
416 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
417 (*current_liboctave_error_handler) ("range error for fill");
418
419 if (r1 > r2) { std::swap (r1, r2); }
420 if (c1 > c2) { std::swap (c1, c2); }
421
422 if (r2 >= r1 && c2 >=c1)
423 {
424 make_unique ();
425
426 for (octave_idx_type j = c1; j <= c2; j++)
427 for (octave_idx_type i = r1; i <= r2; i++)
428 xelem (i, j) = val;
429 }
430
431 return *this;
432 }
433
434 ComplexMatrix
append(const Matrix & a) const435 ComplexMatrix::append (const Matrix& a) const
436 {
437 octave_idx_type nr = rows ();
438 octave_idx_type nc = cols ();
439 if (nr != a.rows ())
440 (*current_liboctave_error_handler) ("row dimension mismatch for append");
441
442 octave_idx_type nc_insert = nc;
443 ComplexMatrix retval (nr, nc + a.cols ());
444 retval.insert (*this, 0, 0);
445 retval.insert (a, 0, nc_insert);
446 return retval;
447 }
448
449 ComplexMatrix
append(const RowVector & a) const450 ComplexMatrix::append (const RowVector& a) const
451 {
452 octave_idx_type nr = rows ();
453 octave_idx_type nc = cols ();
454 if (nr != 1)
455 (*current_liboctave_error_handler) ("row dimension mismatch for append");
456
457 octave_idx_type nc_insert = nc;
458 ComplexMatrix retval (nr, nc + a.numel ());
459 retval.insert (*this, 0, 0);
460 retval.insert (a, 0, nc_insert);
461 return retval;
462 }
463
464 ComplexMatrix
append(const ColumnVector & a) const465 ComplexMatrix::append (const ColumnVector& a) const
466 {
467 octave_idx_type nr = rows ();
468 octave_idx_type nc = cols ();
469 if (nr != a.numel ())
470 (*current_liboctave_error_handler) ("row dimension mismatch for append");
471
472 octave_idx_type nc_insert = nc;
473 ComplexMatrix retval (nr, nc + 1);
474 retval.insert (*this, 0, 0);
475 retval.insert (a, 0, nc_insert);
476 return retval;
477 }
478
479 ComplexMatrix
append(const DiagMatrix & a) const480 ComplexMatrix::append (const DiagMatrix& a) const
481 {
482 octave_idx_type nr = rows ();
483 octave_idx_type nc = cols ();
484 if (nr != a.rows ())
485 (*current_liboctave_error_handler) ("row dimension mismatch for append");
486
487 octave_idx_type nc_insert = nc;
488 ComplexMatrix retval (nr, nc + a.cols ());
489 retval.insert (*this, 0, 0);
490 retval.insert (a, 0, nc_insert);
491 return retval;
492 }
493
494 ComplexMatrix
append(const ComplexMatrix & a) const495 ComplexMatrix::append (const ComplexMatrix& a) const
496 {
497 octave_idx_type nr = rows ();
498 octave_idx_type nc = cols ();
499 if (nr != a.rows ())
500 (*current_liboctave_error_handler) ("row dimension mismatch for append");
501
502 octave_idx_type nc_insert = nc;
503 ComplexMatrix retval (nr, nc + a.cols ());
504 retval.insert (*this, 0, 0);
505 retval.insert (a, 0, nc_insert);
506 return retval;
507 }
508
509 ComplexMatrix
append(const ComplexRowVector & a) const510 ComplexMatrix::append (const ComplexRowVector& a) const
511 {
512 octave_idx_type nr = rows ();
513 octave_idx_type nc = cols ();
514 if (nr != 1)
515 (*current_liboctave_error_handler) ("row dimension mismatch for append");
516
517 octave_idx_type nc_insert = nc;
518 ComplexMatrix retval (nr, nc + a.numel ());
519 retval.insert (*this, 0, 0);
520 retval.insert (a, 0, nc_insert);
521 return retval;
522 }
523
524 ComplexMatrix
append(const ComplexColumnVector & a) const525 ComplexMatrix::append (const ComplexColumnVector& a) const
526 {
527 octave_idx_type nr = rows ();
528 octave_idx_type nc = cols ();
529 if (nr != a.numel ())
530 (*current_liboctave_error_handler) ("row dimension mismatch for append");
531
532 octave_idx_type nc_insert = nc;
533 ComplexMatrix retval (nr, nc + 1);
534 retval.insert (*this, 0, 0);
535 retval.insert (a, 0, nc_insert);
536 return retval;
537 }
538
539 ComplexMatrix
append(const ComplexDiagMatrix & a) const540 ComplexMatrix::append (const ComplexDiagMatrix& a) const
541 {
542 octave_idx_type nr = rows ();
543 octave_idx_type nc = cols ();
544 if (nr != a.rows ())
545 (*current_liboctave_error_handler) ("row dimension mismatch for append");
546
547 octave_idx_type nc_insert = nc;
548 ComplexMatrix retval (nr, nc + a.cols ());
549 retval.insert (*this, 0, 0);
550 retval.insert (a, 0, nc_insert);
551 return retval;
552 }
553
554 ComplexMatrix
stack(const Matrix & a) const555 ComplexMatrix::stack (const Matrix& a) const
556 {
557 octave_idx_type nr = rows ();
558 octave_idx_type nc = cols ();
559 if (nc != a.cols ())
560 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
561
562 octave_idx_type nr_insert = nr;
563 ComplexMatrix retval (nr + a.rows (), nc);
564 retval.insert (*this, 0, 0);
565 retval.insert (a, nr_insert, 0);
566 return retval;
567 }
568
569 ComplexMatrix
stack(const RowVector & a) const570 ComplexMatrix::stack (const RowVector& a) const
571 {
572 octave_idx_type nr = rows ();
573 octave_idx_type nc = cols ();
574 if (nc != a.numel ())
575 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
576
577 octave_idx_type nr_insert = nr;
578 ComplexMatrix retval (nr + 1, nc);
579 retval.insert (*this, 0, 0);
580 retval.insert (a, nr_insert, 0);
581 return retval;
582 }
583
584 ComplexMatrix
stack(const ColumnVector & a) const585 ComplexMatrix::stack (const ColumnVector& a) const
586 {
587 octave_idx_type nr = rows ();
588 octave_idx_type nc = cols ();
589 if (nc != 1)
590 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
591
592 octave_idx_type nr_insert = nr;
593 ComplexMatrix retval (nr + a.numel (), nc);
594 retval.insert (*this, 0, 0);
595 retval.insert (a, nr_insert, 0);
596 return retval;
597 }
598
599 ComplexMatrix
stack(const DiagMatrix & a) const600 ComplexMatrix::stack (const DiagMatrix& a) const
601 {
602 octave_idx_type nr = rows ();
603 octave_idx_type nc = cols ();
604 if (nc != a.cols ())
605 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
606
607 octave_idx_type nr_insert = nr;
608 ComplexMatrix retval (nr + a.rows (), nc);
609 retval.insert (*this, 0, 0);
610 retval.insert (a, nr_insert, 0);
611 return retval;
612 }
613
614 ComplexMatrix
stack(const ComplexMatrix & a) const615 ComplexMatrix::stack (const ComplexMatrix& a) const
616 {
617 octave_idx_type nr = rows ();
618 octave_idx_type nc = cols ();
619 if (nc != a.cols ())
620 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
621
622 octave_idx_type nr_insert = nr;
623 ComplexMatrix retval (nr + a.rows (), nc);
624 retval.insert (*this, 0, 0);
625 retval.insert (a, nr_insert, 0);
626 return retval;
627 }
628
629 ComplexMatrix
stack(const ComplexRowVector & a) const630 ComplexMatrix::stack (const ComplexRowVector& a) const
631 {
632 octave_idx_type nr = rows ();
633 octave_idx_type nc = cols ();
634 if (nc != a.numel ())
635 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
636
637 octave_idx_type nr_insert = nr;
638 ComplexMatrix retval (nr + 1, nc);
639 retval.insert (*this, 0, 0);
640 retval.insert (a, nr_insert, 0);
641 return retval;
642 }
643
644 ComplexMatrix
stack(const ComplexColumnVector & a) const645 ComplexMatrix::stack (const ComplexColumnVector& a) const
646 {
647 octave_idx_type nr = rows ();
648 octave_idx_type nc = cols ();
649 if (nc != 1)
650 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
651
652 octave_idx_type nr_insert = nr;
653 ComplexMatrix retval (nr + a.numel (), nc);
654 retval.insert (*this, 0, 0);
655 retval.insert (a, nr_insert, 0);
656 return retval;
657 }
658
659 ComplexMatrix
stack(const ComplexDiagMatrix & a) const660 ComplexMatrix::stack (const ComplexDiagMatrix& a) const
661 {
662 octave_idx_type nr = rows ();
663 octave_idx_type nc = cols ();
664 if (nc != a.cols ())
665 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
666
667 octave_idx_type nr_insert = nr;
668 ComplexMatrix retval (nr + a.rows (), nc);
669 retval.insert (*this, 0, 0);
670 retval.insert (a, nr_insert, 0);
671 return retval;
672 }
673
674 ComplexMatrix
conj(const ComplexMatrix & a)675 conj (const ComplexMatrix& a)
676 {
677 return do_mx_unary_map<Complex, Complex, std::conj<double>> (a);
678 }
679
680 // resize is the destructive equivalent for this one
681
682 ComplexMatrix
extract(octave_idx_type r1,octave_idx_type c1,octave_idx_type r2,octave_idx_type c2) const683 ComplexMatrix::extract (octave_idx_type r1, octave_idx_type c1,
684 octave_idx_type r2, octave_idx_type c2) const
685 {
686 if (r1 > r2) { std::swap (r1, r2); }
687 if (c1 > c2) { std::swap (c1, c2); }
688
689 return index (idx_vector (r1, r2+1), idx_vector (c1, c2+1));
690 }
691
692 ComplexMatrix
extract_n(octave_idx_type r1,octave_idx_type c1,octave_idx_type nr,octave_idx_type nc) const693 ComplexMatrix::extract_n (octave_idx_type r1, octave_idx_type c1,
694 octave_idx_type nr, octave_idx_type nc) const
695 {
696 return index (idx_vector (r1, r1 + nr), idx_vector (c1, c1 + nc));
697 }
698
699 // extract row or column i.
700
701 ComplexRowVector
row(octave_idx_type i) const702 ComplexMatrix::row (octave_idx_type i) const
703 {
704 return index (idx_vector (i), idx_vector::colon);
705 }
706
707 ComplexColumnVector
column(octave_idx_type i) const708 ComplexMatrix::column (octave_idx_type i) const
709 {
710 return index (idx_vector::colon, idx_vector (i));
711 }
712
713 // Local function to calculate the 1-norm.
714 static
715 double
norm1(const ComplexMatrix & a)716 norm1 (const ComplexMatrix& a)
717 {
718 double anorm = 0.0;
719 RowVector colsum = a.abs ().sum ().row (0);
720
721 for (octave_idx_type i = 0; i < colsum.numel (); i++)
722 {
723 double sum = colsum.xelem (i);
724 if (octave::math::isinf (sum) || octave::math::isnan (sum))
725 {
726 anorm = sum; // Pass Inf or NaN to output
727 break;
728 }
729 else
730 anorm = std::max (anorm, sum);
731 }
732
733 return anorm;
734 }
735
736 ComplexMatrix
inverse(void) const737 ComplexMatrix::inverse (void) const
738 {
739 octave_idx_type info;
740 double rcon;
741 MatrixType mattype (*this);
742 return inverse (mattype, info, rcon, 0, 0);
743 }
744
745 ComplexMatrix
inverse(octave_idx_type & info) const746 ComplexMatrix::inverse (octave_idx_type& info) const
747 {
748 double rcon;
749 MatrixType mattype (*this);
750 return inverse (mattype, info, rcon, 0, 0);
751 }
752
753 ComplexMatrix
inverse(octave_idx_type & info,double & rcon,bool force,bool calc_cond) const754 ComplexMatrix::inverse (octave_idx_type& info, double& rcon, bool force,
755 bool calc_cond) const
756 {
757 MatrixType mattype (*this);
758 return inverse (mattype, info, rcon, force, calc_cond);
759 }
760
761 ComplexMatrix
inverse(MatrixType & mattype) const762 ComplexMatrix::inverse (MatrixType& mattype) const
763 {
764 octave_idx_type info;
765 double rcon;
766 return inverse (mattype, info, rcon, 0, 0);
767 }
768
769 ComplexMatrix
inverse(MatrixType & mattype,octave_idx_type & info) const770 ComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const
771 {
772 double rcon;
773 return inverse (mattype, info, rcon, 0, 0);
774 }
775
776 ComplexMatrix
tinverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const777 ComplexMatrix::tinverse (MatrixType& mattype, octave_idx_type& info,
778 double& rcon, bool force, bool calc_cond) const
779 {
780 ComplexMatrix retval;
781
782 F77_INT nr = octave::to_f77_int (rows ());
783 F77_INT nc = octave::to_f77_int (cols ());
784
785 if (nr != nc || nr == 0 || nc == 0)
786 (*current_liboctave_error_handler) ("inverse requires square matrix");
787
788 int typ = mattype.type ();
789 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
790 char udiag = 'N';
791 retval = *this;
792 Complex *tmp_data = retval.fortran_vec ();
793
794 F77_INT tmp_info = 0;
795
796 F77_XFCN (ztrtri, ZTRTRI,(F77_CONST_CHAR_ARG2 (&uplo, 1),
797 F77_CONST_CHAR_ARG2 (&udiag, 1),
798 nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, tmp_info
799 F77_CHAR_ARG_LEN (1)
800 F77_CHAR_ARG_LEN (1)));
801
802 info = tmp_info;
803
804 // Throw away extra info LAPACK gives so as to not change output.
805 rcon = 0.0;
806 if (info != 0)
807 info = -1;
808 else if (calc_cond)
809 {
810 F77_INT ztrcon_info = 0;
811 char job = '1';
812
813 OCTAVE_LOCAL_BUFFER (Complex, cwork, 2*nr);
814 OCTAVE_LOCAL_BUFFER (double, rwork, nr);
815
816 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
817 F77_CONST_CHAR_ARG2 (&uplo, 1),
818 F77_CONST_CHAR_ARG2 (&udiag, 1),
819 nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, rcon,
820 F77_DBLE_CMPLX_ARG (cwork), rwork, ztrcon_info
821 F77_CHAR_ARG_LEN (1)
822 F77_CHAR_ARG_LEN (1)
823 F77_CHAR_ARG_LEN (1)));
824
825 if (ztrcon_info != 0)
826 info = -1;
827 }
828
829 if (info == -1 && ! force)
830 retval = *this; // Restore matrix contents.
831
832 return retval;
833 }
834
835 ComplexMatrix
finverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const836 ComplexMatrix::finverse (MatrixType& mattype, octave_idx_type& info,
837 double& rcon, bool force, bool calc_cond) const
838 {
839 ComplexMatrix retval;
840
841 F77_INT nr = octave::to_f77_int (rows ());
842 F77_INT nc = octave::to_f77_int (cols ());
843
844 if (nr != nc)
845 (*current_liboctave_error_handler) ("inverse requires square matrix");
846
847 Array<F77_INT> ipvt (dim_vector (nr, 1));
848 F77_INT *pipvt = ipvt.fortran_vec ();
849
850 retval = *this;
851 Complex *tmp_data = retval.fortran_vec ();
852
853 Array<Complex> z (dim_vector (1, 1));
854 F77_INT lwork = -1;
855
856 // Query the optimum work array size.
857
858 F77_INT tmp_info = 0;
859
860 F77_XFCN (zgetri, ZGETRI, (nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, pipvt,
861 F77_DBLE_CMPLX_ARG (z.fortran_vec ()), lwork,
862 tmp_info));
863
864 lwork = static_cast<F77_INT> (std::real (z(0)));
865 lwork = (lwork < 2 * nc ? 2 * nc : lwork);
866 z.resize (dim_vector (lwork, 1));
867 Complex *pz = z.fortran_vec ();
868
869 info = 0;
870 tmp_info = 0;
871
872 // Calculate norm of the matrix (always, see bug #45577) for later use.
873 double anorm = norm1 (retval);
874
875 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
876 // and bug #46330, segfault with matrices containing Inf & NaN
877 if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
878 info = -1;
879 else
880 {
881 F77_XFCN (zgetrf, ZGETRF, (nc, nc, F77_DBLE_CMPLX_ARG (tmp_data), nr,
882 pipvt, tmp_info));
883
884 info = tmp_info;
885 }
886
887 // Throw away extra info LAPACK gives so as to not change output.
888 rcon = 0.0;
889 if (info != 0)
890 info = -1;
891 else if (calc_cond)
892 {
893 F77_INT zgecon_info = 0;
894
895 // Now calculate the condition number for non-singular matrix.
896 char job = '1';
897 Array<double> rz (dim_vector (2 * nc, 1));
898 double *prz = rz.fortran_vec ();
899 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
900 nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
901 rcon, F77_DBLE_CMPLX_ARG (pz), prz, zgecon_info
902 F77_CHAR_ARG_LEN (1)));
903
904 if (zgecon_info != 0)
905 info = -1;
906 }
907
908 if ((info == -1 && ! force)
909 || octave::math::isnan (anorm) || octave::math::isinf (anorm))
910 retval = *this; // Restore contents.
911 else
912 {
913 F77_INT zgetri_info = 0;
914
915 F77_XFCN (zgetri, ZGETRI, (nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, pipvt,
916 F77_DBLE_CMPLX_ARG (pz), lwork, zgetri_info));
917
918 if (zgetri_info != 0)
919 info = -1;
920 }
921
922 if (info != 0)
923 mattype.mark_as_rectangular ();
924
925 return retval;
926 }
927
928 ComplexMatrix
inverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const929 ComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info,
930 double& rcon, bool force, bool calc_cond) const
931 {
932 int typ = mattype.type (false);
933 ComplexMatrix ret;
934
935 if (typ == MatrixType::Unknown)
936 typ = mattype.type (*this);
937
938 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
939 ret = tinverse (mattype, info, rcon, force, calc_cond);
940 else
941 {
942 if (mattype.ishermitian ())
943 {
944 octave::math::chol<ComplexMatrix> chol (*this, info, true, calc_cond);
945 if (info == 0)
946 {
947 if (calc_cond)
948 rcon = chol.rcond ();
949 else
950 rcon = 1.0;
951 ret = chol.inverse ();
952 }
953 else
954 mattype.mark_as_unsymmetric ();
955 }
956
957 if (! mattype.ishermitian ())
958 ret = finverse (mattype, info, rcon, force, calc_cond);
959
960 if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0)
961 {
962 if (numel () == 1)
963 ret = ComplexMatrix (1, 1, 0.0);
964 else
965 ret = ComplexMatrix (rows (), columns (),
966 Complex (octave::numeric_limits<double>::Inf (), 0.0));
967 }
968 }
969
970 return ret;
971 }
972
973 ComplexMatrix
pseudo_inverse(double tol) const974 ComplexMatrix::pseudo_inverse (double tol) const
975 {
976 ComplexMatrix retval;
977
978 octave::math::svd<ComplexMatrix> result (*this,
979 octave::math::svd<ComplexMatrix>::Type::economy);
980
981 DiagMatrix S = result.singular_values ();
982 ComplexMatrix U = result.left_singular_matrix ();
983 ComplexMatrix V = result.right_singular_matrix ();
984
985 ColumnVector sigma = S.extract_diag ();
986
987 octave_idx_type r = sigma.numel () - 1;
988 octave_idx_type nr = rows ();
989 octave_idx_type nc = cols ();
990
991 if (tol <= 0.0)
992 {
993 tol = std::max (nr, nc) * sigma.elem (0)
994 * std::numeric_limits<double>::epsilon ();
995
996 if (tol == 0)
997 tol = std::numeric_limits<double>::min ();
998 }
999
1000 while (r >= 0 && sigma.elem (r) < tol)
1001 r--;
1002
1003 if (r < 0)
1004 retval = ComplexMatrix (nc, nr, 0.0);
1005 else
1006 {
1007 ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
1008 DiagMatrix D = DiagMatrix (sigma.extract (0, r)).inverse ();
1009 ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
1010 retval = Vr * D * Ur.hermitian ();
1011 }
1012
1013 return retval;
1014 }
1015
1016 #if defined (HAVE_FFTW)
1017
1018 ComplexMatrix
fourier(void) const1019 ComplexMatrix::fourier (void) const
1020 {
1021 std::size_t nr = rows ();
1022 std::size_t nc = cols ();
1023
1024 ComplexMatrix retval (nr, nc);
1025
1026 std::size_t npts, nsamples;
1027
1028 if (nr == 1 || nc == 1)
1029 {
1030 npts = (nr > nc ? nr : nc);
1031 nsamples = 1;
1032 }
1033 else
1034 {
1035 npts = nr;
1036 nsamples = nc;
1037 }
1038
1039 const Complex *in (data ());
1040 Complex *out (retval.fortran_vec ());
1041
1042 octave::fftw::fft (in, out, npts, nsamples);
1043
1044 return retval;
1045 }
1046
1047 ComplexMatrix
ifourier(void) const1048 ComplexMatrix::ifourier (void) const
1049 {
1050 std::size_t nr = rows ();
1051 std::size_t nc = cols ();
1052
1053 ComplexMatrix retval (nr, nc);
1054
1055 std::size_t npts, nsamples;
1056
1057 if (nr == 1 || nc == 1)
1058 {
1059 npts = (nr > nc ? nr : nc);
1060 nsamples = 1;
1061 }
1062 else
1063 {
1064 npts = nr;
1065 nsamples = nc;
1066 }
1067
1068 const Complex *in (data ());
1069 Complex *out (retval.fortran_vec ());
1070
1071 octave::fftw::ifft (in, out, npts, nsamples);
1072
1073 return retval;
1074 }
1075
1076 ComplexMatrix
fourier2d(void) const1077 ComplexMatrix::fourier2d (void) const
1078 {
1079 dim_vector dv (rows (), cols ());
1080
1081 ComplexMatrix retval (rows (), cols ());
1082 const Complex *in (data ());
1083 Complex *out (retval.fortran_vec ());
1084
1085 octave::fftw::fftNd (in, out, 2, dv);
1086
1087 return retval;
1088 }
1089
1090 ComplexMatrix
ifourier2d(void) const1091 ComplexMatrix::ifourier2d (void) const
1092 {
1093 dim_vector dv (rows (), cols ());
1094
1095 ComplexMatrix retval (rows (), cols ());
1096 const Complex *in (data ());
1097 Complex *out (retval.fortran_vec ());
1098
1099 octave::fftw::ifftNd (in, out, 2, dv);
1100
1101 return retval;
1102 }
1103
1104 #else
1105
1106 ComplexMatrix
fourier(void) const1107 ComplexMatrix::fourier (void) const
1108 {
1109 (*current_liboctave_error_handler)
1110 ("support for FFTW was unavailable or disabled when liboctave was built");
1111
1112 return ComplexMatrix ();
1113 }
1114
1115 ComplexMatrix
ifourier(void) const1116 ComplexMatrix::ifourier (void) const
1117 {
1118 (*current_liboctave_error_handler)
1119 ("support for FFTW was unavailable or disabled when liboctave was built");
1120
1121 return ComplexMatrix ();
1122 }
1123
1124 ComplexMatrix
fourier2d(void) const1125 ComplexMatrix::fourier2d (void) const
1126 {
1127 (*current_liboctave_error_handler)
1128 ("support for FFTW was unavailable or disabled when liboctave was built");
1129
1130 return ComplexMatrix ();
1131 }
1132
1133 ComplexMatrix
ifourier2d(void) const1134 ComplexMatrix::ifourier2d (void) const
1135 {
1136 (*current_liboctave_error_handler)
1137 ("support for FFTW was unavailable or disabled when liboctave was built");
1138
1139 return ComplexMatrix ();
1140 }
1141
1142 #endif
1143
1144 ComplexDET
determinant(void) const1145 ComplexMatrix::determinant (void) const
1146 {
1147 octave_idx_type info;
1148 double rcon;
1149 return determinant (info, rcon, 0);
1150 }
1151
1152 ComplexDET
determinant(octave_idx_type & info) const1153 ComplexMatrix::determinant (octave_idx_type& info) const
1154 {
1155 double rcon;
1156 return determinant (info, rcon, 0);
1157 }
1158
1159 ComplexDET
determinant(octave_idx_type & info,double & rcon,bool calc_cond) const1160 ComplexMatrix::determinant (octave_idx_type& info, double& rcon,
1161 bool calc_cond) const
1162 {
1163 MatrixType mattype (*this);
1164 return determinant (mattype, info, rcon, calc_cond);
1165 }
1166
1167 ComplexDET
determinant(MatrixType & mattype,octave_idx_type & info,double & rcon,bool calc_cond) const1168 ComplexMatrix::determinant (MatrixType& mattype,
1169 octave_idx_type& info, double& rcon,
1170 bool calc_cond) const
1171 {
1172 ComplexDET retval (1.0);
1173
1174 info = 0;
1175 rcon = 0.0;
1176
1177 F77_INT nr = octave::to_f77_int (rows ());
1178 F77_INT nc = octave::to_f77_int (cols ());
1179
1180 if (nr != nc)
1181 (*current_liboctave_error_handler) ("matrix must be square");
1182
1183 volatile int typ = mattype.type ();
1184
1185 // Even though the matrix is marked as singular (Rectangular), we may
1186 // still get a useful number from the LU factorization, because it always
1187 // completes.
1188
1189 if (typ == MatrixType::Unknown)
1190 typ = mattype.type (*this);
1191 else if (typ == MatrixType::Rectangular)
1192 typ = MatrixType::Full;
1193
1194 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
1195 {
1196 for (F77_INT i = 0; i < nc; i++)
1197 retval *= elem (i,i);
1198 }
1199 else if (typ == MatrixType::Hermitian)
1200 {
1201 ComplexMatrix atmp = *this;
1202 Complex *tmp_data = atmp.fortran_vec ();
1203
1204 double anorm;
1205 if (calc_cond)
1206 anorm = norm1 (*this);
1207
1208 F77_INT tmp_info = 0;
1209
1210 char job = 'L';
1211 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1212 F77_DBLE_CMPLX_ARG (tmp_data), nr, tmp_info
1213 F77_CHAR_ARG_LEN (1)));
1214
1215 info = tmp_info;
1216
1217 if (info != 0)
1218 {
1219 rcon = 0.0;
1220 mattype.mark_as_unsymmetric ();
1221 typ = MatrixType::Full;
1222 }
1223 else
1224 {
1225 if (calc_cond)
1226 {
1227 Array<Complex> z (dim_vector (2 * nc, 1));
1228 Complex *pz = z.fortran_vec ();
1229 Array<double> rz (dim_vector (nc, 1));
1230 double *prz = rz.fortran_vec ();
1231
1232 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1233 nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1234 rcon, F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1235 F77_CHAR_ARG_LEN (1)));
1236
1237 info = tmp_info;
1238
1239 if (info != 0)
1240 rcon = 0.0;
1241 }
1242
1243 for (F77_INT i = 0; i < nc; i++)
1244 retval *= atmp(i,i);
1245
1246 retval = retval.square ();
1247 }
1248 }
1249 else if (typ != MatrixType::Full)
1250 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
1251
1252 if (typ == MatrixType::Full)
1253 {
1254 Array<F77_INT> ipvt (dim_vector (nr, 1));
1255 F77_INT *pipvt = ipvt.fortran_vec ();
1256
1257 ComplexMatrix atmp = *this;
1258 Complex *tmp_data = atmp.fortran_vec ();
1259
1260 info = 0;
1261
1262 // Calculate norm of the matrix (always, see bug #45577) for later use.
1263 double anorm = norm1 (*this);
1264
1265 F77_INT tmp_info = 0;
1266
1267 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1268 if (octave::math::isnan (anorm))
1269 info = -1;
1270 else
1271 {
1272 F77_XFCN (zgetrf, ZGETRF, (nr, nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, pipvt,
1273 tmp_info));
1274
1275 info = tmp_info;
1276 }
1277
1278 // Throw away extra info LAPACK gives so as to not change output.
1279 rcon = 0.0;
1280 if (info != 0)
1281 {
1282 info = -1;
1283 retval = ComplexDET ();
1284 }
1285 else
1286 {
1287 if (calc_cond)
1288 {
1289 // Now calc the condition number for non-singular matrix.
1290 char job = '1';
1291 Array<Complex> z (dim_vector (2 * nc, 1));
1292 Complex *pz = z.fortran_vec ();
1293 Array<double> rz (dim_vector (2 * nc, 1));
1294 double *prz = rz.fortran_vec ();
1295
1296 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1297 nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1298 rcon, F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1299 F77_CHAR_ARG_LEN (1)));
1300
1301 info = tmp_info;
1302 }
1303
1304 if (info != 0)
1305 {
1306 info = -1;
1307 retval = ComplexDET ();
1308 }
1309 else
1310 {
1311 for (F77_INT i = 0; i < nc; i++)
1312 {
1313 Complex c = atmp(i,i);
1314 retval *= (ipvt(i) != (i+1)) ? -c : c;
1315 }
1316 }
1317 }
1318 }
1319
1320 return retval;
1321 }
1322
1323 double
rcond(void) const1324 ComplexMatrix::rcond (void) const
1325 {
1326 MatrixType mattype (*this);
1327 return rcond (mattype);
1328 }
1329
1330 double
rcond(MatrixType & mattype) const1331 ComplexMatrix::rcond (MatrixType& mattype) const
1332 {
1333 double rcon = octave::numeric_limits<double>::NaN ();
1334 F77_INT nr = octave::to_f77_int (rows ());
1335 F77_INT nc = octave::to_f77_int (cols ());
1336
1337 if (nr != nc)
1338 (*current_liboctave_error_handler) ("matrix must be square");
1339
1340 if (nr == 0 || nc == 0)
1341 rcon = octave::numeric_limits<double>::Inf ();
1342 else
1343 {
1344 volatile int typ = mattype.type ();
1345
1346 if (typ == MatrixType::Unknown)
1347 typ = mattype.type (*this);
1348
1349 // Only calculate the condition number for LU/Cholesky
1350 if (typ == MatrixType::Upper)
1351 {
1352 const Complex *tmp_data = fortran_vec ();
1353 F77_INT info = 0;
1354 char norm = '1';
1355 char uplo = 'U';
1356 char dia = 'N';
1357
1358 Array<Complex> z (dim_vector (2 * nc, 1));
1359 Complex *pz = z.fortran_vec ();
1360 Array<double> rz (dim_vector (nc, 1));
1361 double *prz = rz.fortran_vec ();
1362
1363 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1364 F77_CONST_CHAR_ARG2 (&uplo, 1),
1365 F77_CONST_CHAR_ARG2 (&dia, 1),
1366 nr, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr, rcon,
1367 F77_DBLE_CMPLX_ARG (pz), prz, info
1368 F77_CHAR_ARG_LEN (1)
1369 F77_CHAR_ARG_LEN (1)
1370 F77_CHAR_ARG_LEN (1)));
1371
1372 if (info != 0)
1373 rcon = 0;
1374 }
1375 else if (typ == MatrixType::Permuted_Upper)
1376 (*current_liboctave_error_handler)
1377 ("permuted triangular matrix not implemented");
1378 else if (typ == MatrixType::Lower)
1379 {
1380 const Complex *tmp_data = fortran_vec ();
1381 F77_INT info = 0;
1382 char norm = '1';
1383 char uplo = 'L';
1384 char dia = 'N';
1385
1386 Array<Complex> z (dim_vector (2 * nc, 1));
1387 Complex *pz = z.fortran_vec ();
1388 Array<double> rz (dim_vector (nc, 1));
1389 double *prz = rz.fortran_vec ();
1390
1391 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1392 F77_CONST_CHAR_ARG2 (&uplo, 1),
1393 F77_CONST_CHAR_ARG2 (&dia, 1),
1394 nr, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr, rcon,
1395 F77_DBLE_CMPLX_ARG (pz), prz, info
1396 F77_CHAR_ARG_LEN (1)
1397 F77_CHAR_ARG_LEN (1)
1398 F77_CHAR_ARG_LEN (1)));
1399
1400 if (info != 0)
1401 rcon = 0.0;
1402 }
1403 else if (typ == MatrixType::Permuted_Lower)
1404 (*current_liboctave_error_handler)
1405 ("permuted triangular matrix not implemented");
1406 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1407 {
1408 double anorm = -1.0;
1409
1410 if (typ == MatrixType::Hermitian)
1411 {
1412 F77_INT info = 0;
1413 char job = 'L';
1414
1415 ComplexMatrix atmp = *this;
1416 Complex *tmp_data = atmp.fortran_vec ();
1417
1418 anorm = norm1 (atmp);
1419
1420 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1421 F77_DBLE_CMPLX_ARG (tmp_data), nr, info
1422 F77_CHAR_ARG_LEN (1)));
1423
1424 if (info != 0)
1425 {
1426 rcon = 0.0;
1427
1428 mattype.mark_as_unsymmetric ();
1429 typ = MatrixType::Full;
1430 }
1431 else
1432 {
1433 Array<Complex> z (dim_vector (2 * nc, 1));
1434 Complex *pz = z.fortran_vec ();
1435 Array<double> rz (dim_vector (nc, 1));
1436 double *prz = rz.fortran_vec ();
1437
1438 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1439 nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1440 rcon, F77_DBLE_CMPLX_ARG (pz), prz, info
1441 F77_CHAR_ARG_LEN (1)));
1442
1443 if (info != 0)
1444 rcon = 0.0;
1445 }
1446 }
1447
1448 if (typ == MatrixType::Full)
1449 {
1450 F77_INT info = 0;
1451
1452 ComplexMatrix atmp = *this;
1453 Complex *tmp_data = atmp.fortran_vec ();
1454
1455 Array<F77_INT> ipvt (dim_vector (nr, 1));
1456 F77_INT *pipvt = ipvt.fortran_vec ();
1457
1458 if (anorm < 0.0)
1459 anorm = norm1 (atmp);
1460
1461 Array<Complex> z (dim_vector (2 * nc, 1));
1462 Complex *pz = z.fortran_vec ();
1463 Array<double> rz (dim_vector (2 * nc, 1));
1464 double *prz = rz.fortran_vec ();
1465
1466 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1467 if (octave::math::isnan (anorm))
1468 info = -1;
1469 else
1470 F77_XFCN (zgetrf, ZGETRF, (nr, nr,
1471 F77_DBLE_CMPLX_ARG (tmp_data),
1472 nr, pipvt, info));
1473
1474 if (info != 0)
1475 {
1476 rcon = 0.0;
1477 mattype.mark_as_rectangular ();
1478 }
1479 else
1480 {
1481 char job = '1';
1482 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1483 nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1484 rcon, F77_DBLE_CMPLX_ARG (pz), prz, info
1485 F77_CHAR_ARG_LEN (1)));
1486
1487 if (info != 0)
1488 rcon = 0.0;
1489 }
1490 }
1491 }
1492 else
1493 rcon = 0.0;
1494 }
1495
1496 return rcon;
1497 }
1498
1499 ComplexMatrix
utsolve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond,blas_trans_type transt) const1500 ComplexMatrix::utsolve (MatrixType& mattype, const ComplexMatrix& b,
1501 octave_idx_type& info, double& rcon,
1502 solve_singularity_handler sing_handler,
1503 bool calc_cond, blas_trans_type transt) const
1504 {
1505 ComplexMatrix retval;
1506
1507 F77_INT nr = octave::to_f77_int (rows ());
1508 F77_INT nc = octave::to_f77_int (cols ());
1509
1510 F77_INT b_nr = octave::to_f77_int (b.rows ());
1511 F77_INT b_nc = octave::to_f77_int (b.cols ());
1512
1513 if (nr != b_nr)
1514 (*current_liboctave_error_handler)
1515 ("matrix dimension mismatch solution of linear equations");
1516
1517 if (nr == 0 || nc == 0 || b_nc == 0)
1518 retval = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1519 else
1520 {
1521 volatile int typ = mattype.type ();
1522
1523 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1524 (*current_liboctave_error_handler) ("incorrect matrix type");
1525
1526 rcon = 1.0;
1527 info = 0;
1528
1529 if (typ == MatrixType::Permuted_Upper)
1530 (*current_liboctave_error_handler)
1531 ("permuted triangular matrix not implemented");
1532
1533 const Complex *tmp_data = fortran_vec ();
1534
1535 retval = b;
1536 Complex *result = retval.fortran_vec ();
1537
1538 char uplo = 'U';
1539 char trans = get_blas_char (transt);
1540 char dia = 'N';
1541
1542 F77_INT tmp_info = 0;
1543
1544 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1545 F77_CONST_CHAR_ARG2 (&trans, 1),
1546 F77_CONST_CHAR_ARG2 (&dia, 1),
1547 nr, b_nc, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr,
1548 F77_DBLE_CMPLX_ARG (result), nr, tmp_info
1549 F77_CHAR_ARG_LEN (1)
1550 F77_CHAR_ARG_LEN (1)
1551 F77_CHAR_ARG_LEN (1)));
1552
1553 info = tmp_info;
1554
1555 if (calc_cond)
1556 {
1557 char norm = '1';
1558 uplo = 'U';
1559 dia = 'N';
1560
1561 Array<Complex> z (dim_vector (2 * nc, 1));
1562 Complex *pz = z.fortran_vec ();
1563 Array<double> rz (dim_vector (nc, 1));
1564 double *prz = rz.fortran_vec ();
1565
1566 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1567 F77_CONST_CHAR_ARG2 (&uplo, 1),
1568 F77_CONST_CHAR_ARG2 (&dia, 1),
1569 nr, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr, rcon,
1570 F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1571 F77_CHAR_ARG_LEN (1)
1572 F77_CHAR_ARG_LEN (1)
1573 F77_CHAR_ARG_LEN (1)));
1574
1575 info = tmp_info;
1576
1577 if (info != 0)
1578 info = -2;
1579
1580 volatile double rcond_plus_one = rcon + 1.0;
1581
1582 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1583 {
1584 info = -2;
1585
1586 if (sing_handler)
1587 sing_handler (rcon);
1588 else
1589 octave::warn_singular_matrix (rcon);
1590 }
1591 }
1592 }
1593
1594 return retval;
1595 }
1596
1597 ComplexMatrix
ltsolve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond,blas_trans_type transt) const1598 ComplexMatrix::ltsolve (MatrixType& mattype, const ComplexMatrix& b,
1599 octave_idx_type& info, double& rcon,
1600 solve_singularity_handler sing_handler,
1601 bool calc_cond, blas_trans_type transt) const
1602 {
1603 ComplexMatrix retval;
1604
1605 F77_INT nr = octave::to_f77_int (rows ());
1606 F77_INT nc = octave::to_f77_int (cols ());
1607
1608 F77_INT b_nr = octave::to_f77_int (b.rows ());
1609 F77_INT b_nc = octave::to_f77_int (b.cols ());
1610
1611 if (nr != b_nr)
1612 (*current_liboctave_error_handler)
1613 ("matrix dimension mismatch solution of linear equations");
1614
1615 if (nr == 0 || nc == 0 || b_nc == 0)
1616 retval = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1617 else
1618 {
1619 volatile int typ = mattype.type ();
1620
1621 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
1622 (*current_liboctave_error_handler) ("incorrect matrix type");
1623
1624 rcon = 1.0;
1625 info = 0;
1626
1627 if (typ == MatrixType::Permuted_Lower)
1628 (*current_liboctave_error_handler)
1629 ("permuted triangular matrix not implemented");
1630
1631 const Complex *tmp_data = fortran_vec ();
1632
1633 retval = b;
1634 Complex *result = retval.fortran_vec ();
1635
1636 char uplo = 'L';
1637 char trans = get_blas_char (transt);
1638 char dia = 'N';
1639
1640 F77_INT tmp_info = 0;
1641
1642 F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1643 F77_CONST_CHAR_ARG2 (&trans, 1),
1644 F77_CONST_CHAR_ARG2 (&dia, 1),
1645 nr, b_nc, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr,
1646 F77_DBLE_CMPLX_ARG (result), nr, tmp_info
1647 F77_CHAR_ARG_LEN (1)
1648 F77_CHAR_ARG_LEN (1)
1649 F77_CHAR_ARG_LEN (1)));
1650
1651 info = tmp_info;
1652
1653 if (calc_cond)
1654 {
1655 char norm = '1';
1656 uplo = 'L';
1657 dia = 'N';
1658
1659 Array<Complex> z (dim_vector (2 * nc, 1));
1660 Complex *pz = z.fortran_vec ();
1661 Array<double> rz (dim_vector (nc, 1));
1662 double *prz = rz.fortran_vec ();
1663
1664 F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1665 F77_CONST_CHAR_ARG2 (&uplo, 1),
1666 F77_CONST_CHAR_ARG2 (&dia, 1),
1667 nr, F77_CONST_DBLE_CMPLX_ARG (tmp_data), nr, rcon,
1668 F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1669 F77_CHAR_ARG_LEN (1)
1670 F77_CHAR_ARG_LEN (1)
1671 F77_CHAR_ARG_LEN (1)));
1672
1673 info = tmp_info;
1674
1675 if (info != 0)
1676 info = -2;
1677
1678 volatile double rcond_plus_one = rcon + 1.0;
1679
1680 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1681 {
1682 info = -2;
1683
1684 if (sing_handler)
1685 sing_handler (rcon);
1686 else
1687 octave::warn_singular_matrix (rcon);
1688 }
1689 }
1690 }
1691
1692 return retval;
1693 }
1694
1695 ComplexMatrix
fsolve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond) const1696 ComplexMatrix::fsolve (MatrixType& mattype, const ComplexMatrix& b,
1697 octave_idx_type& info, double& rcon,
1698 solve_singularity_handler sing_handler,
1699 bool calc_cond) const
1700 {
1701 ComplexMatrix retval;
1702
1703 F77_INT nr = octave::to_f77_int (rows ());
1704 F77_INT nc = octave::to_f77_int (cols ());
1705
1706 F77_INT b_nr = octave::to_f77_int (b.rows ());
1707 F77_INT b_nc = octave::to_f77_int (b.cols ());
1708
1709 if (nr != nc || nr != b_nr)
1710 (*current_liboctave_error_handler)
1711 ("matrix dimension mismatch solution of linear equations");
1712
1713 if (nr == 0 || b_nc == 0)
1714 retval = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1715 else
1716 {
1717 volatile int typ = mattype.type ();
1718
1719 // Calculate the norm of the matrix for later use when determining rcon.
1720 double anorm = -1.0;
1721
1722 if (typ == MatrixType::Hermitian)
1723 {
1724 info = 0;
1725 char job = 'L';
1726
1727 ComplexMatrix atmp = *this;
1728 Complex *tmp_data = atmp.fortran_vec ();
1729
1730 // The norm of the matrix for later use when determining rcon.
1731 if (calc_cond)
1732 anorm = norm1 (atmp);
1733
1734 F77_INT tmp_info = 0;
1735
1736 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1737 F77_DBLE_CMPLX_ARG (tmp_data), nr, tmp_info
1738 F77_CHAR_ARG_LEN (1)));
1739
1740 info = tmp_info;
1741
1742 // Throw away extra info LAPACK gives so as to not change output.
1743 rcon = 0.0;
1744 if (info != 0)
1745 {
1746 info = -2;
1747
1748 mattype.mark_as_unsymmetric ();
1749 typ = MatrixType::Full;
1750 }
1751 else
1752 {
1753 if (calc_cond)
1754 {
1755 Array<Complex> z (dim_vector (2 * nc, 1));
1756 Complex *pz = z.fortran_vec ();
1757 Array<double> rz (dim_vector (nc, 1));
1758 double *prz = rz.fortran_vec ();
1759
1760 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1761 nr, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1762 rcon, F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1763 F77_CHAR_ARG_LEN (1)));
1764
1765 info = tmp_info;
1766
1767 if (info != 0)
1768 info = -2;
1769
1770 volatile double rcond_plus_one = rcon + 1.0;
1771
1772 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1773 {
1774 info = -2;
1775
1776 if (sing_handler)
1777 sing_handler (rcon);
1778 else
1779 octave::warn_singular_matrix (rcon);
1780 }
1781 }
1782
1783 if (info == 0)
1784 {
1785 retval = b;
1786 Complex *result = retval.fortran_vec ();
1787
1788 F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1789 nr, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), nr,
1790 F77_DBLE_CMPLX_ARG (result), b_nr, tmp_info
1791 F77_CHAR_ARG_LEN (1)));
1792
1793 info = tmp_info;
1794 }
1795 else
1796 {
1797 mattype.mark_as_unsymmetric ();
1798 typ = MatrixType::Full;
1799 }
1800 }
1801 }
1802
1803 if (typ == MatrixType::Full)
1804 {
1805 info = 0;
1806
1807 Array<F77_INT> ipvt (dim_vector (nr, 1));
1808 F77_INT *pipvt = ipvt.fortran_vec ();
1809
1810 ComplexMatrix atmp = *this;
1811 Complex *tmp_data = atmp.fortran_vec ();
1812
1813 Array<Complex> z (dim_vector (2 * nc, 1));
1814 Complex *pz = z.fortran_vec ();
1815 Array<double> rz (dim_vector (2 * nc, 1));
1816 double *prz = rz.fortran_vec ();
1817
1818 // Calculate the norm of the matrix, for later use.
1819 if (calc_cond && anorm < 0.0)
1820 anorm = norm1 (atmp);
1821
1822 F77_INT tmp_info = 0;
1823
1824 // Work around bug #45577, LAPACK crashes Octave if norm is NaN
1825 // and bug #46330, segfault with matrices containing Inf & NaN
1826 if (octave::math::isnan (anorm) || octave::math::isinf (anorm))
1827 info = -2;
1828 else
1829 {
1830 F77_XFCN (zgetrf, ZGETRF, (nr, nr, F77_DBLE_CMPLX_ARG (tmp_data),
1831 nr, pipvt, tmp_info));
1832
1833 info = tmp_info;
1834 }
1835
1836 // Throw away extra info LAPACK gives so as to not change output.
1837 rcon = 0.0;
1838 if (info != 0)
1839 {
1840 info = -2;
1841
1842 if (sing_handler)
1843 sing_handler (rcon);
1844 else
1845 octave::warn_singular_matrix ();
1846
1847 mattype.mark_as_rectangular ();
1848 }
1849 else
1850 {
1851 if (calc_cond)
1852 {
1853 // Calculate the condition number for non-singular matrix.
1854 char job = '1';
1855 F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1856 nc, F77_DBLE_CMPLX_ARG (tmp_data), nr, anorm,
1857 rcon, F77_DBLE_CMPLX_ARG (pz), prz, tmp_info
1858 F77_CHAR_ARG_LEN (1)));
1859
1860 info = tmp_info;
1861
1862 if (info != 0)
1863 info = -2;
1864
1865 volatile double rcond_plus_one = rcon + 1.0;
1866
1867 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1868 {
1869 info = -2;
1870
1871 if (sing_handler)
1872 sing_handler (rcon);
1873 else
1874 octave::warn_singular_matrix (rcon);
1875 }
1876 }
1877
1878 if (info == 0)
1879 {
1880 retval = b;
1881 Complex *result = retval.fortran_vec ();
1882
1883 char job = 'N';
1884 F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1885 nr, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), nr,
1886 pipvt, F77_DBLE_CMPLX_ARG (result), b_nr, tmp_info
1887 F77_CHAR_ARG_LEN (1)));
1888
1889 info = tmp_info;
1890 }
1891 else
1892 mattype.mark_as_rectangular ();
1893 }
1894 }
1895
1896 if (octave::math::isinf (anorm))
1897 {
1898 retval = ComplexMatrix (b_nr, b_nc, Complex (0, 0));
1899 mattype.mark_as_full ();
1900 }
1901 }
1902
1903 return retval;
1904 }
1905
1906 ComplexMatrix
solve(MatrixType & mattype,const Matrix & b) const1907 ComplexMatrix::solve (MatrixType& mattype, const Matrix& b) const
1908 {
1909 octave_idx_type info;
1910 double rcon;
1911 return solve (mattype, b, info, rcon, nullptr);
1912 }
1913
1914 ComplexMatrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info) const1915 ComplexMatrix::solve (MatrixType& mattype, const Matrix& b,
1916 octave_idx_type& info) const
1917 {
1918 double rcon;
1919 return solve (mattype, b, info, rcon, nullptr);
1920 }
1921
1922 ComplexMatrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon) const1923 ComplexMatrix::solve (MatrixType& mattype, const Matrix& b,
1924 octave_idx_type& info, double& rcon) const
1925 {
1926 return solve (mattype, b, info, rcon, nullptr);
1927 }
1928
1929 ComplexMatrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool singular_fallback,blas_trans_type transt) const1930 ComplexMatrix::solve (MatrixType& mattype, const Matrix& b,
1931 octave_idx_type& info, double& rcon,
1932 solve_singularity_handler sing_handler,
1933 bool singular_fallback, blas_trans_type transt) const
1934 {
1935 ComplexMatrix tmp (b);
1936 return solve (mattype, tmp, info, rcon, sing_handler, singular_fallback,
1937 transt);
1938 }
1939
1940 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b) const1941 ComplexMatrix::solve (MatrixType& mattype, const ComplexMatrix& b) const
1942 {
1943 octave_idx_type info;
1944 double rcon;
1945 return solve (mattype, b, info, rcon, nullptr);
1946 }
1947
1948 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info) const1949 ComplexMatrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1950 octave_idx_type& info) const
1951 {
1952 double rcon;
1953 return solve (mattype, b, info, rcon, nullptr);
1954 }
1955
1956 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon) const1957 ComplexMatrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1958 octave_idx_type& info, double& rcon) const
1959 {
1960 return solve (mattype, b, info, rcon, nullptr);
1961 }
1962
1963 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool singular_fallback,blas_trans_type transt) const1964 ComplexMatrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1965 octave_idx_type& info, double& rcon,
1966 solve_singularity_handler sing_handler,
1967 bool singular_fallback, blas_trans_type transt) const
1968 {
1969 ComplexMatrix retval;
1970 int typ = mattype.type ();
1971
1972 if (typ == MatrixType::Unknown)
1973 typ = mattype.type (*this);
1974
1975 // Only calculate the condition number for LU/Cholesky
1976 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1977 retval = utsolve (mattype, b, info, rcon, sing_handler, true, transt);
1978 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1979 retval = ltsolve (mattype, b, info, rcon, sing_handler, true, transt);
1980 else if (transt == blas_trans)
1981 return transpose ().solve (mattype, b, info, rcon, sing_handler,
1982 singular_fallback);
1983 else if (transt == blas_conj_trans)
1984 retval = hermitian ().solve (mattype, b, info, rcon, sing_handler,
1985 singular_fallback);
1986 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1987 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
1988 else if (typ != MatrixType::Rectangular)
1989 (*current_liboctave_error_handler) ("unknown matrix type");
1990
1991 // Rectangular or one of the above solvers flags a singular matrix
1992 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
1993 {
1994 octave_idx_type rank;
1995 retval = lssolve (b, info, rank, rcon);
1996 }
1997
1998 return retval;
1999 }
2000
2001 ComplexColumnVector
solve(MatrixType & mattype,const ColumnVector & b) const2002 ComplexMatrix::solve (MatrixType& mattype, const ColumnVector& b) const
2003 {
2004 octave_idx_type info;
2005 double rcon;
2006 return solve (mattype, ComplexColumnVector (b), info, rcon, nullptr);
2007 }
2008
2009 ComplexColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info) const2010 ComplexMatrix::solve (MatrixType& mattype, const ColumnVector& b,
2011 octave_idx_type& info) const
2012 {
2013 double rcon;
2014 return solve (mattype, ComplexColumnVector (b), info, rcon, nullptr);
2015 }
2016
2017 ComplexColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info,double & rcon) const2018 ComplexMatrix::solve (MatrixType& mattype, const ColumnVector& b,
2019 octave_idx_type& info, double& rcon) const
2020 {
2021 return solve (mattype, ComplexColumnVector (b), info, rcon, nullptr);
2022 }
2023
2024 ComplexColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2025 ComplexMatrix::solve (MatrixType& mattype, const ColumnVector& b,
2026 octave_idx_type& info, double& rcon,
2027 solve_singularity_handler sing_handler,
2028 blas_trans_type transt) const
2029 {
2030 return solve (mattype, ComplexColumnVector (b), info, rcon, sing_handler,
2031 transt);
2032 }
2033
2034 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b) const2035 ComplexMatrix::solve (MatrixType& mattype, const ComplexColumnVector& b) const
2036 {
2037 octave_idx_type info;
2038 double rcon;
2039 return solve (mattype, b, info, rcon, nullptr);
2040 }
2041
2042 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info) const2043 ComplexMatrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
2044 octave_idx_type& info) const
2045 {
2046 double rcon;
2047 return solve (mattype, b, info, rcon, nullptr);
2048 }
2049
2050 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info,double & rcon) const2051 ComplexMatrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
2052 octave_idx_type& info, double& rcon) const
2053 {
2054 return solve (mattype, b, info, rcon, nullptr);
2055 }
2056
2057 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2058 ComplexMatrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
2059 octave_idx_type& info, double& rcon,
2060 solve_singularity_handler sing_handler,
2061 blas_trans_type transt) const
2062 {
2063
2064 ComplexMatrix tmp (b);
2065 tmp = solve (mattype, tmp, info, rcon, sing_handler, true, transt);
2066 return tmp.column (static_cast<octave_idx_type> (0));
2067 }
2068
2069 ComplexMatrix
solve(const Matrix & b) const2070 ComplexMatrix::solve (const Matrix& b) const
2071 {
2072 octave_idx_type info;
2073 double rcon;
2074 return solve (b, info, rcon, nullptr);
2075 }
2076
2077 ComplexMatrix
solve(const Matrix & b,octave_idx_type & info) const2078 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
2079 {
2080 double rcon;
2081 return solve (b, info, rcon, nullptr);
2082 }
2083
2084 ComplexMatrix
solve(const Matrix & b,octave_idx_type & info,double & rcon) const2085 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info,
2086 double& rcon) const
2087 {
2088 return solve (b, info, rcon, nullptr);
2089 }
2090
2091 ComplexMatrix
solve(const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2092 ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcon,
2093 solve_singularity_handler sing_handler,
2094 blas_trans_type transt) const
2095 {
2096 ComplexMatrix tmp (b);
2097 return solve (tmp, info, rcon, sing_handler, transt);
2098 }
2099
2100 ComplexMatrix
solve(const ComplexMatrix & b) const2101 ComplexMatrix::solve (const ComplexMatrix& b) const
2102 {
2103 octave_idx_type info;
2104 double rcon;
2105 return solve (b, info, rcon, nullptr);
2106 }
2107
2108 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info) const2109 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
2110 {
2111 double rcon;
2112 return solve (b, info, rcon, nullptr);
2113 }
2114
2115 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info,double & rcon) const2116 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info,
2117 double& rcon) const
2118 {
2119 return solve (b, info, rcon, nullptr);
2120 }
2121
2122 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2123 ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info,
2124 double& rcon,
2125 solve_singularity_handler sing_handler,
2126 blas_trans_type transt) const
2127 {
2128 MatrixType mattype (*this);
2129 return solve (mattype, b, info, rcon, sing_handler, true, transt);
2130 }
2131
2132 ComplexColumnVector
solve(const ColumnVector & b) const2133 ComplexMatrix::solve (const ColumnVector& b) const
2134 {
2135 octave_idx_type info;
2136 double rcon;
2137 return solve (ComplexColumnVector (b), info, rcon, nullptr);
2138 }
2139
2140 ComplexColumnVector
solve(const ColumnVector & b,octave_idx_type & info) const2141 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const
2142 {
2143 double rcon;
2144 return solve (ComplexColumnVector (b), info, rcon, nullptr);
2145 }
2146
2147 ComplexColumnVector
solve(const ColumnVector & b,octave_idx_type & info,double & rcon) const2148 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2149 double& rcon) const
2150 {
2151 return solve (ComplexColumnVector (b), info, rcon, nullptr);
2152 }
2153
2154 ComplexColumnVector
solve(const ColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2155 ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
2156 double& rcon,
2157 solve_singularity_handler sing_handler,
2158 blas_trans_type transt) const
2159 {
2160 return solve (ComplexColumnVector (b), info, rcon, sing_handler, transt);
2161 }
2162
2163 ComplexColumnVector
solve(const ComplexColumnVector & b) const2164 ComplexMatrix::solve (const ComplexColumnVector& b) const
2165 {
2166 octave_idx_type info;
2167 double rcon;
2168 return solve (b, info, rcon, nullptr);
2169 }
2170
2171 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info) const2172 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const
2173 {
2174 double rcon;
2175 return solve (b, info, rcon, nullptr);
2176 }
2177
2178 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info,double & rcon) const2179 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
2180 double& rcon) const
2181 {
2182 return solve (b, info, rcon, nullptr);
2183 }
2184
2185 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const2186 ComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
2187 double& rcon,
2188 solve_singularity_handler sing_handler,
2189 blas_trans_type transt) const
2190 {
2191 MatrixType mattype (*this);
2192 return solve (mattype, b, info, rcon, sing_handler, transt);
2193 }
2194
2195 ComplexMatrix
lssolve(const Matrix & b) const2196 ComplexMatrix::lssolve (const Matrix& b) const
2197 {
2198 octave_idx_type info;
2199 octave_idx_type rank;
2200 double rcon;
2201 return lssolve (ComplexMatrix (b), info, rank, rcon);
2202 }
2203
2204 ComplexMatrix
lssolve(const Matrix & b,octave_idx_type & info) const2205 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info) const
2206 {
2207 octave_idx_type rank;
2208 double rcon;
2209 return lssolve (ComplexMatrix (b), info, rank, rcon);
2210 }
2211
2212 ComplexMatrix
lssolve(const Matrix & b,octave_idx_type & info,octave_idx_type & rank) const2213 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
2214 octave_idx_type& rank) const
2215 {
2216 double rcon;
2217 return lssolve (ComplexMatrix (b), info, rank, rcon);
2218 }
2219
2220 ComplexMatrix
lssolve(const Matrix & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2221 ComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info,
2222 octave_idx_type& rank, double& rcon) const
2223 {
2224 return lssolve (ComplexMatrix (b), info, rank, rcon);
2225 }
2226
2227 ComplexMatrix
lssolve(const ComplexMatrix & b) const2228 ComplexMatrix::lssolve (const ComplexMatrix& b) const
2229 {
2230 octave_idx_type info;
2231 octave_idx_type rank;
2232 double rcon;
2233 return lssolve (b, info, rank, rcon);
2234 }
2235
2236 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info) const2237 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const
2238 {
2239 octave_idx_type rank;
2240 double rcon;
2241 return lssolve (b, info, rank, rcon);
2242 }
2243
2244 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info,octave_idx_type & rank) const2245 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2246 octave_idx_type& rank) const
2247 {
2248 double rcon;
2249 return lssolve (b, info, rank, rcon);
2250 }
2251
2252 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2253 ComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2254 octave_idx_type& rank, double& rcon) const
2255 {
2256 ComplexMatrix retval;
2257
2258 F77_INT m = octave::to_f77_int (rows ());
2259 F77_INT n = octave::to_f77_int (cols ());
2260
2261 F77_INT b_nr = octave::to_f77_int (b.rows ());
2262 F77_INT b_nc = octave::to_f77_int (b.cols ());
2263 F77_INT nrhs = b_nc; // alias for code readability
2264
2265 if (m != b_nr)
2266 (*current_liboctave_error_handler)
2267 ("matrix dimension mismatch solution of linear equations");
2268
2269 if (m == 0 || n == 0 || b_nc == 0)
2270 retval = ComplexMatrix (n, b_nc, Complex (0.0, 0.0));
2271 else
2272 {
2273 volatile F77_INT minmn = (m < n ? m : n);
2274 F77_INT maxmn = (m > n ? m : n);
2275 rcon = -1.0;
2276
2277 if (m != n)
2278 {
2279 retval = ComplexMatrix (maxmn, nrhs);
2280
2281 for (F77_INT j = 0; j < nrhs; j++)
2282 for (F77_INT i = 0; i < m; i++)
2283 retval.elem (i, j) = b.elem (i, j);
2284 }
2285 else
2286 retval = b;
2287
2288 ComplexMatrix atmp = *this;
2289 Complex *tmp_data = atmp.fortran_vec ();
2290
2291 Complex *pretval = retval.fortran_vec ();
2292 Array<double> s (dim_vector (minmn, 1));
2293 double *ps = s.fortran_vec ();
2294
2295 // Ask ZGELSD what the dimension of WORK should be.
2296 F77_INT lwork = -1;
2297
2298 Array<Complex> work (dim_vector (1, 1));
2299
2300 F77_INT smlsiz;
2301 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6),
2302 F77_CONST_CHAR_ARG2 (" ", 1),
2303 0, 0, 0, 0, smlsiz
2304 F77_CHAR_ARG_LEN (6)
2305 F77_CHAR_ARG_LEN (1));
2306
2307 F77_INT mnthr;
2308 F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("ZGELSD", 6),
2309 F77_CONST_CHAR_ARG2 (" ", 1),
2310 m, n, nrhs, -1, mnthr
2311 F77_CHAR_ARG_LEN (6)
2312 F77_CHAR_ARG_LEN (1));
2313
2314 // We compute the size of rwork and iwork because ZGELSD in
2315 // older versions of LAPACK does not return them on a query
2316 // call.
2317 double dminmn = static_cast<double> (minmn);
2318 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2319 double tmp = octave::math::log2 (dminmn / dsmlsizp1);
2320
2321 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
2322 if (nlvl < 0)
2323 nlvl = 0;
2324
2325 F77_INT lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2326 + 3*smlsiz*nrhs
2327 + std::max ((smlsiz+1)*(smlsiz+1),
2328 n*(1+nrhs) + 2*nrhs);
2329 if (lrwork < 1)
2330 lrwork = 1;
2331 Array<double> rwork (dim_vector (lrwork, 1));
2332 double *prwork = rwork.fortran_vec ();
2333
2334 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
2335 if (liwork < 1)
2336 liwork = 1;
2337 Array<F77_INT> iwork (dim_vector (liwork, 1));
2338 F77_INT *piwork = iwork.fortran_vec ();
2339
2340 F77_INT tmp_info = 0;
2341 F77_INT tmp_rank = 0;
2342
2343 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data), m,
2344 F77_DBLE_CMPLX_ARG (pretval), maxmn,
2345 ps, rcon, tmp_rank, F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
2346 lwork, prwork, piwork, tmp_info));
2347
2348 info = tmp_info;
2349 rank = tmp_rank;
2350
2351 // The workspace query is broken in at least LAPACK 3.0.0
2352 // through 3.1.1 when n >= mnthr. The obtuse formula below
2353 // should provide sufficient workspace for ZGELSD to operate
2354 // efficiently.
2355 if (n > m && n >= mnthr)
2356 {
2357 F77_INT addend = m;
2358
2359 if (2*m-4 > addend)
2360 addend = 2*m-4;
2361
2362 if (nrhs > addend)
2363 addend = nrhs;
2364
2365 if (n-3*m > addend)
2366 addend = n-3*m;
2367
2368 const F77_INT lworkaround = 4*m + m*m + addend;
2369
2370 if (std::real (work(0)) < lworkaround)
2371 work(0) = lworkaround;
2372 }
2373 else if (m >= n)
2374 {
2375 F77_INT lworkaround = 2*m + m*nrhs;
2376
2377 if (std::real (work(0)) < lworkaround)
2378 work(0) = lworkaround;
2379 }
2380
2381 lwork = static_cast<F77_INT> (std::real (work(0)));
2382 work.resize (dim_vector (lwork, 1));
2383
2384 double anorm = norm1 (*this);
2385
2386 if (octave::math::isinf (anorm))
2387 {
2388 rcon = 0.0;
2389 retval = ComplexMatrix (n, b_nc, 0.0);
2390 }
2391 else if (octave::math::isnan (anorm))
2392 {
2393 rcon = octave::numeric_limits<double>::NaN ();
2394 retval = ComplexMatrix (n, b_nc,
2395 octave::numeric_limits<double>::NaN ());
2396 }
2397 else
2398 {
2399 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data),
2400 m, F77_DBLE_CMPLX_ARG (pretval),
2401 maxmn, ps, rcon, tmp_rank,
2402 F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
2403 lwork, prwork, piwork, tmp_info));
2404
2405 info = tmp_info;
2406 rank = tmp_rank;
2407
2408 if (s.elem (0) == 0.0)
2409 rcon = 0.0;
2410 else
2411 rcon = s.elem (minmn - 1) / s.elem (0);
2412
2413 retval.resize (n, nrhs);
2414 }
2415 }
2416
2417 return retval;
2418 }
2419
2420 ComplexColumnVector
lssolve(const ColumnVector & b) const2421 ComplexMatrix::lssolve (const ColumnVector& b) const
2422 {
2423 octave_idx_type info;
2424 octave_idx_type rank;
2425 double rcon;
2426 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2427 }
2428
2429 ComplexColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info) const2430 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
2431 {
2432 octave_idx_type rank;
2433 double rcon;
2434 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2435 }
2436
2437 ComplexColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info,octave_idx_type & rank) const2438 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2439 octave_idx_type& rank) const
2440 {
2441 double rcon;
2442 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2443 }
2444
2445 ComplexColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2446 ComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2447 octave_idx_type& rank, double& rcon) const
2448 {
2449 return lssolve (ComplexColumnVector (b), info, rank, rcon);
2450 }
2451
2452 ComplexColumnVector
lssolve(const ComplexColumnVector & b) const2453 ComplexMatrix::lssolve (const ComplexColumnVector& b) const
2454 {
2455 octave_idx_type info;
2456 octave_idx_type rank;
2457 double rcon;
2458 return lssolve (b, info, rank, rcon);
2459 }
2460
2461 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info) const2462 ComplexMatrix::lssolve (const ComplexColumnVector& b,
2463 octave_idx_type& info) const
2464 {
2465 octave_idx_type rank;
2466 double rcon;
2467 return lssolve (b, info, rank, rcon);
2468 }
2469
2470 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info,octave_idx_type & rank) const2471 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2472 octave_idx_type& rank) const
2473 {
2474 double rcon;
2475 return lssolve (b, info, rank, rcon);
2476
2477 }
2478
2479 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2480 ComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2481 octave_idx_type& rank, double& rcon) const
2482 {
2483 ComplexColumnVector retval;
2484
2485 F77_INT nrhs = 1;
2486
2487 F77_INT m = octave::to_f77_int (rows ());
2488 F77_INT n = octave::to_f77_int (cols ());
2489
2490 F77_INT b_nel = octave::to_f77_int (b.numel ());
2491
2492 if (m != b_nel)
2493 (*current_liboctave_error_handler)
2494 ("matrix dimension mismatch solution of linear equations");
2495
2496 if (m == 0 || n == 0)
2497 retval = ComplexColumnVector (n, Complex (0.0, 0.0));
2498 else
2499 {
2500 volatile F77_INT minmn = (m < n ? m : n);
2501 F77_INT maxmn = (m > n ? m : n);
2502 rcon = -1.0;
2503
2504 if (m != n)
2505 {
2506 retval = ComplexColumnVector (maxmn);
2507
2508 for (F77_INT i = 0; i < m; i++)
2509 retval.elem (i) = b.elem (i);
2510 }
2511 else
2512 retval = b;
2513
2514 ComplexMatrix atmp = *this;
2515 Complex *tmp_data = atmp.fortran_vec ();
2516
2517 Complex *pretval = retval.fortran_vec ();
2518 Array<double> s (dim_vector (minmn, 1));
2519 double *ps = s.fortran_vec ();
2520
2521 // Ask ZGELSD what the dimension of WORK should be.
2522 F77_INT lwork = -1;
2523
2524 Array<Complex> work (dim_vector (1, 1));
2525
2526 F77_INT smlsiz;
2527 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("ZGELSD", 6),
2528 F77_CONST_CHAR_ARG2 (" ", 1),
2529 0, 0, 0, 0, smlsiz
2530 F77_CHAR_ARG_LEN (6)
2531 F77_CHAR_ARG_LEN (1));
2532
2533 // We compute the size of rwork and iwork because ZGELSD in
2534 // older versions of LAPACK does not return them on a query
2535 // call.
2536 double dminmn = static_cast<double> (minmn);
2537 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2538 double tmp = octave::math::log2 (dminmn / dsmlsizp1);
2539
2540 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
2541 if (nlvl < 0)
2542 nlvl = 0;
2543
2544 F77_INT lrwork = minmn*(10 + 2*smlsiz + 8*nlvl)
2545 + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
2546 if (lrwork < 1)
2547 lrwork = 1;
2548 Array<double> rwork (dim_vector (lrwork, 1));
2549 double *prwork = rwork.fortran_vec ();
2550
2551 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
2552 if (liwork < 1)
2553 liwork = 1;
2554 Array<F77_INT> iwork (dim_vector (liwork, 1));
2555 F77_INT *piwork = iwork.fortran_vec ();
2556
2557 F77_INT tmp_info = 0;
2558 F77_INT tmp_rank = 0;
2559
2560 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data), m,
2561 F77_DBLE_CMPLX_ARG (pretval), maxmn,
2562 ps, rcon, tmp_rank, F77_DBLE_CMPLX_ARG (work.fortran_vec ()),
2563 lwork, prwork, piwork, tmp_info));
2564
2565 info = tmp_info;
2566 rank = tmp_rank;
2567
2568 lwork = static_cast<F77_INT> (std::real (work(0)));
2569 work.resize (dim_vector (lwork, 1));
2570 rwork.resize (dim_vector (static_cast<F77_INT> (rwork(0)), 1));
2571 iwork.resize (dim_vector (iwork(0), 1));
2572
2573 F77_XFCN (zgelsd, ZGELSD, (m, n, nrhs, F77_DBLE_CMPLX_ARG (tmp_data), m,
2574 F77_DBLE_CMPLX_ARG (pretval),
2575 maxmn, ps, rcon, tmp_rank,
2576 F77_DBLE_CMPLX_ARG (work.fortran_vec ()), lwork,
2577 prwork, piwork, tmp_info));
2578
2579 info = tmp_info;
2580 rank = tmp_rank;
2581
2582 if (rank < minmn)
2583 {
2584 if (s.elem (0) == 0.0)
2585 rcon = 0.0;
2586 else
2587 rcon = s.elem (minmn - 1) / s.elem (0);
2588
2589 retval.resize (n);
2590 }
2591 }
2592
2593 return retval;
2594 }
2595
2596 // column vector by row vector -> matrix operations
2597
2598 ComplexMatrix
operator *(const ColumnVector & v,const ComplexRowVector & a)2599 operator * (const ColumnVector& v, const ComplexRowVector& a)
2600 {
2601 ComplexColumnVector tmp (v);
2602 return tmp * a;
2603 }
2604
2605 ComplexMatrix
operator *(const ComplexColumnVector & a,const RowVector & b)2606 operator * (const ComplexColumnVector& a, const RowVector& b)
2607 {
2608 ComplexRowVector tmp (b);
2609 return a * tmp;
2610 }
2611
2612 ComplexMatrix
operator *(const ComplexColumnVector & v,const ComplexRowVector & a)2613 operator * (const ComplexColumnVector& v, const ComplexRowVector& a)
2614 {
2615 ComplexMatrix retval;
2616
2617 F77_INT len = octave::to_f77_int (v.numel ());
2618
2619 if (len != 0)
2620 {
2621 F77_INT a_len = octave::to_f77_int (a.numel ());
2622
2623 retval = ComplexMatrix (len, a_len);
2624 Complex *c = retval.fortran_vec ();
2625
2626 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
2627 F77_CONST_CHAR_ARG2 ("N", 1),
2628 len, a_len, 1, 1.0, F77_CONST_DBLE_CMPLX_ARG (v.data ()), len,
2629 F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1, 0.0, F77_DBLE_CMPLX_ARG (c), len
2630 F77_CHAR_ARG_LEN (1)
2631 F77_CHAR_ARG_LEN (1)));
2632 }
2633
2634 return retval;
2635 }
2636
2637 // matrix by diagonal matrix -> matrix operations
2638
2639 ComplexMatrix&
operator +=(const DiagMatrix & a)2640 ComplexMatrix::operator += (const DiagMatrix& a)
2641 {
2642 octave_idx_type nr = rows ();
2643 octave_idx_type nc = cols ();
2644
2645 octave_idx_type a_nr = a.rows ();
2646 octave_idx_type a_nc = a.cols ();
2647
2648 if (nr != a_nr || nc != a_nc)
2649 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2650
2651 for (octave_idx_type i = 0; i < a.length (); i++)
2652 elem (i, i) += a.elem (i, i);
2653
2654 return *this;
2655 }
2656
2657 ComplexMatrix&
operator -=(const DiagMatrix & a)2658 ComplexMatrix::operator -= (const DiagMatrix& a)
2659 {
2660 octave_idx_type nr = rows ();
2661 octave_idx_type nc = cols ();
2662
2663 octave_idx_type a_nr = a.rows ();
2664 octave_idx_type a_nc = a.cols ();
2665
2666 if (nr != a_nr || nc != a_nc)
2667 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2668
2669 for (octave_idx_type i = 0; i < a.length (); i++)
2670 elem (i, i) -= a.elem (i, i);
2671
2672 return *this;
2673 }
2674
2675 ComplexMatrix&
operator +=(const ComplexDiagMatrix & a)2676 ComplexMatrix::operator += (const ComplexDiagMatrix& a)
2677 {
2678 octave_idx_type nr = rows ();
2679 octave_idx_type nc = cols ();
2680
2681 octave_idx_type a_nr = a.rows ();
2682 octave_idx_type a_nc = a.cols ();
2683
2684 if (nr != a_nr || nc != a_nc)
2685 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2686
2687 for (octave_idx_type i = 0; i < a.length (); i++)
2688 elem (i, i) += a.elem (i, i);
2689
2690 return *this;
2691 }
2692
2693 ComplexMatrix&
operator -=(const ComplexDiagMatrix & a)2694 ComplexMatrix::operator -= (const ComplexDiagMatrix& a)
2695 {
2696 octave_idx_type nr = rows ();
2697 octave_idx_type nc = cols ();
2698
2699 octave_idx_type a_nr = a.rows ();
2700 octave_idx_type a_nc = a.cols ();
2701
2702 if (nr != a_nr || nc != a_nc)
2703 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2704
2705 for (octave_idx_type i = 0; i < a.length (); i++)
2706 elem (i, i) -= a.elem (i, i);
2707
2708 return *this;
2709 }
2710
2711 // matrix by matrix -> matrix operations
2712
2713 ComplexMatrix&
operator +=(const Matrix & a)2714 ComplexMatrix::operator += (const Matrix& a)
2715 {
2716 octave_idx_type nr = rows ();
2717 octave_idx_type nc = cols ();
2718
2719 octave_idx_type a_nr = a.rows ();
2720 octave_idx_type a_nc = a.cols ();
2721
2722 if (nr != a_nr || nc != a_nc)
2723 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2724
2725 if (nr == 0 || nc == 0)
2726 return *this;
2727
2728 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2729
2730 mx_inline_add2 (numel (), d, a.data ());
2731 return *this;
2732 }
2733
2734 ComplexMatrix&
operator -=(const Matrix & a)2735 ComplexMatrix::operator -= (const Matrix& a)
2736 {
2737 octave_idx_type nr = rows ();
2738 octave_idx_type nc = cols ();
2739
2740 octave_idx_type a_nr = a.rows ();
2741 octave_idx_type a_nc = a.cols ();
2742
2743 if (nr != a_nr || nc != a_nc)
2744 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2745
2746 if (nr == 0 || nc == 0)
2747 return *this;
2748
2749 Complex *d = fortran_vec (); // Ensures only one reference to my privates!
2750
2751 mx_inline_sub2 (numel (), d, a.data ());
2752 return *this;
2753 }
2754
2755 // other operations
2756
2757 boolMatrix
all(int dim) const2758 ComplexMatrix::all (int dim) const
2759 {
2760 return ComplexNDArray::all (dim);
2761 }
2762
2763 boolMatrix
any(int dim) const2764 ComplexMatrix::any (int dim) const
2765 {
2766 return ComplexNDArray::any (dim);
2767 }
2768
2769 ComplexMatrix
cumprod(int dim) const2770 ComplexMatrix::cumprod (int dim) const
2771 {
2772 return ComplexNDArray::cumprod (dim);
2773 }
2774
2775 ComplexMatrix
cumsum(int dim) const2776 ComplexMatrix::cumsum (int dim) const
2777 {
2778 return ComplexNDArray::cumsum (dim);
2779 }
2780
2781 ComplexMatrix
prod(int dim) const2782 ComplexMatrix::prod (int dim) const
2783 {
2784 return ComplexNDArray::prod (dim);
2785 }
2786
2787 ComplexMatrix
sum(int dim) const2788 ComplexMatrix::sum (int dim) const
2789 {
2790 return ComplexNDArray::sum (dim);
2791 }
2792
2793 ComplexMatrix
sumsq(int dim) const2794 ComplexMatrix::sumsq (int dim) const
2795 {
2796 return ComplexNDArray::sumsq (dim);
2797 }
2798
2799 Matrix
abs(void) const2800 ComplexMatrix::abs (void) const
2801 {
2802 return ComplexNDArray::abs ();
2803 }
2804
2805 ComplexMatrix
diag(octave_idx_type k) const2806 ComplexMatrix::diag (octave_idx_type k) const
2807 {
2808 return ComplexNDArray::diag (k);
2809 }
2810
2811 ComplexDiagMatrix
diag(octave_idx_type m,octave_idx_type n) const2812 ComplexMatrix::diag (octave_idx_type m, octave_idx_type n) const
2813 {
2814 octave_idx_type nr = rows ();
2815 octave_idx_type nc = cols ();
2816
2817 if (nr != 1 && nc != 1)
2818 (*current_liboctave_error_handler) ("diag: expecting vector argument");
2819
2820 return ComplexDiagMatrix (*this, m, n);
2821 }
2822
2823 bool
row_is_real_only(octave_idx_type i) const2824 ComplexMatrix::row_is_real_only (octave_idx_type i) const
2825 {
2826 bool retval = true;
2827
2828 octave_idx_type nc = columns ();
2829
2830 for (octave_idx_type j = 0; j < nc; j++)
2831 {
2832 if (std::imag (elem (i, j)) != 0.0)
2833 {
2834 retval = false;
2835 break;
2836 }
2837 }
2838
2839 return retval;
2840 }
2841
2842 bool
column_is_real_only(octave_idx_type j) const2843 ComplexMatrix::column_is_real_only (octave_idx_type j) const
2844 {
2845 bool retval = true;
2846
2847 octave_idx_type nr = rows ();
2848
2849 for (octave_idx_type i = 0; i < nr; i++)
2850 {
2851 if (std::imag (elem (i, j)) != 0.0)
2852 {
2853 retval = false;
2854 break;
2855 }
2856 }
2857
2858 return retval;
2859 }
2860
2861 ComplexColumnVector
row_min(void) const2862 ComplexMatrix::row_min (void) const
2863 {
2864 Array<octave_idx_type> dummy_idx;
2865 return row_min (dummy_idx);
2866 }
2867
2868 ComplexColumnVector
row_min(Array<octave_idx_type> & idx_arg) const2869 ComplexMatrix::row_min (Array<octave_idx_type>& idx_arg) const
2870 {
2871 ComplexColumnVector result;
2872
2873 octave_idx_type nr = rows ();
2874 octave_idx_type nc = cols ();
2875
2876 if (nr > 0 && nc > 0)
2877 {
2878 result.resize (nr);
2879 idx_arg.resize (dim_vector (nr, 1));
2880
2881 for (octave_idx_type i = 0; i < nr; i++)
2882 {
2883 bool real_only = row_is_real_only (i);
2884
2885 octave_idx_type idx_j;
2886
2887 Complex tmp_min;
2888
2889 double abs_min = octave::numeric_limits<double>::NaN ();
2890
2891 for (idx_j = 0; idx_j < nc; idx_j++)
2892 {
2893 tmp_min = elem (i, idx_j);
2894
2895 if (! octave::math::isnan (tmp_min))
2896 {
2897 abs_min = (real_only ? tmp_min.real ()
2898 : std::abs (tmp_min));
2899 break;
2900 }
2901 }
2902
2903 for (octave_idx_type j = idx_j+1; j < nc; j++)
2904 {
2905 Complex tmp = elem (i, j);
2906
2907 if (octave::math::isnan (tmp))
2908 continue;
2909
2910 double abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
2911
2912 if (abs_tmp < abs_min)
2913 {
2914 idx_j = j;
2915 tmp_min = tmp;
2916 abs_min = abs_tmp;
2917 }
2918 }
2919
2920 if (octave::math::isnan (tmp_min))
2921 {
2922 result.elem (i) = Complex_NaN_result;
2923 idx_arg.elem (i) = 0;
2924 }
2925 else
2926 {
2927 result.elem (i) = tmp_min;
2928 idx_arg.elem (i) = idx_j;
2929 }
2930 }
2931 }
2932
2933 return result;
2934 }
2935
2936 ComplexColumnVector
row_max(void) const2937 ComplexMatrix::row_max (void) const
2938 {
2939 Array<octave_idx_type> dummy_idx;
2940 return row_max (dummy_idx);
2941 }
2942
2943 ComplexColumnVector
row_max(Array<octave_idx_type> & idx_arg) const2944 ComplexMatrix::row_max (Array<octave_idx_type>& idx_arg) const
2945 {
2946 ComplexColumnVector result;
2947
2948 octave_idx_type nr = rows ();
2949 octave_idx_type nc = cols ();
2950
2951 if (nr > 0 && nc > 0)
2952 {
2953 result.resize (nr);
2954 idx_arg.resize (dim_vector (nr, 1));
2955
2956 for (octave_idx_type i = 0; i < nr; i++)
2957 {
2958 bool real_only = row_is_real_only (i);
2959
2960 octave_idx_type idx_j;
2961
2962 Complex tmp_max;
2963
2964 double abs_max = octave::numeric_limits<double>::NaN ();
2965
2966 for (idx_j = 0; idx_j < nc; idx_j++)
2967 {
2968 tmp_max = elem (i, idx_j);
2969
2970 if (! octave::math::isnan (tmp_max))
2971 {
2972 abs_max = (real_only ? tmp_max.real ()
2973 : std::abs (tmp_max));
2974 break;
2975 }
2976 }
2977
2978 for (octave_idx_type j = idx_j+1; j < nc; j++)
2979 {
2980 Complex tmp = elem (i, j);
2981
2982 if (octave::math::isnan (tmp))
2983 continue;
2984
2985 double abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
2986
2987 if (abs_tmp > abs_max)
2988 {
2989 idx_j = j;
2990 tmp_max = tmp;
2991 abs_max = abs_tmp;
2992 }
2993 }
2994
2995 if (octave::math::isnan (tmp_max))
2996 {
2997 result.elem (i) = Complex_NaN_result;
2998 idx_arg.elem (i) = 0;
2999 }
3000 else
3001 {
3002 result.elem (i) = tmp_max;
3003 idx_arg.elem (i) = idx_j;
3004 }
3005 }
3006 }
3007
3008 return result;
3009 }
3010
3011 ComplexRowVector
column_min(void) const3012 ComplexMatrix::column_min (void) const
3013 {
3014 Array<octave_idx_type> dummy_idx;
3015 return column_min (dummy_idx);
3016 }
3017
3018 ComplexRowVector
column_min(Array<octave_idx_type> & idx_arg) const3019 ComplexMatrix::column_min (Array<octave_idx_type>& idx_arg) const
3020 {
3021 ComplexRowVector result;
3022
3023 octave_idx_type nr = rows ();
3024 octave_idx_type nc = cols ();
3025
3026 if (nr > 0 && nc > 0)
3027 {
3028 result.resize (nc);
3029 idx_arg.resize (dim_vector (1, nc));
3030
3031 for (octave_idx_type j = 0; j < nc; j++)
3032 {
3033 bool real_only = column_is_real_only (j);
3034
3035 octave_idx_type idx_i;
3036
3037 Complex tmp_min;
3038
3039 double abs_min = octave::numeric_limits<double>::NaN ();
3040
3041 for (idx_i = 0; idx_i < nr; idx_i++)
3042 {
3043 tmp_min = elem (idx_i, j);
3044
3045 if (! octave::math::isnan (tmp_min))
3046 {
3047 abs_min = (real_only ? tmp_min.real ()
3048 : std::abs (tmp_min));
3049 break;
3050 }
3051 }
3052
3053 for (octave_idx_type i = idx_i+1; i < nr; i++)
3054 {
3055 Complex tmp = elem (i, j);
3056
3057 if (octave::math::isnan (tmp))
3058 continue;
3059
3060 double abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
3061
3062 if (abs_tmp < abs_min)
3063 {
3064 idx_i = i;
3065 tmp_min = tmp;
3066 abs_min = abs_tmp;
3067 }
3068 }
3069
3070 if (octave::math::isnan (tmp_min))
3071 {
3072 result.elem (j) = Complex_NaN_result;
3073 idx_arg.elem (j) = 0;
3074 }
3075 else
3076 {
3077 result.elem (j) = tmp_min;
3078 idx_arg.elem (j) = idx_i;
3079 }
3080 }
3081 }
3082
3083 return result;
3084 }
3085
3086 ComplexRowVector
column_max(void) const3087 ComplexMatrix::column_max (void) const
3088 {
3089 Array<octave_idx_type> dummy_idx;
3090 return column_max (dummy_idx);
3091 }
3092
3093 ComplexRowVector
column_max(Array<octave_idx_type> & idx_arg) const3094 ComplexMatrix::column_max (Array<octave_idx_type>& idx_arg) const
3095 {
3096 ComplexRowVector result;
3097
3098 octave_idx_type nr = rows ();
3099 octave_idx_type nc = cols ();
3100
3101 if (nr > 0 && nc > 0)
3102 {
3103 result.resize (nc);
3104 idx_arg.resize (dim_vector (1, nc));
3105
3106 for (octave_idx_type j = 0; j < nc; j++)
3107 {
3108 bool real_only = column_is_real_only (j);
3109
3110 octave_idx_type idx_i;
3111
3112 Complex tmp_max;
3113
3114 double abs_max = octave::numeric_limits<double>::NaN ();
3115
3116 for (idx_i = 0; idx_i < nr; idx_i++)
3117 {
3118 tmp_max = elem (idx_i, j);
3119
3120 if (! octave::math::isnan (tmp_max))
3121 {
3122 abs_max = (real_only ? tmp_max.real ()
3123 : std::abs (tmp_max));
3124 break;
3125 }
3126 }
3127
3128 for (octave_idx_type i = idx_i+1; i < nr; i++)
3129 {
3130 Complex tmp = elem (i, j);
3131
3132 if (octave::math::isnan (tmp))
3133 continue;
3134
3135 double abs_tmp = (real_only ? tmp.real () : std::abs (tmp));
3136
3137 if (abs_tmp > abs_max)
3138 {
3139 idx_i = i;
3140 tmp_max = tmp;
3141 abs_max = abs_tmp;
3142 }
3143 }
3144
3145 if (octave::math::isnan (tmp_max))
3146 {
3147 result.elem (j) = Complex_NaN_result;
3148 idx_arg.elem (j) = 0;
3149 }
3150 else
3151 {
3152 result.elem (j) = tmp_max;
3153 idx_arg.elem (j) = idx_i;
3154 }
3155 }
3156 }
3157
3158 return result;
3159 }
3160
3161 // i/o
3162
3163 std::ostream&
operator <<(std::ostream & os,const ComplexMatrix & a)3164 operator << (std::ostream& os, const ComplexMatrix& a)
3165 {
3166 for (octave_idx_type i = 0; i < a.rows (); i++)
3167 {
3168 for (octave_idx_type j = 0; j < a.cols (); j++)
3169 {
3170 os << ' ';
3171 octave_write_complex (os, a.elem (i, j));
3172 }
3173 os << "\n";
3174 }
3175 return os;
3176 }
3177
3178 std::istream&
operator >>(std::istream & is,ComplexMatrix & a)3179 operator >> (std::istream& is, ComplexMatrix& a)
3180 {
3181 octave_idx_type nr = a.rows ();
3182 octave_idx_type nc = a.cols ();
3183
3184 if (nr > 0 && nc > 0)
3185 {
3186 Complex tmp;
3187 for (octave_idx_type i = 0; i < nr; i++)
3188 for (octave_idx_type j = 0; j < nc; j++)
3189 {
3190 tmp = octave_read_value<Complex> (is);
3191 if (is)
3192 a.elem (i, j) = tmp;
3193 else
3194 return is;
3195 }
3196 }
3197
3198 return is;
3199 }
3200
3201 ComplexMatrix
Givens(const Complex & x,const Complex & y)3202 Givens (const Complex& x, const Complex& y)
3203 {
3204 double cc;
3205 Complex cs, temp_r;
3206
3207 F77_FUNC (zlartg, ZLARTG) (F77_CONST_DBLE_CMPLX_ARG (&x),
3208 F77_CONST_DBLE_CMPLX_ARG (&y),
3209 cc,
3210 F77_DBLE_CMPLX_ARG (&cs),
3211 F77_DBLE_CMPLX_ARG (&temp_r));
3212
3213 ComplexMatrix g (2, 2);
3214
3215 g.elem (0, 0) = cc;
3216 g.elem (1, 1) = cc;
3217 g.elem (0, 1) = cs;
3218 g.elem (1, 0) = -conj (cs);
3219
3220 return g;
3221 }
3222
3223 ComplexMatrix
Sylvester(const ComplexMatrix & a,const ComplexMatrix & b,const ComplexMatrix & c)3224 Sylvester (const ComplexMatrix& a, const ComplexMatrix& b,
3225 const ComplexMatrix& c)
3226 {
3227 ComplexMatrix retval;
3228
3229 // FIXME: need to check that a, b, and c are all the same size.
3230
3231 // Compute Schur decompositions
3232
3233 octave::math::schur<ComplexMatrix> as (a, "U");
3234 octave::math::schur<ComplexMatrix> bs (b, "U");
3235
3236 // Transform c to new coordinates.
3237
3238 ComplexMatrix ua = as.unitary_matrix ();
3239 ComplexMatrix sch_a = as.schur_matrix ();
3240
3241 ComplexMatrix ub = bs.unitary_matrix ();
3242 ComplexMatrix sch_b = bs.schur_matrix ();
3243
3244 ComplexMatrix cx = ua.hermitian () * c * ub;
3245
3246 // Solve the sylvester equation, back-transform, and return the solution.
3247
3248 F77_INT a_nr = octave::to_f77_int (a.rows ());
3249 F77_INT b_nr = octave::to_f77_int (b.rows ());
3250
3251 double scale;
3252 F77_INT info;
3253
3254 Complex *pa = sch_a.fortran_vec ();
3255 Complex *pb = sch_b.fortran_vec ();
3256 Complex *px = cx.fortran_vec ();
3257
3258 F77_XFCN (ztrsyl, ZTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
3259 F77_CONST_CHAR_ARG2 ("N", 1),
3260 1, a_nr, b_nr, F77_DBLE_CMPLX_ARG (pa), a_nr, F77_DBLE_CMPLX_ARG (pb),
3261 b_nr, F77_DBLE_CMPLX_ARG (px), a_nr, scale, info
3262 F77_CHAR_ARG_LEN (1)
3263 F77_CHAR_ARG_LEN (1)));
3264
3265 // FIXME: check info?
3266
3267 retval = ua * cx * ub.hermitian ();
3268
3269 return retval;
3270 }
3271
3272 ComplexMatrix
operator *(const ComplexMatrix & m,const Matrix & a)3273 operator * (const ComplexMatrix& m, const Matrix& a)
3274 {
3275 if (m.columns () > std::min (m.rows (), a.columns ()) / 10)
3276 return ComplexMatrix (real (m) * a, imag (m) * a);
3277 else
3278 return m * ComplexMatrix (a);
3279 }
3280
3281 ComplexMatrix
operator *(const Matrix & m,const ComplexMatrix & a)3282 operator * (const Matrix& m, const ComplexMatrix& a)
3283 {
3284 if (a.rows () > std::min (m.rows (), a.columns ()) / 10)
3285 return ComplexMatrix (m * real (a), m * imag (a));
3286 else
3287 return ComplexMatrix (m) * a;
3288 }
3289
3290 /*
3291
3292 ## Simple Dot Product, Matrix-Vector, and Matrix-Matrix Unit tests
3293 %!assert ([1+i 2+i 3+i] * [ 4+i ; 5+i ; 6+i], 29+21i, 1e-14)
3294 %!assert ([1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i], [15 + 14i ; 37 + 18i], 1e-14)
3295 %!assert ([1+i 2+i ; 3+i 4+i ] * [5+i 6+i ; 7+i 8+i], [17 + 15i 20 + 17i; 41 + 19i 48 + 21i], 1e-14)
3296 %!assert ([1 i]*[i 0]', -i)
3297
3298 ## Test some simple identities
3299 %!shared M, cv, rv
3300 %! M = randn (10,10) + i*rand (10,10);
3301 %! cv = randn (10,1) + i*rand (10,1);
3302 %! rv = randn (1,10) + i*rand (1,10);
3303 %!assert ([M*cv,M*cv], M*[cv,cv], 1e-14)
3304 %!assert ([M.'*cv,M.'*cv], M.'*[cv,cv], 1e-14)
3305 %!assert ([M'*cv,M'*cv], M'*[cv,cv], 1e-14)
3306 %!assert ([rv*M;rv*M], [rv;rv]*M, 1e-14)
3307 %!assert ([rv*M.';rv*M.'], [rv;rv]*M.', 1e-14)
3308 %!assert ([rv*M';rv*M'], [rv;rv]*M', 1e-14)
3309 %!assert (2*rv*cv, [rv,rv]*[cv;cv], 2e-14)
3310
3311 */
3312
3313 static inline char
get_blas_trans_arg(bool trans,bool conj)3314 get_blas_trans_arg (bool trans, bool conj)
3315 {
3316 return trans ? (conj ? 'C' : 'T') : 'N';
3317 }
3318
3319 // the general GEMM operation
3320
3321 ComplexMatrix
xgemm(const ComplexMatrix & a,const ComplexMatrix & b,blas_trans_type transa,blas_trans_type transb)3322 xgemm (const ComplexMatrix& a, const ComplexMatrix& b,
3323 blas_trans_type transa, blas_trans_type transb)
3324 {
3325 ComplexMatrix retval;
3326
3327 bool tra = transa != blas_no_trans;
3328 bool trb = transb != blas_no_trans;
3329 bool cja = transa == blas_conj_trans;
3330 bool cjb = transb == blas_conj_trans;
3331
3332 F77_INT a_nr = octave::to_f77_int (tra ? a.cols () : a.rows ());
3333 F77_INT a_nc = octave::to_f77_int (tra ? a.rows () : a.cols ());
3334
3335 F77_INT b_nr = octave::to_f77_int (trb ? b.cols () : b.rows ());
3336 F77_INT b_nc = octave::to_f77_int (trb ? b.rows () : b.cols ());
3337
3338 if (a_nc != b_nr)
3339 octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
3340
3341 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
3342 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3343 else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
3344 {
3345 F77_INT lda = octave::to_f77_int (a.rows ());
3346
3347 // FIXME: looking at the reference BLAS, it appears that it
3348 // should not be necessary to initialize the output matrix if
3349 // BETA is 0 in the call to ZHERK, but ATLAS appears to
3350 // use the result matrix before zeroing the elements.
3351
3352 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3353 Complex *c = retval.fortran_vec ();
3354
3355 const char ctra = get_blas_trans_arg (tra, cja);
3356 if (cja || cjb)
3357 {
3358 F77_XFCN (zherk, ZHERK, (F77_CONST_CHAR_ARG2 ("U", 1),
3359 F77_CONST_CHAR_ARG2 (&ctra, 1),
3360 a_nr, a_nc, 1.0,
3361 F77_CONST_DBLE_CMPLX_ARG (a.data ()), lda, 0.0, F77_DBLE_CMPLX_ARG (c), a_nr
3362 F77_CHAR_ARG_LEN (1)
3363 F77_CHAR_ARG_LEN (1)));
3364 for (F77_INT j = 0; j < a_nr; j++)
3365 for (F77_INT i = 0; i < j; i++)
3366 retval.xelem (j,i) = octave::math::conj (retval.xelem (i,j));
3367 }
3368 else
3369 {
3370 F77_XFCN (zsyrk, ZSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),
3371 F77_CONST_CHAR_ARG2 (&ctra, 1),
3372 a_nr, a_nc, 1.0,
3373 F77_CONST_DBLE_CMPLX_ARG (a.data ()), lda, 0.0, F77_DBLE_CMPLX_ARG (c), a_nr
3374 F77_CHAR_ARG_LEN (1)
3375 F77_CHAR_ARG_LEN (1)));
3376 for (F77_INT j = 0; j < a_nr; j++)
3377 for (F77_INT i = 0; i < j; i++)
3378 retval.xelem (j,i) = retval.xelem (i,j);
3379
3380 }
3381
3382 }
3383 else
3384 {
3385 F77_INT lda = octave::to_f77_int (a.rows ());
3386 F77_INT tda = octave::to_f77_int (a.cols ());
3387 F77_INT ldb = octave::to_f77_int (b.rows ());
3388 F77_INT tdb = octave::to_f77_int (b.cols ());
3389
3390 retval = ComplexMatrix (a_nr, b_nc, 0.0);
3391 Complex *c = retval.fortran_vec ();
3392
3393 if (b_nc == 1 && a_nr == 1)
3394 {
3395 if (cja == cjb)
3396 {
3397 F77_FUNC (xzdotu, XZDOTU) (a_nc, F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1,
3398 F77_CONST_DBLE_CMPLX_ARG (b.data ()), 1,
3399 F77_DBLE_CMPLX_ARG (c));
3400 if (cja) *c = octave::math::conj (*c);
3401 }
3402 else if (cja)
3403 F77_FUNC (xzdotc, XZDOTC) (a_nc, F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1,
3404 F77_CONST_DBLE_CMPLX_ARG (b.data ()), 1,
3405 F77_DBLE_CMPLX_ARG (c));
3406 else
3407 F77_FUNC (xzdotc, XZDOTC) (a_nc, F77_CONST_DBLE_CMPLX_ARG (b.data ()), 1,
3408 F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1,
3409 F77_DBLE_CMPLX_ARG (c));
3410 }
3411 else if (b_nc == 1 && ! cjb)
3412 {
3413 const char ctra = get_blas_trans_arg (tra, cja);
3414 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),
3415 lda, tda, 1.0, F77_CONST_DBLE_CMPLX_ARG (a.data ()), lda,
3416 F77_CONST_DBLE_CMPLX_ARG (b.data ()), 1, 0.0, F77_DBLE_CMPLX_ARG (c), 1
3417 F77_CHAR_ARG_LEN (1)));
3418 }
3419 else if (a_nr == 1 && ! cja && ! cjb)
3420 {
3421 const char crevtrb = get_blas_trans_arg (! trb, cjb);
3422 F77_XFCN (zgemv, ZGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),
3423 ldb, tdb, 1.0, F77_CONST_DBLE_CMPLX_ARG (b.data ()), ldb,
3424 F77_CONST_DBLE_CMPLX_ARG (a.data ()), 1, 0.0, F77_DBLE_CMPLX_ARG (c), 1
3425 F77_CHAR_ARG_LEN (1)));
3426 }
3427 else
3428 {
3429 const char ctra = get_blas_trans_arg (tra, cja);
3430 const char ctrb = get_blas_trans_arg (trb, cjb);
3431 F77_XFCN (zgemm, ZGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),
3432 F77_CONST_CHAR_ARG2 (&ctrb, 1),
3433 a_nr, b_nc, a_nc, 1.0, F77_CONST_DBLE_CMPLX_ARG (a.data ()),
3434 lda, F77_CONST_DBLE_CMPLX_ARG (b.data ()), ldb, 0.0, F77_DBLE_CMPLX_ARG (c),
3435 a_nr
3436 F77_CHAR_ARG_LEN (1)
3437 F77_CHAR_ARG_LEN (1)));
3438 }
3439 }
3440
3441 return retval;
3442 }
3443
3444 ComplexMatrix
operator *(const ComplexMatrix & a,const ComplexMatrix & b)3445 operator * (const ComplexMatrix& a, const ComplexMatrix& b)
3446 {
3447 return xgemm (a, b);
3448 }
3449
3450 // FIXME: it would be nice to share code among the min/max functions below.
3451
3452 #define EMPTY_RETURN_CHECK(T) \
3453 if (nr == 0 || nc == 0) \
3454 return T (nr, nc);
3455
3456 ComplexMatrix
min(const Complex & c,const ComplexMatrix & m)3457 min (const Complex& c, const ComplexMatrix& m)
3458 {
3459 octave_idx_type nr = m.rows ();
3460 octave_idx_type nc = m.columns ();
3461
3462 EMPTY_RETURN_CHECK (ComplexMatrix);
3463
3464 ComplexMatrix result (nr, nc);
3465
3466 for (octave_idx_type j = 0; j < nc; j++)
3467 for (octave_idx_type i = 0; i < nr; i++)
3468 {
3469 octave_quit ();
3470 result(i, j) = octave::math::min (c, m(i, j));
3471 }
3472
3473 return result;
3474 }
3475
3476 ComplexMatrix
min(const ComplexMatrix & m,const Complex & c)3477 min (const ComplexMatrix& m, const Complex& c)
3478 {
3479 return min (c, m);
3480 }
3481
3482 ComplexMatrix
min(const ComplexMatrix & a,const ComplexMatrix & b)3483 min (const ComplexMatrix& a, const ComplexMatrix& b)
3484 {
3485 octave_idx_type nr = a.rows ();
3486 octave_idx_type nc = a.columns ();
3487
3488 if (nr != b.rows () || nc != b.columns ())
3489 (*current_liboctave_error_handler)
3490 ("two-arg min requires same size arguments");
3491
3492 EMPTY_RETURN_CHECK (ComplexMatrix);
3493
3494 ComplexMatrix result (nr, nc);
3495
3496 for (octave_idx_type j = 0; j < nc; j++)
3497 {
3498 bool columns_are_real_only = true;
3499 for (octave_idx_type i = 0; i < nr; i++)
3500 {
3501 octave_quit ();
3502 if (std::imag (a(i, j)) != 0.0 || std::imag (b(i, j)) != 0.0)
3503 {
3504 columns_are_real_only = false;
3505 break;
3506 }
3507 }
3508
3509 if (columns_are_real_only)
3510 {
3511 for (octave_idx_type i = 0; i < nr; i++)
3512 result(i, j) = octave::math::min (std::real (a(i, j)),
3513 std::real (b(i, j)));
3514 }
3515 else
3516 {
3517 for (octave_idx_type i = 0; i < nr; i++)
3518 {
3519 octave_quit ();
3520 result(i, j) = octave::math::min (a(i, j), b(i, j));
3521 }
3522 }
3523 }
3524
3525 return result;
3526 }
3527
3528 ComplexMatrix
max(const Complex & c,const ComplexMatrix & m)3529 max (const Complex& c, const ComplexMatrix& m)
3530 {
3531 octave_idx_type nr = m.rows ();
3532 octave_idx_type nc = m.columns ();
3533
3534 EMPTY_RETURN_CHECK (ComplexMatrix);
3535
3536 ComplexMatrix result (nr, nc);
3537
3538 for (octave_idx_type j = 0; j < nc; j++)
3539 for (octave_idx_type i = 0; i < nr; i++)
3540 {
3541 octave_quit ();
3542 result(i, j) = octave::math::max (c, m(i, j));
3543 }
3544
3545 return result;
3546 }
3547
3548 ComplexMatrix
max(const ComplexMatrix & m,const Complex & c)3549 max (const ComplexMatrix& m, const Complex& c)
3550 {
3551 return max (c, m);
3552 }
3553
3554 ComplexMatrix
max(const ComplexMatrix & a,const ComplexMatrix & b)3555 max (const ComplexMatrix& a, const ComplexMatrix& b)
3556 {
3557 octave_idx_type nr = a.rows ();
3558 octave_idx_type nc = a.columns ();
3559
3560 if (nr != b.rows () || nc != b.columns ())
3561 (*current_liboctave_error_handler)
3562 ("two-arg max requires same size arguments");
3563
3564 EMPTY_RETURN_CHECK (ComplexMatrix);
3565
3566 ComplexMatrix result (nr, nc);
3567
3568 for (octave_idx_type j = 0; j < nc; j++)
3569 {
3570 bool columns_are_real_only = true;
3571 for (octave_idx_type i = 0; i < nr; i++)
3572 {
3573 octave_quit ();
3574 if (std::imag (a(i, j)) != 0.0 || std::imag (b(i, j)) != 0.0)
3575 {
3576 columns_are_real_only = false;
3577 break;
3578 }
3579 }
3580
3581 // FIXME: is it so much faster?
3582 if (columns_are_real_only)
3583 {
3584 for (octave_idx_type i = 0; i < nr; i++)
3585 {
3586 octave_quit ();
3587 result(i, j) = octave::math::max (std::real (a(i, j)),
3588 std::real (b(i, j)));
3589 }
3590 }
3591 else
3592 {
3593 for (octave_idx_type i = 0; i < nr; i++)
3594 {
3595 octave_quit ();
3596 result(i, j) = octave::math::max (a(i, j), b(i, j));
3597 }
3598 }
3599 }
3600
3601 return result;
3602 }
3603
linspace(const ComplexColumnVector & x1,const ComplexColumnVector & x2,octave_idx_type n)3604 ComplexMatrix linspace (const ComplexColumnVector& x1,
3605 const ComplexColumnVector& x2,
3606 octave_idx_type n)
3607 {
3608 octave_idx_type m = x1.numel ();
3609
3610 if (x2.numel () != m)
3611 (*current_liboctave_error_handler)
3612 ("linspace: vectors must be of equal length");
3613
3614 ComplexMatrix retval;
3615
3616 if (n < 1)
3617 {
3618 retval.clear (m, 0);
3619 return retval;
3620 }
3621
3622 retval.clear (m, n);
3623 for (octave_idx_type i = 0; i < m; i++)
3624 retval.xelem (i, 0) = x1(i);
3625
3626 // The last column is unused so temporarily store delta there
3627 Complex *delta = &retval.xelem (0, n-1);
3628 for (octave_idx_type i = 0; i < m; i++)
3629 delta[i] = (x1(i) == x2(i)) ? 0 : (x2(i) - x1(i)) / (n - 1.0);
3630
3631 for (octave_idx_type j = 1; j < n-1; j++)
3632 for (octave_idx_type i = 0; i < m; i++)
3633 retval.xelem (i, j) = x1(i) + static_cast<double> (j)*delta[i];
3634
3635 for (octave_idx_type i = 0; i < m; i++)
3636 retval.xelem (i, n-1) = x2(i);
3637
3638 return retval;
3639 }
3640
3641 MS_CMP_OPS (ComplexMatrix, Complex)
3642 MS_BOOL_OPS (ComplexMatrix, Complex)
3643
3644 SM_CMP_OPS (Complex, ComplexMatrix)
3645 SM_BOOL_OPS (Complex, ComplexMatrix)
3646
3647 MM_CMP_OPS (ComplexMatrix, ComplexMatrix)
3648 MM_BOOL_OPS (ComplexMatrix, ComplexMatrix)
3649