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