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 <istream>
32 #include <limits>
33 #include <ostream>
34
35 #include "Array-util.h"
36 #include "CColVector.h"
37 #include "CMatrix.h"
38 #include "DET.h"
39 #include "PermMatrix.h"
40 #include "boolMatrix.h"
41 #include "byte-swap.h"
42 #include "chMatrix.h"
43 #include "chol.h"
44 #include "dColVector.h"
45 #include "dDiagMatrix.h"
46 #include "dMatrix.h"
47 #include "dRowVector.h"
48 #include "lo-blas-proto.h"
49 #include "lo-error.h"
50 #include "lo-ieee.h"
51 #include "lo-lapack-proto.h"
52 #include "lo-mappers.h"
53 #include "lo-utils.h"
54 #include "mx-dm-m.h"
55 #include "mx-inlines.cc"
56 #include "mx-m-dm.h"
57 #include "mx-op-defs.h"
58 #include "oct-cmplx.h"
59 #include "oct-fftw.h"
60 #include "oct-locbuf.h"
61 #include "oct-norm.h"
62 #include "quit.h"
63 #include "schur.h"
64 #include "svd.h"
65
66 // Matrix class.
67
Matrix(const RowVector & rv)68 Matrix::Matrix (const RowVector& rv)
69 : NDArray (rv)
70 { }
71
Matrix(const ColumnVector & cv)72 Matrix::Matrix (const ColumnVector& cv)
73 : NDArray (cv)
74 { }
75
Matrix(const DiagMatrix & a)76 Matrix::Matrix (const DiagMatrix& a)
77 : NDArray (a.dims (), 0.0)
78 {
79 for (octave_idx_type i = 0; i < a.length (); i++)
80 elem (i, i) = a.elem (i, i);
81 }
82
Matrix(const MDiagArray2<double> & a)83 Matrix::Matrix (const MDiagArray2<double>& a)
84 : NDArray (a.dims (), 0.0)
85 {
86 for (octave_idx_type i = 0; i < a.length (); i++)
87 elem (i, i) = a.elem (i, i);
88 }
89
Matrix(const DiagArray2<double> & a)90 Matrix::Matrix (const DiagArray2<double>& a)
91 : NDArray (a.dims (), 0.0)
92 {
93 for (octave_idx_type i = 0; i < a.length (); i++)
94 elem (i, i) = a.elem (i, i);
95 }
96
Matrix(const PermMatrix & a)97 Matrix::Matrix (const PermMatrix& a)
98 : NDArray (a.dims (), 0.0)
99 {
100 const Array<octave_idx_type> ia (a.col_perm_vec ());
101 octave_idx_type len = a.rows ();
102 for (octave_idx_type i = 0; i < len; i++)
103 elem (ia(i), i) = 1.0;
104 }
105
106 // FIXME: could we use a templated mixed-type copy function here?
107
Matrix(const boolMatrix & a)108 Matrix::Matrix (const boolMatrix& a)
109 : NDArray (a)
110 { }
111
Matrix(const charMatrix & a)112 Matrix::Matrix (const charMatrix& a)
113 : NDArray (a.dims ())
114 {
115 for (octave_idx_type i = 0; i < a.rows (); i++)
116 for (octave_idx_type j = 0; j < a.cols (); j++)
117 elem (i, j) = static_cast<unsigned char> (a.elem (i, j));
118 }
119
120 bool
operator ==(const Matrix & a) const121 Matrix::operator == (const Matrix& a) const
122 {
123 if (rows () != a.rows () || cols () != a.cols ())
124 return false;
125
126 return mx_inline_equal (numel (), data (), a.data ());
127 }
128
129 bool
operator !=(const Matrix & a) const130 Matrix::operator != (const Matrix& a) const
131 {
132 return !(*this == a);
133 }
134
135 bool
issymmetric(void) const136 Matrix::issymmetric (void) const
137 {
138 if (issquare () && rows () > 0)
139 {
140 for (octave_idx_type i = 0; i < rows (); i++)
141 for (octave_idx_type j = i+1; j < cols (); j++)
142 if (elem (i, j) != elem (j, i))
143 return false;
144
145 return true;
146 }
147
148 return false;
149 }
150
151 Matrix&
insert(const Matrix & a,octave_idx_type r,octave_idx_type c)152 Matrix::insert (const Matrix& a, octave_idx_type r, octave_idx_type c)
153 {
154 Array<double>::insert (a, r, c);
155 return *this;
156 }
157
158 Matrix&
insert(const RowVector & a,octave_idx_type r,octave_idx_type c)159 Matrix::insert (const RowVector& a, octave_idx_type r, octave_idx_type c)
160 {
161 octave_idx_type a_len = a.numel ();
162
163 if (r < 0 || r >= rows () || c < 0 || c + a_len > cols ())
164 (*current_liboctave_error_handler) ("range error for insert");
165
166 if (a_len > 0)
167 {
168 make_unique ();
169
170 for (octave_idx_type i = 0; i < a_len; i++)
171 xelem (r, c+i) = a.elem (i);
172 }
173
174 return *this;
175 }
176
177 Matrix&
insert(const ColumnVector & a,octave_idx_type r,octave_idx_type c)178 Matrix::insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c)
179 {
180 octave_idx_type a_len = a.numel ();
181
182 if (r < 0 || r + a_len > rows () || c < 0 || c >= cols ())
183 (*current_liboctave_error_handler) ("range error for insert");
184
185 if (a_len > 0)
186 {
187 make_unique ();
188
189 for (octave_idx_type i = 0; i < a_len; i++)
190 xelem (r+i, c) = a.elem (i);
191 }
192
193 return *this;
194 }
195
196 Matrix&
insert(const DiagMatrix & a,octave_idx_type r,octave_idx_type c)197 Matrix::insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c)
198 {
199 octave_idx_type a_nr = a.rows ();
200 octave_idx_type a_nc = a.cols ();
201
202 if (r < 0 || r + a_nr > rows () || c < 0 || c + a_nc > cols ())
203 (*current_liboctave_error_handler) ("range error for insert");
204
205 fill (0.0, r, c, r + a_nr - 1, c + a_nc - 1);
206
207 octave_idx_type a_len = a.length ();
208
209 if (a_len > 0)
210 {
211 make_unique ();
212
213 for (octave_idx_type i = 0; i < a_len; i++)
214 xelem (r+i, c+i) = a.elem (i, i);
215 }
216
217 return *this;
218 }
219
220 Matrix&
fill(double val)221 Matrix::fill (double val)
222 {
223 octave_idx_type nr = rows ();
224 octave_idx_type nc = cols ();
225
226 if (nr > 0 && nc > 0)
227 {
228 make_unique ();
229
230 for (octave_idx_type j = 0; j < nc; j++)
231 for (octave_idx_type i = 0; i < nr; i++)
232 xelem (i, j) = val;
233 }
234
235 return *this;
236 }
237
238 Matrix&
fill(double val,octave_idx_type r1,octave_idx_type c1,octave_idx_type r2,octave_idx_type c2)239 Matrix::fill (double val, octave_idx_type r1, octave_idx_type c1,
240 octave_idx_type r2, octave_idx_type c2)
241 {
242 octave_idx_type nr = rows ();
243 octave_idx_type nc = cols ();
244
245 if (r1 < 0 || r2 < 0 || c1 < 0 || c2 < 0
246 || r1 >= nr || r2 >= nr || c1 >= nc || c2 >= nc)
247 (*current_liboctave_error_handler) ("range error for fill");
248
249 if (r1 > r2) { std::swap (r1, r2); }
250 if (c1 > c2) { std::swap (c1, c2); }
251
252 if (r2 >= r1 && c2 >= c1)
253 {
254 make_unique ();
255
256 for (octave_idx_type j = c1; j <= c2; j++)
257 for (octave_idx_type i = r1; i <= r2; i++)
258 xelem (i, j) = val;
259 }
260
261 return *this;
262 }
263
264 Matrix
append(const Matrix & a) const265 Matrix::append (const Matrix& a) const
266 {
267 octave_idx_type nr = rows ();
268 octave_idx_type nc = cols ();
269 if (nr != a.rows ())
270 (*current_liboctave_error_handler) ("row dimension mismatch for append");
271
272 octave_idx_type nc_insert = nc;
273 Matrix retval (nr, nc + a.cols ());
274 retval.insert (*this, 0, 0);
275 retval.insert (a, 0, nc_insert);
276 return retval;
277 }
278
279 Matrix
append(const RowVector & a) const280 Matrix::append (const RowVector& a) const
281 {
282 octave_idx_type nr = rows ();
283 octave_idx_type nc = cols ();
284 if (nr != 1)
285 (*current_liboctave_error_handler) ("row dimension mismatch for append");
286
287 octave_idx_type nc_insert = nc;
288 Matrix retval (nr, nc + a.numel ());
289 retval.insert (*this, 0, 0);
290 retval.insert (a, 0, nc_insert);
291 return retval;
292 }
293
294 Matrix
append(const ColumnVector & a) const295 Matrix::append (const ColumnVector& a) const
296 {
297 octave_idx_type nr = rows ();
298 octave_idx_type nc = cols ();
299 if (nr != a.numel ())
300 (*current_liboctave_error_handler) ("row dimension mismatch for append");
301
302 octave_idx_type nc_insert = nc;
303 Matrix retval (nr, nc + 1);
304 retval.insert (*this, 0, 0);
305 retval.insert (a, 0, nc_insert);
306 return retval;
307 }
308
309 Matrix
append(const DiagMatrix & a) const310 Matrix::append (const DiagMatrix& a) const
311 {
312 octave_idx_type nr = rows ();
313 octave_idx_type nc = cols ();
314 if (nr != a.rows ())
315 (*current_liboctave_error_handler) ("row dimension mismatch for append");
316
317 octave_idx_type nc_insert = nc;
318 Matrix retval (nr, nc + a.cols ());
319 retval.insert (*this, 0, 0);
320 retval.insert (a, 0, nc_insert);
321 return retval;
322 }
323
324 Matrix
stack(const Matrix & a) const325 Matrix::stack (const Matrix& a) const
326 {
327 octave_idx_type nr = rows ();
328 octave_idx_type nc = cols ();
329 if (nc != a.cols ())
330 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
331
332 octave_idx_type nr_insert = nr;
333 Matrix retval (nr + a.rows (), nc);
334 retval.insert (*this, 0, 0);
335 retval.insert (a, nr_insert, 0);
336 return retval;
337 }
338
339 Matrix
stack(const RowVector & a) const340 Matrix::stack (const RowVector& a) const
341 {
342 octave_idx_type nr = rows ();
343 octave_idx_type nc = cols ();
344 if (nc != a.numel ())
345 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
346
347 octave_idx_type nr_insert = nr;
348 Matrix retval (nr + 1, nc);
349 retval.insert (*this, 0, 0);
350 retval.insert (a, nr_insert, 0);
351 return retval;
352 }
353
354 Matrix
stack(const ColumnVector & a) const355 Matrix::stack (const ColumnVector& a) const
356 {
357 octave_idx_type nr = rows ();
358 octave_idx_type nc = cols ();
359 if (nc != 1)
360 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
361
362 octave_idx_type nr_insert = nr;
363 Matrix retval (nr + a.numel (), nc);
364 retval.insert (*this, 0, 0);
365 retval.insert (a, nr_insert, 0);
366 return retval;
367 }
368
369 Matrix
stack(const DiagMatrix & a) const370 Matrix::stack (const DiagMatrix& a) const
371 {
372 octave_idx_type nr = rows ();
373 octave_idx_type nc = cols ();
374 if (nc != a.cols ())
375 (*current_liboctave_error_handler) ("column dimension mismatch for stack");
376
377 octave_idx_type nr_insert = nr;
378 Matrix retval (nr + a.rows (), nc);
379 retval.insert (*this, 0, 0);
380 retval.insert (a, nr_insert, 0);
381 return retval;
382 }
383
384 Matrix
real(const ComplexMatrix & a)385 real (const ComplexMatrix& a)
386 {
387 return do_mx_unary_op<double, Complex> (a, mx_inline_real);
388 }
389
390 Matrix
imag(const ComplexMatrix & a)391 imag (const ComplexMatrix& a)
392 {
393 return do_mx_unary_op<double, Complex> (a, mx_inline_imag);
394 }
395
396 Matrix
extract(octave_idx_type r1,octave_idx_type c1,octave_idx_type r2,octave_idx_type c2) const397 Matrix::extract (octave_idx_type r1, octave_idx_type c1,
398 octave_idx_type r2, octave_idx_type c2) const
399 {
400 if (r1 > r2) { std::swap (r1, r2); }
401 if (c1 > c2) { std::swap (c1, c2); }
402
403 return index (idx_vector (r1, r2+1), idx_vector (c1, c2+1));
404 }
405
406 Matrix
extract_n(octave_idx_type r1,octave_idx_type c1,octave_idx_type nr,octave_idx_type nc) const407 Matrix::extract_n (octave_idx_type r1, octave_idx_type c1, octave_idx_type nr,
408 octave_idx_type nc) const
409 {
410 return index (idx_vector (r1, r1 + nr), idx_vector (c1, c1 + nc));
411 }
412
413 // extract row or column i.
414
415 RowVector
row(octave_idx_type i) const416 Matrix::row (octave_idx_type i) const
417 {
418 return index (idx_vector (i), idx_vector::colon);
419 }
420
421 ColumnVector
column(octave_idx_type i) const422 Matrix::column (octave_idx_type i) const
423 {
424 return index (idx_vector::colon, idx_vector (i));
425 }
426
427 // Local function to calculate the 1-norm.
428 static
429 double
norm1(const Matrix & a)430 norm1 (const Matrix& a)
431 {
432 double anorm = 0.0;
433 RowVector colsum = a.abs ().sum ().row (0);
434
435 for (octave_idx_type i = 0; i < colsum.numel (); i++)
436 {
437 double sum = colsum.xelem (i);
438 if (octave::math::isinf (sum) || octave::math::isnan (sum))
439 {
440 anorm = sum; // Pass Inf or NaN to output
441 break;
442 }
443 else
444 anorm = std::max (anorm, sum);
445 }
446
447 return anorm;
448 }
449
450 Matrix
inverse(void) const451 Matrix::inverse (void) const
452 {
453 octave_idx_type info;
454 double rcon;
455 MatrixType mattype (*this);
456 return inverse (mattype, info, rcon, 0, 0);
457 }
458
459 Matrix
inverse(octave_idx_type & info) const460 Matrix::inverse (octave_idx_type& info) const
461 {
462 double rcon;
463 MatrixType mattype (*this);
464 return inverse (mattype, info, rcon, 0, 0);
465 }
466
467 Matrix
inverse(octave_idx_type & info,double & rcon,bool force,bool calc_cond) const468 Matrix::inverse (octave_idx_type& info, double& rcon, bool force,
469 bool calc_cond) const
470 {
471 MatrixType mattype (*this);
472 return inverse (mattype, info, rcon, force, calc_cond);
473 }
474
475 Matrix
inverse(MatrixType & mattype) const476 Matrix::inverse (MatrixType& mattype) const
477 {
478 octave_idx_type info;
479 double rcon;
480 return inverse (mattype, info, rcon, 0, 0);
481 }
482
483 Matrix
inverse(MatrixType & mattype,octave_idx_type & info) const484 Matrix::inverse (MatrixType& mattype, octave_idx_type& info) const
485 {
486 double rcon;
487 return inverse (mattype, info, rcon, 0, 0);
488 }
489
490 Matrix
tinverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const491 Matrix::tinverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
492 bool force, bool calc_cond) const
493 {
494 Matrix retval;
495
496 F77_INT nr = octave::to_f77_int (rows ());
497 F77_INT nc = octave::to_f77_int (cols ());
498
499 if (nr != nc || nr == 0 || nc == 0)
500 (*current_liboctave_error_handler) ("inverse requires square matrix");
501
502 int typ = mattype.type ();
503 char uplo = (typ == MatrixType::Lower ? 'L' : 'U');
504 char udiag = 'N';
505 retval = *this;
506 double *tmp_data = retval.fortran_vec ();
507
508 F77_INT tmp_info = 0;
509
510 F77_XFCN (dtrtri, DTRTRI, (F77_CONST_CHAR_ARG2 (&uplo, 1),
511 F77_CONST_CHAR_ARG2 (&udiag, 1),
512 nr, tmp_data, nr, tmp_info
513 F77_CHAR_ARG_LEN (1)
514 F77_CHAR_ARG_LEN (1)));
515
516 info = tmp_info;
517
518 // Throw away extra info LAPACK gives so as to not change output.
519 rcon = 0.0;
520 if (info != 0)
521 info = -1;
522 else if (calc_cond)
523 {
524 F77_INT dtrcon_info = 0;
525 char job = '1';
526
527 OCTAVE_LOCAL_BUFFER (double, work, 3 * nr);
528 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, nr);
529
530 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&job, 1),
531 F77_CONST_CHAR_ARG2 (&uplo, 1),
532 F77_CONST_CHAR_ARG2 (&udiag, 1),
533 nr, tmp_data, nr, rcon,
534 work, iwork, dtrcon_info
535 F77_CHAR_ARG_LEN (1)
536 F77_CHAR_ARG_LEN (1)
537 F77_CHAR_ARG_LEN (1)));
538
539 if (dtrcon_info != 0)
540 info = -1;
541 }
542
543 if (info == -1 && ! force)
544 retval = *this; // Restore matrix contents.
545
546 return retval;
547 }
548
549 Matrix
finverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const550 Matrix::finverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
551 bool force, bool calc_cond) const
552 {
553 Matrix retval;
554
555 F77_INT nr = octave::to_f77_int (rows ());
556 F77_INT nc = octave::to_f77_int (cols ());
557
558 if (nr != nc || nr == 0 || nc == 0)
559 (*current_liboctave_error_handler) ("inverse requires square matrix");
560
561 Array<F77_INT> ipvt (dim_vector (nr, 1));
562 F77_INT *pipvt = ipvt.fortran_vec ();
563
564 retval = *this;
565 double *tmp_data = retval.fortran_vec ();
566
567 Array<double> z (dim_vector (1, 1));
568 F77_INT lwork = -1;
569
570 F77_INT tmp_info = 0;
571
572 // Query the optimum work array size.
573 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
574 z.fortran_vec (), lwork, tmp_info));
575
576 lwork = static_cast<F77_INT> (z(0));
577 lwork = (lwork < 4 * nc ? 4 * nc : lwork);
578 z.resize (dim_vector (lwork, 1));
579 double *pz = z.fortran_vec ();
580
581 info = 0;
582 tmp_info = 0;
583
584 // Calculate the norm of the matrix for later use when determining rcon.
585 double anorm;
586 if (calc_cond)
587 anorm = norm1 (retval);
588
589 F77_XFCN (dgetrf, DGETRF, (nc, nc, tmp_data, nr, pipvt, tmp_info));
590
591 info = tmp_info;
592
593 // Throw away extra info LAPACK gives so as to not change output.
594 rcon = 0.0;
595 if (info != 0)
596 info = -1;
597 else if (calc_cond)
598 {
599 F77_INT dgecon_info = 0;
600
601 // Now calculate the condition number for non-singular matrix.
602 char job = '1';
603 Array<F77_INT> iz (dim_vector (nc, 1));
604 F77_INT *piz = iz.fortran_vec ();
605 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
606 nc, tmp_data, nr, anorm,
607 rcon, pz, piz, dgecon_info
608 F77_CHAR_ARG_LEN (1)));
609
610 if (dgecon_info != 0)
611 info = -1;
612 }
613
614 if (info == -1 && ! force)
615 retval = *this; // Restore matrix contents.
616 else
617 {
618 F77_INT dgetri_info = 0;
619
620 F77_XFCN (dgetri, DGETRI, (nc, tmp_data, nr, pipvt,
621 pz, lwork, dgetri_info));
622
623 if (dgetri_info != 0)
624 info = -1;
625 }
626
627 if (info != 0)
628 mattype.mark_as_rectangular ();
629
630 return retval;
631 }
632
633 Matrix
inverse(MatrixType & mattype,octave_idx_type & info,double & rcon,bool force,bool calc_cond) const634 Matrix::inverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
635 bool force, bool calc_cond) const
636 {
637 int typ = mattype.type (false);
638 Matrix ret;
639
640 if (typ == MatrixType::Unknown)
641 typ = mattype.type (*this);
642
643 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
644 ret = tinverse (mattype, info, rcon, force, calc_cond);
645 else
646 {
647 if (mattype.ishermitian ())
648 {
649 octave::math::chol<Matrix> chol (*this, info, true, calc_cond);
650 if (info == 0)
651 {
652 if (calc_cond)
653 rcon = chol.rcond ();
654 else
655 rcon = 1.0;
656 ret = chol.inverse ();
657 }
658 else
659 mattype.mark_as_unsymmetric ();
660 }
661
662 if (! mattype.ishermitian ())
663 ret = finverse (mattype, info, rcon, force, calc_cond);
664
665 if ((calc_cond || mattype.ishermitian ()) && rcon == 0.0
666 && (numel () != 1))
667 ret = Matrix (rows (), columns (),
668 octave::numeric_limits<double>::Inf ());
669 }
670
671 return ret;
672 }
673
674 Matrix
pseudo_inverse(double tol) const675 Matrix::pseudo_inverse (double tol) const
676 {
677 octave::math::svd<Matrix> result (*this,
678 octave::math::svd<Matrix>::Type::economy);
679
680 DiagMatrix S = result.singular_values ();
681 Matrix U = result.left_singular_matrix ();
682 Matrix V = result.right_singular_matrix ();
683
684 ColumnVector sigma = S.extract_diag ();
685
686 octave_idx_type r = sigma.numel () - 1;
687 octave_idx_type nr = rows ();
688 octave_idx_type nc = cols ();
689
690 if (tol <= 0.0)
691 {
692 tol = std::max (nr, nc) * sigma.elem (0)
693 * std::numeric_limits<double>::epsilon ();
694
695 if (tol == 0)
696 tol = std::numeric_limits<double>::min ();
697 }
698
699 while (r >= 0 && sigma.elem (r) < tol)
700 r--;
701
702 if (r < 0)
703 return Matrix (nc, nr, 0.0);
704 else
705 {
706 Matrix Ur = U.extract (0, 0, nr-1, r);
707 DiagMatrix D = DiagMatrix (sigma.extract (0, r)).inverse ();
708 Matrix Vr = V.extract (0, 0, nc-1, r);
709 return Vr * D * Ur.transpose ();
710 }
711 }
712
713 #if defined (HAVE_FFTW)
714
715 ComplexMatrix
fourier(void) const716 Matrix::fourier (void) const
717 {
718 std::size_t nr = rows ();
719 std::size_t nc = cols ();
720
721 ComplexMatrix retval (nr, nc);
722
723 std::size_t npts, nsamples;
724
725 if (nr == 1 || nc == 1)
726 {
727 npts = (nr > nc ? nr : nc);
728 nsamples = 1;
729 }
730 else
731 {
732 npts = nr;
733 nsamples = nc;
734 }
735
736 const double *in (fortran_vec ());
737 Complex *out (retval.fortran_vec ());
738
739 octave::fftw::fft (in, out, npts, nsamples);
740
741 return retval;
742 }
743
744 ComplexMatrix
ifourier(void) const745 Matrix::ifourier (void) const
746 {
747 std::size_t nr = rows ();
748 std::size_t nc = cols ();
749
750 ComplexMatrix retval (nr, nc);
751
752 std::size_t npts, nsamples;
753
754 if (nr == 1 || nc == 1)
755 {
756 npts = (nr > nc ? nr : nc);
757 nsamples = 1;
758 }
759 else
760 {
761 npts = nr;
762 nsamples = nc;
763 }
764
765 ComplexMatrix tmp (*this);
766 Complex *in (tmp.fortran_vec ());
767 Complex *out (retval.fortran_vec ());
768
769 octave::fftw::ifft (in, out, npts, nsamples);
770
771 return retval;
772 }
773
774 ComplexMatrix
fourier2d(void) const775 Matrix::fourier2d (void) const
776 {
777 dim_vector dv (rows (), cols ());
778
779 const double *in = fortran_vec ();
780 ComplexMatrix retval (rows (), cols ());
781 octave::fftw::fftNd (in, retval.fortran_vec (), 2, dv);
782
783 return retval;
784 }
785
786 ComplexMatrix
ifourier2d(void) const787 Matrix::ifourier2d (void) const
788 {
789 dim_vector dv (rows (), cols ());
790
791 ComplexMatrix retval (*this);
792 Complex *out (retval.fortran_vec ());
793
794 octave::fftw::ifftNd (out, out, 2, dv);
795
796 return retval;
797 }
798
799 #else
800
801 ComplexMatrix
fourier(void) const802 Matrix::fourier (void) const
803 {
804 (*current_liboctave_error_handler)
805 ("support for FFTW was unavailable or disabled when liboctave was built");
806
807 return ComplexMatrix ();
808 }
809
810 ComplexMatrix
ifourier(void) const811 Matrix::ifourier (void) const
812 {
813 (*current_liboctave_error_handler)
814 ("support for FFTW was unavailable or disabled when liboctave was built");
815
816 return ComplexMatrix ();
817 }
818
819 ComplexMatrix
fourier2d(void) const820 Matrix::fourier2d (void) const
821 {
822 (*current_liboctave_error_handler)
823 ("support for FFTW was unavailable or disabled when liboctave was built");
824
825 return ComplexMatrix ();
826 }
827
828 ComplexMatrix
ifourier2d(void) const829 Matrix::ifourier2d (void) const
830 {
831 (*current_liboctave_error_handler)
832 ("support for FFTW was unavailable or disabled when liboctave was built");
833
834 return ComplexMatrix ();
835 }
836
837 #endif
838
839 DET
determinant(void) const840 Matrix::determinant (void) const
841 {
842 octave_idx_type info;
843 double rcon;
844 return determinant (info, rcon, 0);
845 }
846
847 DET
determinant(octave_idx_type & info) const848 Matrix::determinant (octave_idx_type& info) const
849 {
850 double rcon;
851 return determinant (info, rcon, 0);
852 }
853
854 DET
determinant(octave_idx_type & info,double & rcon,bool calc_cond) const855 Matrix::determinant (octave_idx_type& info, double& rcon, bool calc_cond) const
856 {
857 MatrixType mattype (*this);
858 return determinant (mattype, info, rcon, calc_cond);
859 }
860
861 DET
determinant(MatrixType & mattype,octave_idx_type & info,double & rcon,bool calc_cond) const862 Matrix::determinant (MatrixType& mattype,
863 octave_idx_type& info, double& rcon, bool calc_cond) const
864 {
865 DET retval (1.0);
866
867 info = 0;
868 rcon = 0.0;
869
870 F77_INT nr = octave::to_f77_int (rows ());
871 F77_INT nc = octave::to_f77_int (cols ());
872
873 if (nr != nc)
874 (*current_liboctave_error_handler) ("matrix must be square");
875
876 volatile int typ = mattype.type ();
877
878 // Even though the matrix is marked as singular (Rectangular), we may still
879 // get a useful number from the LU factorization, because it always completes.
880
881 if (typ == MatrixType::Unknown)
882 typ = mattype.type (*this);
883 else if (typ == MatrixType::Rectangular)
884 typ = MatrixType::Full;
885
886 if (typ == MatrixType::Lower || typ == MatrixType::Upper)
887 {
888 for (F77_INT i = 0; i < nc; i++)
889 retval *= elem (i,i);
890 }
891 else if (typ == MatrixType::Hermitian)
892 {
893 Matrix atmp = *this;
894 double *tmp_data = atmp.fortran_vec ();
895
896 // Calculate the norm of the matrix for later use when determining rcon.
897 double anorm;
898 if (calc_cond)
899 anorm = norm1 (*this);
900
901 F77_INT tmp_info = 0;
902
903 char job = 'L';
904 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
905 tmp_data, nr, tmp_info
906 F77_CHAR_ARG_LEN (1)));
907
908 info = tmp_info;
909
910 if (info != 0)
911 {
912 rcon = 0.0;
913 mattype.mark_as_unsymmetric ();
914 typ = MatrixType::Full;
915 }
916 else
917 {
918 if (calc_cond)
919 {
920 Array<double> z (dim_vector (3 * nc, 1));
921 double *pz = z.fortran_vec ();
922 Array<F77_INT> iz (dim_vector (nc, 1));
923 F77_INT *piz = iz.fortran_vec ();
924
925 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
926 nr, tmp_data, nr, anorm,
927 rcon, pz, piz, tmp_info
928 F77_CHAR_ARG_LEN (1)));
929
930 info = tmp_info;
931
932 if (info != 0)
933 rcon = 0.0;
934 }
935
936 for (F77_INT i = 0; i < nc; i++)
937 retval *= atmp(i,i);
938
939 retval = retval.square ();
940 }
941 }
942 else if (typ != MatrixType::Full)
943 (*current_liboctave_error_handler) ("det: invalid dense matrix type");
944
945 if (typ == MatrixType::Full)
946 {
947 Array<F77_INT> ipvt (dim_vector (nr, 1));
948 F77_INT *pipvt = ipvt.fortran_vec ();
949
950 Matrix atmp = *this;
951 double *tmp_data = atmp.fortran_vec ();
952
953 info = 0;
954 F77_INT tmp_info = 0;
955
956 // Calculate the norm of the matrix for later use when determining rcon.
957 double anorm;
958 if (calc_cond)
959 anorm = norm1 (*this);
960
961 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, tmp_info));
962
963 info = tmp_info;
964
965 // Throw away extra info LAPACK gives so as to not change output.
966 rcon = 0.0;
967 if (info != 0)
968 {
969 info = -1;
970 retval = DET ();
971 }
972 else
973 {
974 if (calc_cond)
975 {
976 // Now calc the condition number for non-singular matrix.
977 char job = '1';
978 Array<double> z (dim_vector (4 * nc, 1));
979 double *pz = z.fortran_vec ();
980 Array<F77_INT> iz (dim_vector (nc, 1));
981 F77_INT *piz = iz.fortran_vec ();
982
983 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
984 nc, tmp_data, nr, anorm,
985 rcon, pz, piz, tmp_info
986 F77_CHAR_ARG_LEN (1)));
987
988 info = tmp_info;
989 }
990
991 if (info != 0)
992 {
993 info = -1;
994 retval = DET ();
995 }
996 else
997 {
998 for (F77_INT i = 0; i < nc; i++)
999 {
1000 double c = atmp(i,i);
1001 retval *= (ipvt(i) != (i+1)) ? -c : c;
1002 }
1003 }
1004 }
1005 }
1006
1007 return retval;
1008 }
1009
1010 double
rcond(void) const1011 Matrix::rcond (void) const
1012 {
1013 MatrixType mattype (*this);
1014 return rcond (mattype);
1015 }
1016
1017 double
rcond(MatrixType & mattype) const1018 Matrix::rcond (MatrixType& mattype) const
1019 {
1020 double rcon = octave::numeric_limits<double>::NaN ();
1021 F77_INT nr = octave::to_f77_int (rows ());
1022 F77_INT nc = octave::to_f77_int (cols ());
1023
1024 if (nr != nc)
1025 (*current_liboctave_error_handler) ("matrix must be square");
1026
1027 if (nr == 0 || nc == 0)
1028 rcon = octave::numeric_limits<double>::Inf ();
1029 else
1030 {
1031 volatile int typ = mattype.type ();
1032
1033 if (typ == MatrixType::Unknown)
1034 typ = mattype.type (*this);
1035
1036 // Only calculate the condition number for LU/Cholesky
1037 if (typ == MatrixType::Upper)
1038 {
1039 const double *tmp_data = fortran_vec ();
1040 F77_INT info = 0;
1041 char norm = '1';
1042 char uplo = 'U';
1043 char dia = 'N';
1044
1045 Array<double> z (dim_vector (3 * nc, 1));
1046 double *pz = z.fortran_vec ();
1047 Array<F77_INT> iz (dim_vector (nc, 1));
1048 F77_INT *piz = iz.fortran_vec ();
1049
1050 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1051 F77_CONST_CHAR_ARG2 (&uplo, 1),
1052 F77_CONST_CHAR_ARG2 (&dia, 1),
1053 nr, tmp_data, nr, rcon,
1054 pz, piz, info
1055 F77_CHAR_ARG_LEN (1)
1056 F77_CHAR_ARG_LEN (1)
1057 F77_CHAR_ARG_LEN (1)));
1058
1059 if (info != 0)
1060 rcon = 0.0;
1061 }
1062 else if (typ == MatrixType::Permuted_Upper)
1063 (*current_liboctave_error_handler)
1064 ("permuted triangular matrix not implemented");
1065 else if (typ == MatrixType::Lower)
1066 {
1067 const double *tmp_data = fortran_vec ();
1068 F77_INT info = 0;
1069 char norm = '1';
1070 char uplo = 'L';
1071 char dia = 'N';
1072
1073 Array<double> z (dim_vector (3 * nc, 1));
1074 double *pz = z.fortran_vec ();
1075 Array<F77_INT> iz (dim_vector (nc, 1));
1076 F77_INT *piz = iz.fortran_vec ();
1077
1078 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1079 F77_CONST_CHAR_ARG2 (&uplo, 1),
1080 F77_CONST_CHAR_ARG2 (&dia, 1),
1081 nr, tmp_data, nr, rcon,
1082 pz, piz, info
1083 F77_CHAR_ARG_LEN (1)
1084 F77_CHAR_ARG_LEN (1)
1085 F77_CHAR_ARG_LEN (1)));
1086
1087 if (info != 0)
1088 rcon = 0.0;
1089 }
1090 else if (typ == MatrixType::Permuted_Lower)
1091 (*current_liboctave_error_handler)
1092 ("permuted triangular matrix not implemented");
1093 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1094 {
1095 double anorm = -1.0;
1096
1097 if (typ == MatrixType::Hermitian)
1098 {
1099 F77_INT info = 0;
1100 char job = 'L';
1101
1102 Matrix atmp = *this;
1103 double *tmp_data = atmp.fortran_vec ();
1104
1105 anorm = norm1 (atmp);
1106
1107 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1108 tmp_data, nr, info
1109 F77_CHAR_ARG_LEN (1)));
1110
1111 if (info != 0)
1112 {
1113 rcon = 0.0;
1114 mattype.mark_as_unsymmetric ();
1115 typ = MatrixType::Full;
1116 }
1117 else
1118 {
1119 Array<double> z (dim_vector (3 * nc, 1));
1120 double *pz = z.fortran_vec ();
1121 Array<F77_INT> iz (dim_vector (nc, 1));
1122 F77_INT *piz = iz.fortran_vec ();
1123
1124 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1125 nr, tmp_data, nr, anorm,
1126 rcon, pz, piz, info
1127 F77_CHAR_ARG_LEN (1)));
1128
1129 if (info != 0)
1130 rcon = 0.0;
1131 }
1132 }
1133
1134 if (typ == MatrixType::Full)
1135 {
1136 F77_INT info = 0;
1137
1138 Matrix atmp = *this;
1139 double *tmp_data = atmp.fortran_vec ();
1140
1141 Array<F77_INT> ipvt (dim_vector (nr, 1));
1142 F77_INT *pipvt = ipvt.fortran_vec ();
1143
1144 if (anorm < 0.0)
1145 anorm = norm1 (atmp);
1146
1147 Array<double> z (dim_vector (4 * nc, 1));
1148 double *pz = z.fortran_vec ();
1149 Array<F77_INT> iz (dim_vector (nc, 1));
1150 F77_INT *piz = iz.fortran_vec ();
1151
1152 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
1153
1154 if (info != 0)
1155 {
1156 rcon = 0.0;
1157 mattype.mark_as_rectangular ();
1158 }
1159 else
1160 {
1161 char job = '1';
1162 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1163 nc, tmp_data, nr, anorm,
1164 rcon, pz, piz, info
1165 F77_CHAR_ARG_LEN (1)));
1166
1167 if (info != 0)
1168 rcon = 0.0;
1169 }
1170 }
1171 }
1172 else
1173 rcon = 0.0;
1174 }
1175
1176 return rcon;
1177 }
1178
1179 Matrix
utsolve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond,blas_trans_type transt) const1180 Matrix::utsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
1181 double& rcon, solve_singularity_handler sing_handler,
1182 bool calc_cond, blas_trans_type transt) const
1183 {
1184 Matrix retval;
1185
1186 F77_INT nr = octave::to_f77_int (rows ());
1187 F77_INT nc = octave::to_f77_int (cols ());
1188
1189 F77_INT b_nr = octave::to_f77_int (b.rows ());
1190 F77_INT b_nc = octave::to_f77_int (b.cols ());
1191
1192 if (nr != b_nr)
1193 (*current_liboctave_error_handler)
1194 ("matrix dimension mismatch solution of linear equations");
1195
1196 if (nr == 0 || nc == 0 || b_nc == 0)
1197 retval = Matrix (nc, b_nc, 0.0);
1198 else
1199 {
1200 volatile int typ = mattype.type ();
1201
1202 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1203 (*current_liboctave_error_handler) ("incorrect matrix type");
1204
1205 rcon = 1.0;
1206 info = 0;
1207
1208 if (typ == MatrixType::Permuted_Upper)
1209 (*current_liboctave_error_handler)
1210 ("permuted triangular matrix not implemented");
1211
1212 const double *tmp_data = fortran_vec ();
1213
1214 retval = b;
1215 double *result = retval.fortran_vec ();
1216
1217 char uplo = 'U';
1218 char trans = get_blas_char (transt);
1219 char dia = 'N';
1220
1221 F77_INT tmp_info = 0;
1222
1223 F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1224 F77_CONST_CHAR_ARG2 (&trans, 1),
1225 F77_CONST_CHAR_ARG2 (&dia, 1),
1226 nr, b_nc, tmp_data, nr,
1227 result, nr, tmp_info
1228 F77_CHAR_ARG_LEN (1)
1229 F77_CHAR_ARG_LEN (1)
1230 F77_CHAR_ARG_LEN (1)));
1231
1232 info = tmp_info;
1233
1234 if (calc_cond)
1235 {
1236 char norm = '1';
1237 uplo = 'U';
1238 dia = 'N';
1239
1240 Array<double> z (dim_vector (3 * nc, 1));
1241 double *pz = z.fortran_vec ();
1242 Array<F77_INT> iz (dim_vector (nc, 1));
1243 F77_INT *piz = iz.fortran_vec ();
1244
1245 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1246 F77_CONST_CHAR_ARG2 (&uplo, 1),
1247 F77_CONST_CHAR_ARG2 (&dia, 1),
1248 nr, tmp_data, nr, rcon,
1249 pz, piz, tmp_info
1250 F77_CHAR_ARG_LEN (1)
1251 F77_CHAR_ARG_LEN (1)
1252 F77_CHAR_ARG_LEN (1)));
1253
1254 info = tmp_info;
1255
1256 if (info != 0)
1257 info = -2;
1258
1259 // FIXME: Why calculate this, rather than just compare to 0.0?
1260 volatile double rcond_plus_one = rcon + 1.0;
1261
1262 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1263 {
1264 info = -2;
1265
1266 if (sing_handler)
1267 sing_handler (rcon);
1268 else
1269 octave::warn_singular_matrix (rcon);
1270 }
1271 }
1272 }
1273
1274 return retval;
1275 }
1276
1277 Matrix
ltsolve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond,blas_trans_type transt) const1278 Matrix::ltsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
1279 double& rcon, solve_singularity_handler sing_handler,
1280 bool calc_cond, blas_trans_type transt) const
1281 {
1282 Matrix retval;
1283
1284 F77_INT nr = octave::to_f77_int (rows ());
1285 F77_INT nc = octave::to_f77_int (cols ());
1286
1287 F77_INT b_nr = octave::to_f77_int (b.rows ());
1288 F77_INT b_nc = octave::to_f77_int (b.cols ());
1289
1290 if (nr != b_nr)
1291 (*current_liboctave_error_handler)
1292 ("matrix dimension mismatch solution of linear equations");
1293
1294 if (nr == 0 || nc == 0 || b_nc == 0)
1295 retval = Matrix (nc, b_nc, 0.0);
1296 else
1297 {
1298 volatile int typ = mattype.type ();
1299
1300 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
1301 (*current_liboctave_error_handler) ("incorrect matrix type");
1302
1303 rcon = 1.0;
1304 info = 0;
1305
1306 if (typ == MatrixType::Permuted_Lower)
1307 (*current_liboctave_error_handler)
1308 ("permuted triangular matrix not implemented");
1309
1310 const double *tmp_data = fortran_vec ();
1311
1312 retval = b;
1313 double *result = retval.fortran_vec ();
1314
1315 char uplo = 'L';
1316 char trans = get_blas_char (transt);
1317 char dia = 'N';
1318
1319 F77_INT tmp_info = 0;
1320
1321 F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
1322 F77_CONST_CHAR_ARG2 (&trans, 1),
1323 F77_CONST_CHAR_ARG2 (&dia, 1),
1324 nr, b_nc, tmp_data, nr,
1325 result, nr, tmp_info
1326 F77_CHAR_ARG_LEN (1)
1327 F77_CHAR_ARG_LEN (1)
1328 F77_CHAR_ARG_LEN (1)));
1329
1330 info = tmp_info;
1331
1332 if (calc_cond)
1333 {
1334 char norm = '1';
1335 uplo = 'L';
1336 dia = 'N';
1337
1338 Array<double> z (dim_vector (3 * nc, 1));
1339 double *pz = z.fortran_vec ();
1340 Array<F77_INT> iz (dim_vector (nc, 1));
1341 F77_INT *piz = iz.fortran_vec ();
1342
1343 F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
1344 F77_CONST_CHAR_ARG2 (&uplo, 1),
1345 F77_CONST_CHAR_ARG2 (&dia, 1),
1346 nr, tmp_data, nr, rcon,
1347 pz, piz, tmp_info
1348 F77_CHAR_ARG_LEN (1)
1349 F77_CHAR_ARG_LEN (1)
1350 F77_CHAR_ARG_LEN (1)));
1351
1352 info = tmp_info;
1353
1354 if (info != 0)
1355 info = -2;
1356
1357 volatile double rcond_plus_one = rcon + 1.0;
1358
1359 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1360 {
1361 info = -2;
1362
1363 if (sing_handler)
1364 sing_handler (rcon);
1365 else
1366 octave::warn_singular_matrix (rcon);
1367 }
1368 }
1369 }
1370
1371 return retval;
1372 }
1373
1374 Matrix
fsolve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool calc_cond) const1375 Matrix::fsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
1376 double& rcon, solve_singularity_handler sing_handler,
1377 bool calc_cond) const
1378 {
1379 Matrix retval;
1380
1381 F77_INT nr = octave::to_f77_int (rows ());
1382 F77_INT nc = octave::to_f77_int (cols ());
1383
1384 if (nr != nc || nr != b.rows ())
1385 (*current_liboctave_error_handler)
1386 ("matrix dimension mismatch solution of linear equations");
1387
1388 if (nr == 0 || b.cols () == 0)
1389 retval = Matrix (nc, b.cols (), 0.0);
1390 else
1391 {
1392 volatile int typ = mattype.type ();
1393
1394 // Calculate the norm of the matrix for later use when determining rcon.
1395 double anorm = -1.0;
1396
1397 if (typ == MatrixType::Hermitian)
1398 {
1399 info = 0;
1400 char job = 'L';
1401
1402 Matrix atmp = *this;
1403 double *tmp_data = atmp.fortran_vec ();
1404
1405 // The norm of the matrix for later use when determining rcon.
1406 if (calc_cond)
1407 anorm = norm1 (atmp);
1408
1409 F77_INT tmp_info = 0;
1410
1411 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
1412 tmp_data, nr, tmp_info
1413 F77_CHAR_ARG_LEN (1)));
1414
1415 info = tmp_info;
1416
1417 // Throw away extra info LAPACK gives so as to not change output.
1418 rcon = 0.0;
1419 if (info != 0)
1420 {
1421 info = -2;
1422
1423 mattype.mark_as_unsymmetric ();
1424 typ = MatrixType::Full;
1425 }
1426 else
1427 {
1428 if (calc_cond)
1429 {
1430 Array<double> z (dim_vector (3 * nc, 1));
1431 double *pz = z.fortran_vec ();
1432 Array<F77_INT> iz (dim_vector (nc, 1));
1433 F77_INT *piz = iz.fortran_vec ();
1434
1435 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
1436 nr, tmp_data, nr, anorm,
1437 rcon, pz, piz, tmp_info
1438 F77_CHAR_ARG_LEN (1)));
1439
1440 info = tmp_info;
1441
1442 if (info != 0)
1443 info = -2;
1444
1445 volatile double rcond_plus_one = rcon + 1.0;
1446
1447 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1448 {
1449 info = -2;
1450
1451 if (sing_handler)
1452 sing_handler (rcon);
1453 else
1454 octave::warn_singular_matrix (rcon);
1455 }
1456 }
1457
1458 if (info == 0)
1459 {
1460 retval = b;
1461 double *result = retval.fortran_vec ();
1462
1463 F77_INT b_nr = octave::to_f77_int (b.rows ());
1464 F77_INT b_nc = octave::to_f77_int (b.cols ());
1465
1466 F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1467 nr, b_nc, tmp_data, nr,
1468 result, b_nr, tmp_info
1469 F77_CHAR_ARG_LEN (1)));
1470
1471 info = tmp_info;
1472 }
1473 else
1474 {
1475 mattype.mark_as_unsymmetric ();
1476 typ = MatrixType::Full;
1477 }
1478 }
1479 }
1480
1481 if (typ == MatrixType::Full)
1482 {
1483 info = 0;
1484
1485 Array<F77_INT> ipvt (dim_vector (nr, 1));
1486 F77_INT *pipvt = ipvt.fortran_vec ();
1487
1488 Matrix atmp = *this;
1489 double *tmp_data = atmp.fortran_vec ();
1490
1491 if (calc_cond && anorm < 0.0)
1492 anorm = norm1 (atmp);
1493
1494 Array<double> z (dim_vector (4 * nc, 1));
1495 double *pz = z.fortran_vec ();
1496 Array<F77_INT> iz (dim_vector (nc, 1));
1497 F77_INT *piz = iz.fortran_vec ();
1498
1499 F77_INT tmp_info = 0;
1500
1501 F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, tmp_info));
1502
1503 info = tmp_info;
1504
1505 // Throw away extra info LAPACK gives so as to not change output.
1506 rcon = 0.0;
1507 if (info != 0)
1508 {
1509 info = -2;
1510
1511 if (sing_handler)
1512 sing_handler (rcon);
1513 else
1514 octave::warn_singular_matrix ();
1515
1516 mattype.mark_as_rectangular ();
1517 }
1518 else
1519 {
1520 if (calc_cond)
1521 {
1522 // Calculate the condition number for non-singular matrix.
1523 char job = '1';
1524 F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
1525 nc, tmp_data, nr, anorm,
1526 rcon, pz, piz, tmp_info
1527 F77_CHAR_ARG_LEN (1)));
1528
1529 info = tmp_info;
1530
1531 if (info != 0)
1532 info = -2;
1533
1534 volatile double rcond_plus_one = rcon + 1.0;
1535
1536 if (rcond_plus_one == 1.0 || octave::math::isnan (rcon))
1537 {
1538 info = -2;
1539
1540 if (sing_handler)
1541 sing_handler (rcon);
1542 else
1543 octave::warn_singular_matrix (rcon);
1544 }
1545 }
1546
1547 if (info == 0)
1548 {
1549 retval = b;
1550 double *result = retval.fortran_vec ();
1551
1552 F77_INT b_nr = octave::to_f77_int (b.rows ());
1553 F77_INT b_nc = octave::to_f77_int (b.cols ());
1554
1555 char job = 'N';
1556 F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
1557 nr, b_nc, tmp_data, nr,
1558 pipvt, result, b_nr, tmp_info
1559 F77_CHAR_ARG_LEN (1)));
1560
1561 info = tmp_info;
1562 }
1563 else
1564 mattype.mark_as_rectangular ();
1565 }
1566 }
1567 else if (typ != MatrixType::Hermitian)
1568 (*current_liboctave_error_handler) ("incorrect matrix type");
1569 }
1570
1571 return retval;
1572 }
1573
1574 Matrix
solve(MatrixType & mattype,const Matrix & b) const1575 Matrix::solve (MatrixType& mattype, const Matrix& b) const
1576 {
1577 octave_idx_type info;
1578 double rcon;
1579 return solve (mattype, b, info, rcon, nullptr);
1580 }
1581
1582 Matrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info) const1583 Matrix::solve (MatrixType& mattype, const Matrix& b,
1584 octave_idx_type& info) const
1585 {
1586 double rcon;
1587 return solve (mattype, b, info, rcon, nullptr);
1588 }
1589
1590 Matrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon) const1591 Matrix::solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
1592 double& rcon) const
1593 {
1594 return solve (mattype, b, info, rcon, nullptr);
1595 }
1596
1597 Matrix
solve(MatrixType & mattype,const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool singular_fallback,blas_trans_type transt) const1598 Matrix::solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
1599 double& rcon, solve_singularity_handler sing_handler,
1600 bool singular_fallback, blas_trans_type transt) const
1601 {
1602 Matrix retval;
1603 int typ = mattype.type ();
1604
1605 if (typ == MatrixType::Unknown)
1606 typ = mattype.type (*this);
1607
1608 // Only calculate the condition number for LU/Cholesky
1609 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
1610 retval = utsolve (mattype, b, info, rcon, sing_handler, true, transt);
1611 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
1612 retval = ltsolve (mattype, b, info, rcon, sing_handler, true, transt);
1613 else if (transt == blas_trans || transt == blas_conj_trans)
1614 return transpose ().solve (mattype, b, info, rcon, sing_handler,
1615 singular_fallback);
1616 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
1617 retval = fsolve (mattype, b, info, rcon, sing_handler, true);
1618 else if (typ != MatrixType::Rectangular)
1619 (*current_liboctave_error_handler) ("unknown matrix type");
1620
1621 // Rectangular or one of the above solvers flags a singular matrix
1622 if (singular_fallback && mattype.type () == MatrixType::Rectangular)
1623 {
1624 octave_idx_type rank;
1625 retval = lssolve (b, info, rank, rcon);
1626 }
1627
1628 return retval;
1629 }
1630
1631 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b) const1632 Matrix::solve (MatrixType& mattype, const ComplexMatrix& b) const
1633 {
1634 octave_idx_type info;
1635 double rcon;
1636 return solve (mattype, b, info, rcon, nullptr);
1637 }
1638
1639 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info) const1640 Matrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1641 octave_idx_type& info) const
1642 {
1643 double rcon;
1644 return solve (mattype, b, info, rcon, nullptr);
1645 }
1646
1647 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon) const1648 Matrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1649 octave_idx_type& info, double& rcon) const
1650 {
1651 return solve (mattype, b, info, rcon, nullptr);
1652 }
1653
1654 static Matrix
stack_complex_matrix(const ComplexMatrix & cm)1655 stack_complex_matrix (const ComplexMatrix& cm)
1656 {
1657 octave_idx_type m = cm.rows ();
1658 octave_idx_type n = cm.cols ();
1659 octave_idx_type nel = m*n;
1660 Matrix retval (m, 2*n);
1661 const Complex *cmd = cm.data ();
1662 double *rd = retval.fortran_vec ();
1663 for (octave_idx_type i = 0; i < nel; i++)
1664 {
1665 rd[i] = std::real (cmd[i]);
1666 rd[nel+i] = std::imag (cmd[i]);
1667 }
1668 return retval;
1669 }
1670
1671 static ComplexMatrix
unstack_complex_matrix(const Matrix & sm)1672 unstack_complex_matrix (const Matrix& sm)
1673 {
1674 octave_idx_type m = sm.rows ();
1675 octave_idx_type n = sm.cols () / 2;
1676 octave_idx_type nel = m*n;
1677 ComplexMatrix retval (m, n);
1678 const double *smd = sm.data ();
1679 Complex *rd = retval.fortran_vec ();
1680 for (octave_idx_type i = 0; i < nel; i++)
1681 rd[i] = Complex (smd[i], smd[nel+i]);
1682 return retval;
1683 }
1684
1685 ComplexMatrix
solve(MatrixType & mattype,const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,bool singular_fallback,blas_trans_type transt) const1686 Matrix::solve (MatrixType& mattype, const ComplexMatrix& b,
1687 octave_idx_type& info, double& rcon,
1688 solve_singularity_handler sing_handler,
1689 bool singular_fallback, blas_trans_type transt) const
1690 {
1691 Matrix tmp = stack_complex_matrix (b);
1692 tmp = solve (mattype, tmp, info, rcon, sing_handler, singular_fallback,
1693 transt);
1694 return unstack_complex_matrix (tmp);
1695 }
1696
1697 ColumnVector
solve(MatrixType & mattype,const ColumnVector & b) const1698 Matrix::solve (MatrixType& mattype, const ColumnVector& b) const
1699 {
1700 octave_idx_type info; double rcon;
1701 return solve (mattype, b, info, rcon);
1702 }
1703
1704 ColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info) const1705 Matrix::solve (MatrixType& mattype, const ColumnVector& b,
1706 octave_idx_type& info) const
1707 {
1708 double rcon;
1709 return solve (mattype, b, info, rcon);
1710 }
1711
1712 ColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info,double & rcon) const1713 Matrix::solve (MatrixType& mattype, const ColumnVector& b,
1714 octave_idx_type& info, double& rcon) const
1715 {
1716 return solve (mattype, b, info, rcon, nullptr);
1717 }
1718
1719 ColumnVector
solve(MatrixType & mattype,const ColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1720 Matrix::solve (MatrixType& mattype, const ColumnVector& b,
1721 octave_idx_type& info, double& rcon,
1722 solve_singularity_handler sing_handler,
1723 blas_trans_type transt) const
1724 {
1725 Matrix tmp (b);
1726 tmp = solve (mattype, tmp, info, rcon, sing_handler, true, transt);
1727 return tmp.column (static_cast<octave_idx_type> (0));
1728 }
1729
1730 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b) const1731 Matrix::solve (MatrixType& mattype, const ComplexColumnVector& b) const
1732 {
1733 ComplexMatrix tmp (*this);
1734 return tmp.solve (mattype, b);
1735 }
1736
1737 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info) const1738 Matrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
1739 octave_idx_type& info) const
1740 {
1741 ComplexMatrix tmp (*this);
1742 return tmp.solve (mattype, b, info);
1743 }
1744
1745 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info,double & rcon) const1746 Matrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
1747 octave_idx_type& info, double& rcon) const
1748 {
1749 ComplexMatrix tmp (*this);
1750 return tmp.solve (mattype, b, info, rcon);
1751 }
1752
1753 ComplexColumnVector
solve(MatrixType & mattype,const ComplexColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1754 Matrix::solve (MatrixType& mattype, const ComplexColumnVector& b,
1755 octave_idx_type& info, double& rcon,
1756 solve_singularity_handler sing_handler,
1757 blas_trans_type transt) const
1758 {
1759 ComplexMatrix tmp (*this);
1760 return tmp.solve (mattype, b, info, rcon, sing_handler, transt);
1761 }
1762
1763 Matrix
solve(const Matrix & b) const1764 Matrix::solve (const Matrix& b) const
1765 {
1766 octave_idx_type info;
1767 double rcon;
1768 return solve (b, info, rcon, nullptr);
1769 }
1770
1771 Matrix
solve(const Matrix & b,octave_idx_type & info) const1772 Matrix::solve (const Matrix& b, octave_idx_type& info) const
1773 {
1774 double rcon;
1775 return solve (b, info, rcon, nullptr);
1776 }
1777
1778 Matrix
solve(const Matrix & b,octave_idx_type & info,double & rcon) const1779 Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcon) const
1780 {
1781 return solve (b, info, rcon, nullptr);
1782 }
1783
1784 Matrix
solve(const Matrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1785 Matrix::solve (const Matrix& b, octave_idx_type& info,
1786 double& rcon, solve_singularity_handler sing_handler,
1787 blas_trans_type transt) const
1788 {
1789 MatrixType mattype (*this);
1790 return solve (mattype, b, info, rcon, sing_handler, true, transt);
1791 }
1792
1793 ComplexMatrix
solve(const ComplexMatrix & b) const1794 Matrix::solve (const ComplexMatrix& b) const
1795 {
1796 ComplexMatrix tmp (*this);
1797 return tmp.solve (b);
1798 }
1799
1800 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info) const1801 Matrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
1802 {
1803 ComplexMatrix tmp (*this);
1804 return tmp.solve (b, info);
1805 }
1806
1807 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info,double & rcon) const1808 Matrix::solve (const ComplexMatrix& b, octave_idx_type& info,
1809 double& rcon) const
1810 {
1811 ComplexMatrix tmp (*this);
1812 return tmp.solve (b, info, rcon);
1813 }
1814
1815 ComplexMatrix
solve(const ComplexMatrix & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1816 Matrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon,
1817 solve_singularity_handler sing_handler,
1818 blas_trans_type transt) const
1819 {
1820 ComplexMatrix tmp (*this);
1821 return tmp.solve (b, info, rcon, sing_handler, transt);
1822 }
1823
1824 ColumnVector
solve(const ColumnVector & b) const1825 Matrix::solve (const ColumnVector& b) const
1826 {
1827 octave_idx_type info; double rcon;
1828 return solve (b, info, rcon);
1829 }
1830
1831 ColumnVector
solve(const ColumnVector & b,octave_idx_type & info) const1832 Matrix::solve (const ColumnVector& b, octave_idx_type& info) const
1833 {
1834 double rcon;
1835 return solve (b, info, rcon);
1836 }
1837
1838 ColumnVector
solve(const ColumnVector & b,octave_idx_type & info,double & rcon) const1839 Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const
1840 {
1841 return solve (b, info, rcon, nullptr);
1842 }
1843
1844 ColumnVector
solve(const ColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1845 Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcon,
1846 solve_singularity_handler sing_handler,
1847 blas_trans_type transt) const
1848 {
1849 MatrixType mattype (*this);
1850 return solve (mattype, b, info, rcon, sing_handler, transt);
1851 }
1852
1853 ComplexColumnVector
solve(const ComplexColumnVector & b) const1854 Matrix::solve (const ComplexColumnVector& b) const
1855 {
1856 ComplexMatrix tmp (*this);
1857 return tmp.solve (b);
1858 }
1859
1860 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info) const1861 Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const
1862 {
1863 ComplexMatrix tmp (*this);
1864 return tmp.solve (b, info);
1865 }
1866
1867 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info,double & rcon) const1868 Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
1869 double& rcon) const
1870 {
1871 ComplexMatrix tmp (*this);
1872 return tmp.solve (b, info, rcon);
1873 }
1874
1875 ComplexColumnVector
solve(const ComplexColumnVector & b,octave_idx_type & info,double & rcon,solve_singularity_handler sing_handler,blas_trans_type transt) const1876 Matrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
1877 double& rcon,
1878 solve_singularity_handler sing_handler,
1879 blas_trans_type transt) const
1880 {
1881 ComplexMatrix tmp (*this);
1882 return tmp.solve (b, info, rcon, sing_handler, transt);
1883 }
1884
1885 Matrix
lssolve(const Matrix & b) const1886 Matrix::lssolve (const Matrix& b) const
1887 {
1888 octave_idx_type info;
1889 octave_idx_type rank;
1890 double rcon;
1891 return lssolve (b, info, rank, rcon);
1892 }
1893
1894 Matrix
lssolve(const Matrix & b,octave_idx_type & info) const1895 Matrix::lssolve (const Matrix& b, octave_idx_type& info) const
1896 {
1897 octave_idx_type rank;
1898 double rcon;
1899 return lssolve (b, info, rank, rcon);
1900 }
1901
1902 Matrix
lssolve(const Matrix & b,octave_idx_type & info,octave_idx_type & rank) const1903 Matrix::lssolve (const Matrix& b, octave_idx_type& info,
1904 octave_idx_type& rank) const
1905 {
1906 double rcon;
1907 return lssolve (b, info, rank, rcon);
1908 }
1909
1910 Matrix
lssolve(const Matrix & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const1911 Matrix::lssolve (const Matrix& b, octave_idx_type& info,
1912 octave_idx_type& rank, double& rcon) const
1913 {
1914 Matrix retval;
1915
1916 F77_INT m = octave::to_f77_int (rows ());
1917 F77_INT n = octave::to_f77_int (cols ());
1918
1919 F77_INT b_nr = octave::to_f77_int (b.rows ());
1920 F77_INT b_nc = octave::to_f77_int (b.cols ());
1921 F77_INT nrhs = b_nc; // alias for code readability
1922
1923 if (m != b_nr)
1924 (*current_liboctave_error_handler)
1925 ("matrix dimension mismatch solution of linear equations");
1926
1927 if (m == 0 || n == 0 || b_nc == 0)
1928 retval = Matrix (n, b_nc, 0.0);
1929 else
1930 {
1931 volatile F77_INT minmn = (m < n ? m : n);
1932 F77_INT maxmn = (m > n ? m : n);
1933 rcon = -1.0;
1934 if (m != n)
1935 {
1936 retval = Matrix (maxmn, nrhs, 0.0);
1937
1938 for (F77_INT j = 0; j < nrhs; j++)
1939 for (F77_INT i = 0; i < m; i++)
1940 retval.elem (i, j) = b.elem (i, j);
1941 }
1942 else
1943 retval = b;
1944
1945 Matrix atmp = *this;
1946 double *tmp_data = atmp.fortran_vec ();
1947
1948 double *pretval = retval.fortran_vec ();
1949 Array<double> s (dim_vector (minmn, 1));
1950 double *ps = s.fortran_vec ();
1951
1952 // Ask DGELSD what the dimension of WORK should be.
1953 F77_INT lwork = -1;
1954
1955 Array<double> work (dim_vector (1, 1));
1956
1957 F77_INT smlsiz;
1958 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
1959 F77_CONST_CHAR_ARG2 (" ", 1),
1960 0, 0, 0, 0, smlsiz
1961 F77_CHAR_ARG_LEN (6)
1962 F77_CHAR_ARG_LEN (1));
1963
1964 F77_INT mnthr;
1965 F77_FUNC (xilaenv, XILAENV) (6, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
1966 F77_CONST_CHAR_ARG2 (" ", 1),
1967 m, n, nrhs, -1, mnthr
1968 F77_CHAR_ARG_LEN (6)
1969 F77_CHAR_ARG_LEN (1));
1970
1971 // We compute the size of iwork because DGELSD in older versions
1972 // of LAPACK does not return it on a query call.
1973 double dminmn = static_cast<double> (minmn);
1974 double dsmlsizp1 = static_cast<double> (smlsiz+1);
1975 double tmp = octave::math::log2 (dminmn / dsmlsizp1);
1976
1977 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
1978 if (nlvl < 0)
1979 nlvl = 0;
1980
1981 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
1982 if (liwork < 1)
1983 liwork = 1;
1984 Array<F77_INT> iwork (dim_vector (liwork, 1));
1985 F77_INT *piwork = iwork.fortran_vec ();
1986
1987 F77_INT tmp_info = 0;
1988 F77_INT tmp_rank = 0;
1989
1990 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
1991 ps, rcon, tmp_rank, work.fortran_vec (),
1992 lwork, piwork, tmp_info));
1993
1994 info = tmp_info;
1995 rank = tmp_rank;
1996
1997 // The workspace query is broken in at least LAPACK 3.0.0
1998 // through 3.1.1 when n >= mnthr. The obtuse formula below
1999 // should provide sufficient workspace for DGELSD to operate
2000 // efficiently.
2001 if (n > m && n >= mnthr)
2002 {
2003 const F77_INT wlalsd
2004 = 9*m + 2*m*smlsiz + 8*m*nlvl + m*nrhs + (smlsiz+1)*(smlsiz+1);
2005
2006 F77_INT addend = m;
2007
2008 if (2*m-4 > addend)
2009 addend = 2*m-4;
2010
2011 if (nrhs > addend)
2012 addend = nrhs;
2013
2014 if (n-3*m > addend)
2015 addend = n-3*m;
2016
2017 if (wlalsd > addend)
2018 addend = wlalsd;
2019
2020 const F77_INT lworkaround = 4*m + m*m + addend;
2021
2022 if (work(0) < lworkaround)
2023 work(0) = lworkaround;
2024 }
2025 else if (m >= n)
2026 {
2027 F77_INT lworkaround
2028 = 12*n + 2*n*smlsiz + 8*n*nlvl + n*nrhs + (smlsiz+1)*(smlsiz+1);
2029
2030 if (work(0) < lworkaround)
2031 work(0) = lworkaround;
2032 }
2033
2034 lwork = static_cast<F77_INT> (work(0));
2035 work.resize (dim_vector (lwork, 1));
2036
2037 double anorm = norm1 (*this);
2038
2039 if (octave::math::isinf (anorm))
2040 {
2041 rcon = 0.0;
2042 retval = Matrix (n, b_nc, 0.0);
2043 }
2044 else if (octave::math::isnan (anorm))
2045 {
2046 rcon = octave::numeric_limits<double>::NaN ();
2047 retval = Matrix (n, b_nc, octave::numeric_limits<double>::NaN ());
2048 }
2049 else
2050 {
2051 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,
2052 maxmn, ps, rcon, tmp_rank,
2053 work.fortran_vec (), lwork,
2054 piwork, tmp_info));
2055
2056 info = tmp_info;
2057 rank = tmp_rank;
2058
2059 if (s.elem (0) == 0.0)
2060 rcon = 0.0;
2061 else
2062 rcon = s.elem (minmn - 1) / s.elem (0);
2063
2064 retval.resize (n, nrhs);
2065 }
2066 }
2067
2068 return retval;
2069 }
2070
2071 ComplexMatrix
lssolve(const ComplexMatrix & b) const2072 Matrix::lssolve (const ComplexMatrix& b) const
2073 {
2074 ComplexMatrix tmp (*this);
2075 octave_idx_type info;
2076 octave_idx_type rank;
2077 double rcon;
2078 return tmp.lssolve (b, info, rank, rcon);
2079 }
2080
2081 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info) const2082 Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info) const
2083 {
2084 ComplexMatrix tmp (*this);
2085 octave_idx_type rank;
2086 double rcon;
2087 return tmp.lssolve (b, info, rank, rcon);
2088 }
2089
2090 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info,octave_idx_type & rank) const2091 Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2092 octave_idx_type& rank) const
2093 {
2094 ComplexMatrix tmp (*this);
2095 double rcon;
2096 return tmp.lssolve (b, info, rank, rcon);
2097 }
2098
2099 ComplexMatrix
lssolve(const ComplexMatrix & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2100 Matrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
2101 octave_idx_type& rank, double& rcon) const
2102 {
2103 ComplexMatrix tmp (*this);
2104 return tmp.lssolve (b, info, rank, rcon);
2105 }
2106
2107 ColumnVector
lssolve(const ColumnVector & b) const2108 Matrix::lssolve (const ColumnVector& b) const
2109 {
2110 octave_idx_type info;
2111 octave_idx_type rank;
2112 double rcon;
2113 return lssolve (b, info, rank, rcon);
2114 }
2115
2116 ColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info) const2117 Matrix::lssolve (const ColumnVector& b, octave_idx_type& info) const
2118 {
2119 octave_idx_type rank;
2120 double rcon;
2121 return lssolve (b, info, rank, rcon);
2122 }
2123
2124 ColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info,octave_idx_type & rank) const2125 Matrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2126 octave_idx_type& rank) const
2127 {
2128 double rcon;
2129 return lssolve (b, info, rank, rcon);
2130 }
2131
2132 ColumnVector
lssolve(const ColumnVector & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2133 Matrix::lssolve (const ColumnVector& b, octave_idx_type& info,
2134 octave_idx_type& rank, double& rcon) const
2135 {
2136 ColumnVector retval;
2137
2138 F77_INT nrhs = 1;
2139
2140 F77_INT m = octave::to_f77_int (rows ());
2141 F77_INT n = octave::to_f77_int (cols ());
2142
2143 F77_INT b_nel = octave::to_f77_int (b.numel ());
2144
2145 if (m != b_nel)
2146 (*current_liboctave_error_handler)
2147 ("matrix dimension mismatch solution of linear equations");
2148
2149 if (m == 0 || n == 0)
2150 retval = ColumnVector (n, 0.0);
2151 else
2152 {
2153 volatile F77_INT minmn = (m < n ? m : n);
2154 F77_INT maxmn = (m > n ? m : n);
2155 rcon = -1.0;
2156
2157 if (m != n)
2158 {
2159 retval = ColumnVector (maxmn, 0.0);
2160
2161 for (F77_INT i = 0; i < m; i++)
2162 retval.elem (i) = b.elem (i);
2163 }
2164 else
2165 retval = b;
2166
2167 Matrix atmp = *this;
2168 double *tmp_data = atmp.fortran_vec ();
2169
2170 double *pretval = retval.fortran_vec ();
2171 Array<double> s (dim_vector (minmn, 1));
2172 double *ps = s.fortran_vec ();
2173
2174 // Ask DGELSD what the dimension of WORK should be.
2175 F77_INT lwork = -1;
2176
2177 Array<double> work (dim_vector (1, 1));
2178
2179 F77_INT smlsiz;
2180 F77_FUNC (xilaenv, XILAENV) (9, F77_CONST_CHAR_ARG2 ("DGELSD", 6),
2181 F77_CONST_CHAR_ARG2 (" ", 1),
2182 0, 0, 0, 0, smlsiz
2183 F77_CHAR_ARG_LEN (6)
2184 F77_CHAR_ARG_LEN (1));
2185
2186 // We compute the size of iwork because DGELSD in older versions
2187 // of LAPACK does not return it on a query call.
2188 double dminmn = static_cast<double> (minmn);
2189 double dsmlsizp1 = static_cast<double> (smlsiz+1);
2190 double tmp = octave::math::log2 (dminmn / dsmlsizp1);
2191
2192 F77_INT nlvl = static_cast<F77_INT> (tmp) + 1;
2193 if (nlvl < 0)
2194 nlvl = 0;
2195
2196 F77_INT liwork = 3 * minmn * nlvl + 11 * minmn;
2197 if (liwork < 1)
2198 liwork = 1;
2199 Array<F77_INT> iwork (dim_vector (liwork, 1));
2200 F77_INT *piwork = iwork.fortran_vec ();
2201
2202 F77_INT tmp_info = 0;
2203 F77_INT tmp_rank = 0;
2204
2205 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval, maxmn,
2206 ps, rcon, tmp_rank, work.fortran_vec (),
2207 lwork, piwork, tmp_info));
2208
2209 info = tmp_info;
2210 rank = tmp_rank;
2211
2212 lwork = static_cast<F77_INT> (work(0));
2213 work.resize (dim_vector (lwork, 1));
2214
2215 F77_XFCN (dgelsd, DGELSD, (m, n, nrhs, tmp_data, m, pretval,
2216 maxmn, ps, rcon, tmp_rank,
2217 work.fortran_vec (), lwork,
2218 piwork, tmp_info));
2219
2220 info = tmp_info;
2221 rank = tmp_rank;
2222
2223 if (rank < minmn)
2224 {
2225 if (s.elem (0) == 0.0)
2226 rcon = 0.0;
2227 else
2228 rcon = s.elem (minmn - 1) / s.elem (0);
2229 }
2230
2231 retval.resize (n);
2232 }
2233
2234 return retval;
2235 }
2236
2237 ComplexColumnVector
lssolve(const ComplexColumnVector & b) const2238 Matrix::lssolve (const ComplexColumnVector& b) const
2239 {
2240 ComplexMatrix tmp (*this);
2241 octave_idx_type info;
2242 octave_idx_type rank;
2243 double rcon;
2244 return tmp.lssolve (b, info, rank, rcon);
2245 }
2246
2247 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info) const2248 Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info) const
2249 {
2250 ComplexMatrix tmp (*this);
2251 octave_idx_type rank;
2252 double rcon;
2253 return tmp.lssolve (b, info, rank, rcon);
2254 }
2255
2256 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info,octave_idx_type & rank) const2257 Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2258 octave_idx_type& rank) const
2259 {
2260 ComplexMatrix tmp (*this);
2261 double rcon;
2262 return tmp.lssolve (b, info, rank, rcon);
2263 }
2264
2265 ComplexColumnVector
lssolve(const ComplexColumnVector & b,octave_idx_type & info,octave_idx_type & rank,double & rcon) const2266 Matrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
2267 octave_idx_type& rank, double& rcon) const
2268 {
2269 ComplexMatrix tmp (*this);
2270 return tmp.lssolve (b, info, rank, rcon);
2271 }
2272
2273 Matrix&
operator +=(const DiagMatrix & a)2274 Matrix::operator += (const DiagMatrix& a)
2275 {
2276 octave_idx_type nr = rows ();
2277 octave_idx_type nc = cols ();
2278
2279 octave_idx_type a_nr = a.rows ();
2280 octave_idx_type a_nc = a.cols ();
2281
2282 if (nr != a_nr || nc != a_nc)
2283 octave::err_nonconformant ("operator +=", nr, nc, a_nr, a_nc);
2284
2285 for (octave_idx_type i = 0; i < a.length (); i++)
2286 elem (i, i) += a.elem (i, i);
2287
2288 return *this;
2289 }
2290
2291 Matrix&
operator -=(const DiagMatrix & a)2292 Matrix::operator -= (const DiagMatrix& a)
2293 {
2294 octave_idx_type nr = rows ();
2295 octave_idx_type nc = cols ();
2296
2297 octave_idx_type a_nr = a.rows ();
2298 octave_idx_type a_nc = a.cols ();
2299
2300 if (nr != a_nr || nc != a_nc)
2301 octave::err_nonconformant ("operator -=", nr, nc, a_nr, a_nc);
2302
2303 for (octave_idx_type i = 0; i < a.length (); i++)
2304 elem (i, i) -= a.elem (i, i);
2305
2306 return *this;
2307 }
2308
2309 // unary operations
2310
2311 // column vector by row vector -> matrix operations
2312
2313 Matrix
operator *(const ColumnVector & v,const RowVector & a)2314 operator * (const ColumnVector& v, const RowVector& a)
2315 {
2316 Matrix retval;
2317
2318 F77_INT len = octave::to_f77_int (v.numel ());
2319
2320 if (len != 0)
2321 {
2322 F77_INT a_len = octave::to_f77_int (a.numel ());
2323
2324 retval = Matrix (len, a_len);
2325 double *c = retval.fortran_vec ();
2326
2327 F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 ("N", 1),
2328 F77_CONST_CHAR_ARG2 ("N", 1),
2329 len, a_len, 1, 1.0, v.data (), len,
2330 a.data (), 1, 0.0, c, len
2331 F77_CHAR_ARG_LEN (1)
2332 F77_CHAR_ARG_LEN (1)));
2333 }
2334
2335 return retval;
2336 }
2337
2338 // other operations.
2339
2340 // FIXME: Do these really belong here? Maybe they should be in a base class?
2341
2342 boolMatrix
all(int dim) const2343 Matrix::all (int dim) const
2344 {
2345 return NDArray::all (dim);
2346 }
2347
2348 boolMatrix
any(int dim) const2349 Matrix::any (int dim) const
2350 {
2351 return NDArray::any (dim);
2352 }
2353
2354 Matrix
cumprod(int dim) const2355 Matrix::cumprod (int dim) const
2356 {
2357 return NDArray::cumprod (dim);
2358 }
2359
2360 Matrix
cumsum(int dim) const2361 Matrix::cumsum (int dim) const
2362 {
2363 return NDArray::cumsum (dim);
2364 }
2365
2366 Matrix
prod(int dim) const2367 Matrix::prod (int dim) const
2368 {
2369 return NDArray::prod (dim);
2370 }
2371
2372 Matrix
sum(int dim) const2373 Matrix::sum (int dim) const
2374 {
2375 return NDArray::sum (dim);
2376 }
2377
2378 Matrix
sumsq(int dim) const2379 Matrix::sumsq (int dim) const
2380 {
2381 return NDArray::sumsq (dim);
2382 }
2383
2384 Matrix
abs(void) const2385 Matrix::abs (void) const
2386 {
2387 return NDArray::abs ();
2388 }
2389
2390 Matrix
diag(octave_idx_type k) const2391 Matrix::diag (octave_idx_type k) const
2392 {
2393 return NDArray::diag (k);
2394 }
2395
2396 DiagMatrix
diag(octave_idx_type m,octave_idx_type n) const2397 Matrix::diag (octave_idx_type m, octave_idx_type n) const
2398 {
2399 DiagMatrix retval;
2400
2401 octave_idx_type nr = rows ();
2402 octave_idx_type nc = cols ();
2403
2404 if (nr == 1 || nc == 1)
2405 retval = DiagMatrix (*this, m, n);
2406 else
2407 (*current_liboctave_error_handler) ("diag: expecting vector argument");
2408
2409 return retval;
2410 }
2411
2412 ColumnVector
row_min(void) const2413 Matrix::row_min (void) const
2414 {
2415 Array<octave_idx_type> dummy_idx;
2416 return row_min (dummy_idx);
2417 }
2418
2419 ColumnVector
row_min(Array<octave_idx_type> & idx_arg) const2420 Matrix::row_min (Array<octave_idx_type>& idx_arg) const
2421 {
2422 ColumnVector result;
2423
2424 octave_idx_type nr = rows ();
2425 octave_idx_type nc = cols ();
2426
2427 if (nr > 0 && nc > 0)
2428 {
2429 result.resize (nr);
2430 idx_arg.resize (dim_vector (nr, 1));
2431
2432 for (octave_idx_type i = 0; i < nr; i++)
2433 {
2434 octave_idx_type idx_j;
2435
2436 double tmp_min = octave::numeric_limits<double>::NaN ();
2437
2438 for (idx_j = 0; idx_j < nc; idx_j++)
2439 {
2440 tmp_min = elem (i, idx_j);
2441
2442 if (! octave::math::isnan (tmp_min))
2443 break;
2444 }
2445
2446 for (octave_idx_type j = idx_j+1; j < nc; j++)
2447 {
2448 double tmp = elem (i, j);
2449
2450 if (octave::math::isnan (tmp))
2451 continue;
2452 else if (tmp < tmp_min)
2453 {
2454 idx_j = j;
2455 tmp_min = tmp;
2456 }
2457 }
2458
2459 result.elem (i) = tmp_min;
2460 idx_arg.elem (i) = (octave::math::isnan (tmp_min) ? 0 : idx_j);
2461 }
2462 }
2463
2464 return result;
2465 }
2466
2467 ColumnVector
row_max(void) const2468 Matrix::row_max (void) const
2469 {
2470 Array<octave_idx_type> dummy_idx;
2471 return row_max (dummy_idx);
2472 }
2473
2474 ColumnVector
row_max(Array<octave_idx_type> & idx_arg) const2475 Matrix::row_max (Array<octave_idx_type>& idx_arg) const
2476 {
2477 ColumnVector result;
2478
2479 octave_idx_type nr = rows ();
2480 octave_idx_type nc = cols ();
2481
2482 if (nr > 0 && nc > 0)
2483 {
2484 result.resize (nr);
2485 idx_arg.resize (dim_vector (nr, 1));
2486
2487 for (octave_idx_type i = 0; i < nr; i++)
2488 {
2489 octave_idx_type idx_j;
2490
2491 double tmp_max = octave::numeric_limits<double>::NaN ();
2492
2493 for (idx_j = 0; idx_j < nc; idx_j++)
2494 {
2495 tmp_max = elem (i, idx_j);
2496
2497 if (! octave::math::isnan (tmp_max))
2498 break;
2499 }
2500
2501 for (octave_idx_type j = idx_j+1; j < nc; j++)
2502 {
2503 double tmp = elem (i, j);
2504
2505 if (octave::math::isnan (tmp))
2506 continue;
2507 else if (tmp > tmp_max)
2508 {
2509 idx_j = j;
2510 tmp_max = tmp;
2511 }
2512 }
2513
2514 result.elem (i) = tmp_max;
2515 idx_arg.elem (i) = (octave::math::isnan (tmp_max) ? 0 : idx_j);
2516 }
2517 }
2518
2519 return result;
2520 }
2521
2522 RowVector
column_min(void) const2523 Matrix::column_min (void) const
2524 {
2525 Array<octave_idx_type> dummy_idx;
2526 return column_min (dummy_idx);
2527 }
2528
2529 RowVector
column_min(Array<octave_idx_type> & idx_arg) const2530 Matrix::column_min (Array<octave_idx_type>& idx_arg) const
2531 {
2532 RowVector result;
2533
2534 octave_idx_type nr = rows ();
2535 octave_idx_type nc = cols ();
2536
2537 if (nr > 0 && nc > 0)
2538 {
2539 result.resize (nc);
2540 idx_arg.resize (dim_vector (1, nc));
2541
2542 for (octave_idx_type j = 0; j < nc; j++)
2543 {
2544 octave_idx_type idx_i;
2545
2546 double tmp_min = octave::numeric_limits<double>::NaN ();
2547
2548 for (idx_i = 0; idx_i < nr; idx_i++)
2549 {
2550 tmp_min = elem (idx_i, j);
2551
2552 if (! octave::math::isnan (tmp_min))
2553 break;
2554 }
2555
2556 for (octave_idx_type i = idx_i+1; i < nr; i++)
2557 {
2558 double tmp = elem (i, j);
2559
2560 if (octave::math::isnan (tmp))
2561 continue;
2562 else if (tmp < tmp_min)
2563 {
2564 idx_i = i;
2565 tmp_min = tmp;
2566 }
2567 }
2568
2569 result.elem (j) = tmp_min;
2570 idx_arg.elem (j) = (octave::math::isnan (tmp_min) ? 0 : idx_i);
2571 }
2572 }
2573
2574 return result;
2575 }
2576
2577 RowVector
column_max(void) const2578 Matrix::column_max (void) const
2579 {
2580 Array<octave_idx_type> dummy_idx;
2581 return column_max (dummy_idx);
2582 }
2583
2584 RowVector
column_max(Array<octave_idx_type> & idx_arg) const2585 Matrix::column_max (Array<octave_idx_type>& idx_arg) const
2586 {
2587 RowVector result;
2588
2589 octave_idx_type nr = rows ();
2590 octave_idx_type nc = cols ();
2591
2592 if (nr > 0 && nc > 0)
2593 {
2594 result.resize (nc);
2595 idx_arg.resize (dim_vector (1, nc));
2596
2597 for (octave_idx_type j = 0; j < nc; j++)
2598 {
2599 octave_idx_type idx_i;
2600
2601 double tmp_max = octave::numeric_limits<double>::NaN ();
2602
2603 for (idx_i = 0; idx_i < nr; idx_i++)
2604 {
2605 tmp_max = elem (idx_i, j);
2606
2607 if (! octave::math::isnan (tmp_max))
2608 break;
2609 }
2610
2611 for (octave_idx_type i = idx_i+1; i < nr; i++)
2612 {
2613 double tmp = elem (i, j);
2614
2615 if (octave::math::isnan (tmp))
2616 continue;
2617 else if (tmp > tmp_max)
2618 {
2619 idx_i = i;
2620 tmp_max = tmp;
2621 }
2622 }
2623
2624 result.elem (j) = tmp_max;
2625 idx_arg.elem (j) = (octave::math::isnan (tmp_max) ? 0 : idx_i);
2626 }
2627 }
2628
2629 return result;
2630 }
2631
2632 std::ostream&
operator <<(std::ostream & os,const Matrix & a)2633 operator << (std::ostream& os, const Matrix& a)
2634 {
2635 for (octave_idx_type i = 0; i < a.rows (); i++)
2636 {
2637 for (octave_idx_type j = 0; j < a.cols (); j++)
2638 {
2639 os << ' ';
2640 octave_write_double (os, a.elem (i, j));
2641 }
2642 os << "\n";
2643 }
2644 return os;
2645 }
2646
2647 std::istream&
operator >>(std::istream & is,Matrix & a)2648 operator >> (std::istream& is, Matrix& a)
2649 {
2650 octave_idx_type nr = a.rows ();
2651 octave_idx_type nc = a.cols ();
2652
2653 if (nr > 0 && nc > 0)
2654 {
2655 double tmp;
2656 for (octave_idx_type i = 0; i < nr; i++)
2657 for (octave_idx_type j = 0; j < nc; j++)
2658 {
2659 tmp = octave_read_value<double> (is);
2660 if (is)
2661 a.elem (i, j) = tmp;
2662 else
2663 return is;
2664 }
2665 }
2666
2667 return is;
2668 }
2669
2670 Matrix
Givens(double x,double y)2671 Givens (double x, double y)
2672 {
2673 double cc, s, temp_r;
2674
2675 F77_FUNC (dlartg, DLARTG) (x, y, cc, s, temp_r);
2676
2677 Matrix g (2, 2);
2678
2679 g.elem (0, 0) = cc;
2680 g.elem (1, 1) = cc;
2681 g.elem (0, 1) = s;
2682 g.elem (1, 0) = -s;
2683
2684 return g;
2685 }
2686
2687 Matrix
Sylvester(const Matrix & a,const Matrix & b,const Matrix & c)2688 Sylvester (const Matrix& a, const Matrix& b, const Matrix& c)
2689 {
2690 Matrix retval;
2691
2692 // FIXME: need to check that a, b, and c are all the same size.
2693
2694 // Compute Schur decompositions.
2695
2696 octave::math::schur<Matrix> as (a, "U");
2697 octave::math::schur<Matrix> bs (b, "U");
2698
2699 // Transform c to new coordinates.
2700
2701 Matrix ua = as.unitary_matrix ();
2702 Matrix sch_a = as.schur_matrix ();
2703
2704 Matrix ub = bs.unitary_matrix ();
2705 Matrix sch_b = bs.schur_matrix ();
2706
2707 Matrix cx = ua.transpose () * c * ub;
2708
2709 // Solve the sylvester equation, back-transform, and return the solution.
2710
2711 F77_INT a_nr = octave::to_f77_int (a.rows ());
2712 F77_INT b_nr = octave::to_f77_int (b.rows ());
2713
2714 double scale;
2715 F77_INT info;
2716
2717 double *pa = sch_a.fortran_vec ();
2718 double *pb = sch_b.fortran_vec ();
2719 double *px = cx.fortran_vec ();
2720
2721 F77_XFCN (dtrsyl, DTRSYL, (F77_CONST_CHAR_ARG2 ("N", 1),
2722 F77_CONST_CHAR_ARG2 ("N", 1),
2723 1, a_nr, b_nr, pa, a_nr, pb,
2724 b_nr, px, a_nr, scale, info
2725 F77_CHAR_ARG_LEN (1)
2726 F77_CHAR_ARG_LEN (1)));
2727
2728 // FIXME: check info?
2729
2730 retval = ua*cx*ub.transpose ();
2731
2732 return retval;
2733 }
2734
2735 // matrix by matrix -> matrix operations
2736
2737 /*
2738
2739 ## Simple Dot Product, Matrix-Vector and Matrix-Matrix Unit tests
2740 %!assert ([1 2 3] * [ 4 ; 5 ; 6], 32, 1e-14)
2741 %!assert ([1 2 ; 3 4 ] * [5 ; 6], [17 ; 39 ], 1e-14)
2742 %!assert ([1 2 ; 3 4 ] * [5 6 ; 7 8], [19 22; 43 50], 1e-14)
2743
2744 ## Test some simple identities
2745 %!shared M, cv, rv, Mt, rvt
2746 %! M = randn (10,10) + 100*eye (10,10);
2747 %! Mt = M';
2748 %! cv = randn (10,1);
2749 %! rv = randn (1,10);
2750 %! rvt = rv';
2751 %!assert ([M*cv,M*cv], M*[cv,cv], 2e-13)
2752 %!assert ([M'*cv,M'*cv], M'*[cv,cv], 2e-13)
2753 %!assert ([rv*M;rv*M], [rv;rv]*M, 2e-13)
2754 %!assert ([rv*M';rv*M'], [rv;rv]*M', 2e-13)
2755 %!assert (2*rv*cv, [rv,rv]*[cv;cv], 2e-13)
2756 %!assert (M'\cv, Mt\cv, 1e-14)
2757 %!assert (M'\rv', Mt\rvt, 1e-14)
2758
2759 */
2760
2761 static inline char
get_blas_trans_arg(bool trans)2762 get_blas_trans_arg (bool trans)
2763 {
2764 return trans ? 'T' : 'N';
2765 }
2766
2767 // the general GEMM operation
2768
2769 Matrix
xgemm(const Matrix & a,const Matrix & b,blas_trans_type transa,blas_trans_type transb)2770 xgemm (const Matrix& a, const Matrix& b,
2771 blas_trans_type transa, blas_trans_type transb)
2772 {
2773 Matrix retval;
2774
2775 bool tra = transa != blas_no_trans;
2776 bool trb = transb != blas_no_trans;
2777
2778 F77_INT a_nr = octave::to_f77_int (tra ? a.cols () : a.rows ());
2779 F77_INT a_nc = octave::to_f77_int (tra ? a.rows () : a.cols ());
2780
2781 F77_INT b_nr = octave::to_f77_int (trb ? b.cols () : b.rows ());
2782 F77_INT b_nc = octave::to_f77_int (trb ? b.rows () : b.cols ());
2783
2784 if (a_nc != b_nr)
2785 octave::err_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
2786
2787 if (a_nr == 0 || a_nc == 0 || b_nc == 0)
2788 retval = Matrix (a_nr, b_nc, 0.0);
2789 else if (a.data () == b.data () && a_nr == b_nc && tra != trb)
2790 {
2791 F77_INT lda = octave::to_f77_int (a.rows ());
2792
2793 retval = Matrix (a_nr, b_nc);
2794 double *c = retval.fortran_vec ();
2795
2796 const char ctra = get_blas_trans_arg (tra);
2797 F77_XFCN (dsyrk, DSYRK, (F77_CONST_CHAR_ARG2 ("U", 1),
2798 F77_CONST_CHAR_ARG2 (&ctra, 1),
2799 a_nr, a_nc, 1.0,
2800 a.data (), lda, 0.0, c, a_nr
2801 F77_CHAR_ARG_LEN (1)
2802 F77_CHAR_ARG_LEN (1)));
2803 for (int j = 0; j < a_nr; j++)
2804 for (int i = 0; i < j; i++)
2805 retval.xelem (j,i) = retval.xelem (i,j);
2806
2807 }
2808 else
2809 {
2810 F77_INT lda = octave::to_f77_int (a.rows ());
2811 F77_INT tda = octave::to_f77_int (a.cols ());
2812 F77_INT ldb = octave::to_f77_int (b.rows ());
2813 F77_INT tdb = octave::to_f77_int (b.cols ());
2814
2815 retval = Matrix (a_nr, b_nc);
2816 double *c = retval.fortran_vec ();
2817
2818 if (b_nc == 1)
2819 {
2820 if (a_nr == 1)
2821 F77_FUNC (xddot, XDDOT) (a_nc, a.data (), 1, b.data (), 1, *c);
2822 else
2823 {
2824 const char ctra = get_blas_trans_arg (tra);
2825 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 (&ctra, 1),
2826 lda, tda, 1.0, a.data (), lda,
2827 b.data (), 1, 0.0, c, 1
2828 F77_CHAR_ARG_LEN (1)));
2829 }
2830 }
2831 else if (a_nr == 1)
2832 {
2833 const char crevtrb = get_blas_trans_arg (! trb);
2834 F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 (&crevtrb, 1),
2835 ldb, tdb, 1.0, b.data (), ldb,
2836 a.data (), 1, 0.0, c, 1
2837 F77_CHAR_ARG_LEN (1)));
2838 }
2839 else
2840 {
2841 const char ctra = get_blas_trans_arg (tra);
2842 const char ctrb = get_blas_trans_arg (trb);
2843 F77_XFCN (dgemm, DGEMM, (F77_CONST_CHAR_ARG2 (&ctra, 1),
2844 F77_CONST_CHAR_ARG2 (&ctrb, 1),
2845 a_nr, b_nc, a_nc, 1.0, a.data (),
2846 lda, b.data (), ldb, 0.0, c, a_nr
2847 F77_CHAR_ARG_LEN (1)
2848 F77_CHAR_ARG_LEN (1)));
2849 }
2850 }
2851
2852 return retval;
2853 }
2854
2855 Matrix
operator *(const Matrix & a,const Matrix & b)2856 operator * (const Matrix& a, const Matrix& b)
2857 {
2858 return xgemm (a, b);
2859 }
2860
2861 // FIXME: it would be nice to share code among the min/max functions below.
2862
2863 #define EMPTY_RETURN_CHECK(T) \
2864 if (nr == 0 || nc == 0) \
2865 return T (nr, nc);
2866
2867 Matrix
min(double d,const Matrix & m)2868 min (double d, const Matrix& m)
2869 {
2870 octave_idx_type nr = m.rows ();
2871 octave_idx_type nc = m.columns ();
2872
2873 EMPTY_RETURN_CHECK (Matrix);
2874
2875 Matrix result (nr, nc);
2876
2877 for (octave_idx_type j = 0; j < nc; j++)
2878 for (octave_idx_type i = 0; i < nr; i++)
2879 {
2880 octave_quit ();
2881 result(i, j) = octave::math::min (d, m(i, j));
2882 }
2883
2884 return result;
2885 }
2886
2887 Matrix
min(const Matrix & m,double d)2888 min (const Matrix& m, double d)
2889 {
2890 octave_idx_type nr = m.rows ();
2891 octave_idx_type nc = m.columns ();
2892
2893 EMPTY_RETURN_CHECK (Matrix);
2894
2895 Matrix result (nr, nc);
2896
2897 for (octave_idx_type j = 0; j < nc; j++)
2898 for (octave_idx_type i = 0; i < nr; i++)
2899 {
2900 octave_quit ();
2901 result(i, j) = octave::math::min (m(i, j), d);
2902 }
2903
2904 return result;
2905 }
2906
2907 Matrix
min(const Matrix & a,const Matrix & b)2908 min (const Matrix& a, const Matrix& b)
2909 {
2910 octave_idx_type nr = a.rows ();
2911 octave_idx_type nc = a.columns ();
2912
2913 if (nr != b.rows () || nc != b.columns ())
2914 (*current_liboctave_error_handler)
2915 ("two-arg min requires same size arguments");
2916
2917 EMPTY_RETURN_CHECK (Matrix);
2918
2919 Matrix result (nr, nc);
2920
2921 for (octave_idx_type j = 0; j < nc; j++)
2922 for (octave_idx_type i = 0; i < nr; i++)
2923 {
2924 octave_quit ();
2925 result(i, j) = octave::math::min (a(i, j), b(i, j));
2926 }
2927
2928 return result;
2929 }
2930
2931 Matrix
max(double d,const Matrix & m)2932 max (double d, const Matrix& m)
2933 {
2934 octave_idx_type nr = m.rows ();
2935 octave_idx_type nc = m.columns ();
2936
2937 EMPTY_RETURN_CHECK (Matrix);
2938
2939 Matrix result (nr, nc);
2940
2941 for (octave_idx_type j = 0; j < nc; j++)
2942 for (octave_idx_type i = 0; i < nr; i++)
2943 {
2944 octave_quit ();
2945 result(i, j) = octave::math::max (d, m(i, j));
2946 }
2947
2948 return result;
2949 }
2950
2951 Matrix
max(const Matrix & m,double d)2952 max (const Matrix& m, double d)
2953 {
2954 octave_idx_type nr = m.rows ();
2955 octave_idx_type nc = m.columns ();
2956
2957 EMPTY_RETURN_CHECK (Matrix);
2958
2959 Matrix result (nr, nc);
2960
2961 for (octave_idx_type j = 0; j < nc; j++)
2962 for (octave_idx_type i = 0; i < nr; i++)
2963 {
2964 octave_quit ();
2965 result(i, j) = octave::math::max (m(i, j), d);
2966 }
2967
2968 return result;
2969 }
2970
2971 Matrix
max(const Matrix & a,const Matrix & b)2972 max (const Matrix& a, const Matrix& b)
2973 {
2974 octave_idx_type nr = a.rows ();
2975 octave_idx_type nc = a.columns ();
2976
2977 if (nr != b.rows () || nc != b.columns ())
2978 (*current_liboctave_error_handler)
2979 ("two-arg max requires same size arguments");
2980
2981 EMPTY_RETURN_CHECK (Matrix);
2982
2983 Matrix result (nr, nc);
2984
2985 for (octave_idx_type j = 0; j < nc; j++)
2986 for (octave_idx_type i = 0; i < nr; i++)
2987 {
2988 octave_quit ();
2989 result(i, j) = octave::math::max (a(i, j), b(i, j));
2990 }
2991
2992 return result;
2993 }
2994
linspace(const ColumnVector & x1,const ColumnVector & x2,octave_idx_type n)2995 Matrix linspace (const ColumnVector& x1,
2996 const ColumnVector& x2,
2997 octave_idx_type n)
2998
2999 {
3000 octave_idx_type m = x1.numel ();
3001
3002 if (x2.numel () != m)
3003 (*current_liboctave_error_handler)
3004 ("linspace: vectors must be of equal length");
3005
3006 Matrix retval;
3007
3008 if (n < 1)
3009 {
3010 retval.clear (m, 0);
3011 return retval;
3012 }
3013
3014 retval.clear (m, n);
3015 for (octave_idx_type i = 0; i < m; i++)
3016 retval.xelem (i, 0) = x1(i);
3017
3018 // The last column is unused so temporarily store delta there
3019 double *delta = &retval.xelem (0, n-1);
3020 for (octave_idx_type i = 0; i < m; i++)
3021 delta[i] = (x1(i) == x2(i)) ? 0 : (x2(i) - x1(i)) / (n - 1);
3022
3023 for (octave_idx_type j = 1; j < n-1; j++)
3024 for (octave_idx_type i = 0; i < m; i++)
3025 retval.xelem (i, j) = x1(i) + j*delta[i];
3026
3027 for (octave_idx_type i = 0; i < m; i++)
3028 retval.xelem (i, n-1) = x2(i);
3029
3030 return retval;
3031 }
3032
3033 MS_CMP_OPS (Matrix, double)
3034 MS_BOOL_OPS (Matrix, double)
3035
3036 SM_CMP_OPS (double, Matrix)
3037 SM_BOOL_OPS (double, Matrix)
3038
3039 MM_CMP_OPS (Matrix, Matrix)
3040 MM_BOOL_OPS (Matrix, Matrix)
3041