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 (octave_CMatrix_h)
27 #define octave_CMatrix_h 1
28 
29 #include "octave-config.h"
30 
31 #include "CNDArray.h"
32 #include "DET.h"
33 #include "MArray.h"
34 #include "MDiagArray2.h"
35 #include "MatrixType.h"
36 #include "mx-defs.h"
37 #include "mx-op-decl.h"
38 #include "oct-cmplx.h"
39 
40 class
41 OCTAVE_API
42 ComplexMatrix : public ComplexNDArray
43 {
44 public:
45 
46   typedef ComplexColumnVector column_vector_type;
47   typedef ComplexRowVector row_vector_type;
48 
49   typedef ColumnVector real_column_vector_type;
50   typedef RowVector real_row_vector_type;
51 
52   typedef Matrix real_matrix_type;
53   typedef ComplexMatrix complex_matrix_type;
54 
55   typedef DiagMatrix real_diag_matrix_type;
56   typedef ComplexDiagMatrix complex_diag_matrix_type;
57 
58   typedef double real_elt_type;
59   typedef Complex complex_elt_type;
60 
61   typedef void (*solve_singularity_handler) (double rcon);
62 
63   ComplexMatrix (void) = default;
64 
65   ComplexMatrix (const ComplexMatrix& a) = default;
66 
67   ComplexMatrix& operator = (const ComplexMatrix& a) = default;
68 
69   ~ComplexMatrix (void) = default;
70 
ComplexMatrix(octave_idx_type r,octave_idx_type c)71   ComplexMatrix (octave_idx_type r, octave_idx_type c)
72     : ComplexNDArray (dim_vector (r, c)) { }
73 
ComplexMatrix(octave_idx_type r,octave_idx_type c,const Complex & val)74   ComplexMatrix (octave_idx_type r, octave_idx_type c, const Complex& val)
75     : ComplexNDArray (dim_vector (r, c), val) { }
76 
ComplexMatrix(const dim_vector & dv)77   ComplexMatrix (const dim_vector& dv) : ComplexNDArray (dv.redim (2)) { }
78 
ComplexMatrix(const dim_vector & dv,const Complex & val)79   ComplexMatrix (const dim_vector& dv, const Complex& val)
80     : ComplexNDArray (dv.redim (2), val) { }
81 
82   template <typename U>
ComplexMatrix(const MArray<U> & a)83   ComplexMatrix (const MArray<U>& a) : ComplexNDArray (a.as_matrix ()) { }
84 
85   template <typename U>
ComplexMatrix(const Array<U> & a)86   ComplexMatrix (const Array<U>& a) : ComplexNDArray (a.as_matrix ()) { }
87 
88   ComplexMatrix (const Matrix& re, const Matrix& im);
89 
90   explicit ComplexMatrix (const Matrix& a);
91 
92   explicit ComplexMatrix (const RowVector& rv);
93 
94   explicit ComplexMatrix (const ColumnVector& cv);
95 
96   explicit ComplexMatrix (const DiagMatrix& a);
97 
98   explicit ComplexMatrix (const MDiagArray2<double>& a);
99 
100   explicit ComplexMatrix (const DiagArray2<double>& a);
101 
102   explicit ComplexMatrix (const ComplexRowVector& rv);
103 
104   explicit ComplexMatrix (const ComplexColumnVector& cv);
105 
106   explicit ComplexMatrix (const ComplexDiagMatrix& a);
107 
108   explicit ComplexMatrix (const MDiagArray2<Complex>& a);
109 
110   explicit ComplexMatrix (const DiagArray2<Complex>& a);
111 
112   explicit ComplexMatrix (const boolMatrix& a);
113 
114   explicit ComplexMatrix (const charMatrix& a);
115 
116   bool operator == (const ComplexMatrix& a) const;
117   bool operator != (const ComplexMatrix& a) const;
118 
119   bool ishermitian (void) const;
120 
121   // destructive insert/delete/reorder operations
122 
123   ComplexMatrix& insert (const Matrix& a, octave_idx_type r, octave_idx_type c);
124   ComplexMatrix& insert (const RowVector& a,
125                          octave_idx_type r, octave_idx_type c);
126   ComplexMatrix& insert (const ColumnVector& a,
127                          octave_idx_type r, octave_idx_type c);
128   ComplexMatrix& insert (const DiagMatrix& a,
129                          octave_idx_type r, octave_idx_type c);
130 
131   ComplexMatrix& insert (const ComplexMatrix& a,
132                          octave_idx_type r, octave_idx_type c);
133   ComplexMatrix& insert (const ComplexRowVector& a,
134                          octave_idx_type r, octave_idx_type c);
135   ComplexMatrix& insert (const ComplexColumnVector& a,
136                          octave_idx_type r, octave_idx_type c);
137   ComplexMatrix& insert (const ComplexDiagMatrix& a,
138                          octave_idx_type r, octave_idx_type c);
139 
140   ComplexMatrix& fill (double val);
141   ComplexMatrix& fill (const Complex& val);
142   ComplexMatrix& fill (double val,
143                        octave_idx_type r1, octave_idx_type c1,
144                        octave_idx_type r2, octave_idx_type c2);
145   ComplexMatrix& fill (const Complex& val,
146                        octave_idx_type r1, octave_idx_type c1,
147                        octave_idx_type r2, octave_idx_type c2);
148 
149   ComplexMatrix append (const Matrix& a) const;
150   ComplexMatrix append (const RowVector& a) const;
151   ComplexMatrix append (const ColumnVector& a) const;
152   ComplexMatrix append (const DiagMatrix& a) const;
153 
154   ComplexMatrix append (const ComplexMatrix& a) const;
155   ComplexMatrix append (const ComplexRowVector& a) const;
156   ComplexMatrix append (const ComplexColumnVector& a) const;
157   ComplexMatrix append (const ComplexDiagMatrix& a) const;
158 
159   ComplexMatrix stack (const Matrix& a) const;
160   ComplexMatrix stack (const RowVector& a) const;
161   ComplexMatrix stack (const ColumnVector& a) const;
162   ComplexMatrix stack (const DiagMatrix& a) const;
163 
164   ComplexMatrix stack (const ComplexMatrix& a) const;
165   ComplexMatrix stack (const ComplexRowVector& a) const;
166   ComplexMatrix stack (const ComplexColumnVector& a) const;
167   ComplexMatrix stack (const ComplexDiagMatrix& a) const;
168 
hermitian(void)169   ComplexMatrix hermitian (void) const
170   { return MArray<Complex>::hermitian (std::conj); }
transpose(void)171   ComplexMatrix transpose (void) const
172   { return MArray<Complex>::transpose (); }
173 
174   friend OCTAVE_API ComplexMatrix conj (const ComplexMatrix& a);
175 
176   // resize is the destructive equivalent for this one
177 
178   ComplexMatrix extract (octave_idx_type r1, octave_idx_type c1,
179                          octave_idx_type r2, octave_idx_type c2) const;
180 
181   ComplexMatrix extract_n (octave_idx_type r1, octave_idx_type c1,
182                            octave_idx_type nr, octave_idx_type nc) const;
183 
184   // extract row or column i.
185 
186   ComplexRowVector row (octave_idx_type i) const;
187 
188   ComplexColumnVector column (octave_idx_type i) const;
189 
190   void resize (octave_idx_type nr, octave_idx_type nc,
191                const Complex& rfv = Complex (0))
192   {
193     MArray<Complex>::resize (dim_vector (nr, nc), rfv);
194   }
195 
196 private:
197   ComplexMatrix tinverse (MatrixType& mattype, octave_idx_type& info,
198                           double& rcon, bool force, bool calc_cond) const;
199 
200   ComplexMatrix finverse (MatrixType& mattype, octave_idx_type& info,
201                           double& rcon, bool force, bool calc_cond) const;
202 
203 public:
204   ComplexMatrix inverse (void) const;
205   ComplexMatrix inverse (octave_idx_type& info) const;
206   ComplexMatrix inverse (octave_idx_type& info, double& rcon,
207                          bool force = false, bool calc_cond = true) const;
208 
209   ComplexMatrix inverse (MatrixType& mattype) const;
210   ComplexMatrix inverse (MatrixType& mattype, octave_idx_type& info) const;
211   ComplexMatrix inverse (MatrixType& mattype, octave_idx_type& info,
212                          double& rcon, bool force = false,
213                          bool calc_cond = true) const;
214 
215   ComplexMatrix pseudo_inverse (double tol = 0.0) const;
216 
217   ComplexMatrix fourier (void) const;
218   ComplexMatrix ifourier (void) const;
219 
220   ComplexMatrix fourier2d (void) const;
221   ComplexMatrix ifourier2d (void) const;
222 
223   ComplexDET determinant (void) const;
224   ComplexDET determinant (octave_idx_type& info) const;
225   ComplexDET determinant (octave_idx_type& info, double& rcon,
226                           bool calc_cond = true) const;
227   ComplexDET determinant (MatrixType& mattype, octave_idx_type& info,
228                           double& rcon, bool calc_cond = true) const;
229 
230   double rcond (void) const;
231   double rcond (MatrixType& mattype) const;
232 
233 private:
234   // Upper triangular matrix solvers
235   ComplexMatrix utsolve (MatrixType& mattype, const ComplexMatrix& b,
236                          octave_idx_type& info, double& rcon,
237                          solve_singularity_handler sing_handler,
238                          bool calc_cond = false,
239                          blas_trans_type transt = blas_no_trans) const;
240 
241   // Lower triangular matrix solvers
242   ComplexMatrix ltsolve (MatrixType& mattype, const ComplexMatrix& b,
243                          octave_idx_type& info, double& rcon,
244                          solve_singularity_handler sing_handler,
245                          bool calc_cond = false,
246                          blas_trans_type transt = blas_no_trans) const;
247 
248   // Full matrix solvers (umfpack/cholesky)
249   ComplexMatrix fsolve (MatrixType& mattype, const ComplexMatrix& b,
250                         octave_idx_type& info, double& rcon,
251                         solve_singularity_handler sing_handler,
252                         bool calc_cond = false) const;
253 
254 public:
255   // Generic interface to solver with no probing of type
256   ComplexMatrix solve (MatrixType& mattype, const Matrix& b) const;
257   ComplexMatrix solve (MatrixType& mattype, const Matrix& b,
258                        octave_idx_type& info) const;
259   ComplexMatrix solve (MatrixType& mattype, const Matrix& b,
260                        octave_idx_type& info, double& rcon) const;
261   ComplexMatrix solve (MatrixType& mattype, const Matrix& b,
262                        octave_idx_type& info, double& rcon,
263                        solve_singularity_handler sing_handler,
264                        bool singular_fallback = true,
265                        blas_trans_type transt = blas_no_trans) const;
266 
267   ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b) const;
268   ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b,
269                        octave_idx_type& info) const;
270   ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b,
271                        octave_idx_type& info, double& rcon) const;
272   ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b,
273                        octave_idx_type& info, double& rcon,
274                        solve_singularity_handler sing_handler,
275                        bool singular_fallback = true,
276                        blas_trans_type transt = blas_no_trans) const;
277 
278   ComplexColumnVector solve (MatrixType& mattype, const ColumnVector& b) const;
279   ComplexColumnVector solve (MatrixType& mattype, const ColumnVector& b,
280                              octave_idx_type& info) const;
281   ComplexColumnVector solve (MatrixType& mattype, const ColumnVector& b,
282                              octave_idx_type& info, double& rcon) const;
283   ComplexColumnVector solve (MatrixType& mattype, const ColumnVector& b,
284                              octave_idx_type& info, double& rcon,
285                              solve_singularity_handler sing_handler,
286                              blas_trans_type transt = blas_no_trans) const;
287 
288   ComplexColumnVector solve (MatrixType& mattype,
289                              const ComplexColumnVector& b) const;
290   ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b,
291                              octave_idx_type& info) const;
292   ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b,
293                              octave_idx_type& info, double& rcon) const;
294   ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b,
295                              octave_idx_type& info, double& rcon,
296                              solve_singularity_handler sing_handler,
297                              blas_trans_type transt = blas_no_trans) const;
298 
299   // Generic interface to solver with probing of type
300   ComplexMatrix solve (const Matrix& b) const;
301   ComplexMatrix solve (const Matrix& b, octave_idx_type& info) const;
302   ComplexMatrix solve (const Matrix& b, octave_idx_type& info,
303                        double& rcon) const;
304   ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcon,
305                        solve_singularity_handler sing_handler,
306                        blas_trans_type transt = blas_no_trans) const;
307 
308   ComplexMatrix solve (const ComplexMatrix& b) const;
309   ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const;
310   ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info,
311                        double& rcon) const;
312   ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info,
313                        double& rcon,
314                        solve_singularity_handler sing_handler,
315                        blas_trans_type transt = blas_no_trans) const;
316 
317   ComplexColumnVector solve (const ColumnVector& b) const;
318   ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info) const;
319   ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info,
320                              double& rcon) const;
321   ComplexColumnVector solve (const ColumnVector& b, octave_idx_type& info,
322                              double& rcon,
323                              solve_singularity_handler sing_handler,
324                              blas_trans_type transt = blas_no_trans) const;
325 
326   ComplexColumnVector solve (const ComplexColumnVector& b) const;
327   ComplexColumnVector solve (const ComplexColumnVector& b,
328                              octave_idx_type& info) const;
329   ComplexColumnVector solve (const ComplexColumnVector& b,
330                              octave_idx_type& info,
331                              double& rcon) const;
332   ComplexColumnVector solve (const ComplexColumnVector& b,
333                              octave_idx_type& info,
334                              double& rcon,
335                              solve_singularity_handler sing_handler,
336                              blas_trans_type transt = blas_no_trans) const;
337 
338   ComplexMatrix lssolve (const Matrix& b) const;
339   ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info) const;
340   ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info,
341                          octave_idx_type& rank) const;
342   ComplexMatrix lssolve (const Matrix& b, octave_idx_type& info,
343                          octave_idx_type& rank, double& rcon) const;
344 
345   ComplexMatrix lssolve (const ComplexMatrix& b) const;
346   ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info) const;
347   ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info,
348                          octave_idx_type& rank) const;
349   ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info,
350                          octave_idx_type& rank, double& rcon) const;
351 
352   ComplexColumnVector lssolve (const ColumnVector& b) const;
353   ComplexColumnVector lssolve (const ColumnVector& b,
354                                octave_idx_type& info) const;
355   ComplexColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
356                                octave_idx_type& rank) const;
357   ComplexColumnVector lssolve (const ColumnVector& b, octave_idx_type& info,
358                                octave_idx_type& rank, double& rcon) const;
359 
360   ComplexColumnVector lssolve (const ComplexColumnVector& b) const;
361   ComplexColumnVector lssolve (const ComplexColumnVector& b,
362                                octave_idx_type& info) const;
363   ComplexColumnVector lssolve (const ComplexColumnVector& b,
364                                octave_idx_type& info,
365                                octave_idx_type& rank) const;
366   ComplexColumnVector lssolve (const ComplexColumnVector& b,
367                                octave_idx_type& info,
368                                octave_idx_type& rank, double& rcon) const;
369 
370   // matrix by diagonal matrix -> matrix operations
371 
372   ComplexMatrix& operator += (const DiagMatrix& a);
373   ComplexMatrix& operator -= (const DiagMatrix& a);
374 
375   ComplexMatrix& operator += (const ComplexDiagMatrix& a);
376   ComplexMatrix& operator -= (const ComplexDiagMatrix& a);
377 
378   // matrix by matrix -> matrix operations
379 
380   ComplexMatrix& operator += (const Matrix& a);
381   ComplexMatrix& operator -= (const Matrix& a);
382 
383   // other operations
384 
385   boolMatrix all (int dim = -1) const;
386   boolMatrix any (int dim = -1) const;
387 
388   ComplexMatrix cumprod (int dim = -1) const;
389   ComplexMatrix cumsum (int dim = -1) const;
390   ComplexMatrix prod (int dim = -1) const;
391   ComplexMatrix sum (int dim = -1) const;
392   ComplexMatrix sumsq (int dim = -1) const;
393   Matrix abs (void) const;
394 
395   ComplexMatrix diag (octave_idx_type k = 0) const;
396 
397   ComplexDiagMatrix diag (octave_idx_type m, octave_idx_type n) const;
398 
399   bool row_is_real_only (octave_idx_type) const;
400   bool column_is_real_only (octave_idx_type) const;
401 
402   ComplexColumnVector row_min (void) const;
403   ComplexColumnVector row_max (void) const;
404 
405   ComplexColumnVector row_min (Array<octave_idx_type>& index) const;
406   ComplexColumnVector row_max (Array<octave_idx_type>& index) const;
407 
408   ComplexRowVector column_min (void) const;
409   ComplexRowVector column_max (void) const;
410 
411   ComplexRowVector column_min (Array<octave_idx_type>& index) const;
412   ComplexRowVector column_max (Array<octave_idx_type>& index) const;
413 
414   // i/o
415 
416   friend OCTAVE_API std::ostream& operator << (std::ostream& os,
417                                                const ComplexMatrix& a);
418   friend OCTAVE_API std::istream& operator >> (std::istream& is,
419                                                ComplexMatrix& a);
420 };
421 
422 extern OCTAVE_API ComplexMatrix conj (const ComplexMatrix& a);
423 
424 // column vector by row vector -> matrix operations
425 
426 extern OCTAVE_API ComplexMatrix
427 operator * (const ColumnVector& a, const ComplexRowVector& b);
428 
429 extern OCTAVE_API ComplexMatrix
430 operator * (const ComplexColumnVector& a, const RowVector& b);
431 
432 extern OCTAVE_API ComplexMatrix
433 operator * (const ComplexColumnVector& a, const ComplexRowVector& b);
434 
435 extern OCTAVE_API ComplexMatrix
436 Givens (const Complex&, const Complex&);
437 
438 extern OCTAVE_API ComplexMatrix
439 Sylvester (const ComplexMatrix&, const ComplexMatrix&, const ComplexMatrix&);
440 
441 extern OCTAVE_API ComplexMatrix
442 xgemm (const ComplexMatrix& a, const ComplexMatrix& b,
443        blas_trans_type transa = blas_no_trans,
444        blas_trans_type transb = blas_no_trans);
445 
446 extern OCTAVE_API ComplexMatrix operator * (const Matrix&,
447                                             const ComplexMatrix&);
448 extern OCTAVE_API ComplexMatrix operator * (const ComplexMatrix&,
449                                             const Matrix&);
450 extern OCTAVE_API ComplexMatrix operator * (const ComplexMatrix&,
451                                             const ComplexMatrix&);
452 
453 extern OCTAVE_API ComplexMatrix min (const Complex& c, const ComplexMatrix& m);
454 extern OCTAVE_API ComplexMatrix min (const ComplexMatrix& m, const Complex& c);
455 extern OCTAVE_API ComplexMatrix min (const ComplexMatrix& a,
456                                      const ComplexMatrix& b);
457 
458 extern OCTAVE_API ComplexMatrix max (const Complex& c, const ComplexMatrix& m);
459 extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& m, const Complex& c);
460 extern OCTAVE_API ComplexMatrix max (const ComplexMatrix& a,
461                                      const ComplexMatrix& b);
462 
463 extern OCTAVE_API ComplexMatrix linspace (const ComplexColumnVector& x1,
464                                           const ComplexColumnVector& x2,
465                                           octave_idx_type n);
466 
467 MS_CMP_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API)
468 MS_BOOL_OP_DECLS (ComplexMatrix, Complex, OCTAVE_API)
469 
470 SM_CMP_OP_DECLS (Complex, ComplexMatrix, OCTAVE_API)
471 SM_BOOL_OP_DECLS (Complex, ComplexMatrix, OCTAVE_API)
472 
473 MM_CMP_OP_DECLS (ComplexMatrix, ComplexMatrix, OCTAVE_API)
474 MM_BOOL_OP_DECLS (ComplexMatrix, ComplexMatrix, OCTAVE_API)
475 
476 MARRAY_FORWARD_DEFS (MArray, ComplexMatrix, Complex)
477 
478 #endif
479